I have an existing arcpy script which uses PolygonToRaster_conversion to convert a polygon shapefile to raster (tiff). I am now exploring GDAL/OGR and have used the cookbook approach to carry out the same task but find two key differences in the output rasters:
NoDataVal = -9999 # Open the data source and read in the extent shpDS = ogr.Open(inPolygonShp) shpLayer = shpDS.GetLayer() # Create the destination data source xMin, xMax, yMin, yMax = shpLayer.GetExtent() xRes = int((xMax - xMin) / inGridSize) yRes = int((yMax - yMin) / inGridSize) rasterDS = gdal.GetDriverByName('GTiff').Create(outputRaster, xRes, yRes, 1, gdal.GDT_Byte) # Define spatial reference rasterDS.SetProjection(shpLayer.GetSpatialRef().ExportToWkt()) rasterDS.SetGeoTransform((xMin, inGridSize, 0, yMax, 0, -inGridSize)) rBand = rasterDS.GetRasterBand(1) rBand.SetNoDataValue(NoDataVal) rBand.Fill(NoDataVal) # Rasterize err = gdal.RasterizeLayer(rasterDS, [1], shpLayer, burn_values=[1], options = ["ALL_TOUCHED=TRUE"]) # http://www.gdal.org/gdal_rasterize.html and using arcpy:
idFld = arcpy.Describe(inPolygonShp).OIDFieldName addMsgAndPrint("OID Field used: "+ idFld) arcpy.PolygonToRaster_conversion(inPolygonShp, idFld, outputRaster, "CELL_CENTER", "NONE", inGridSize)Any advice appreciated.Thanks
أكثر...
- The extent of the GDAL output raster is different from the arcpy output. Xmin is the same in each. However, the number of rows and cols are one less each in the gdal rasters, and the Y values differ by an odd amount (which is not a multiple of grid pixel size). I have tried the rasterize function with and without the ALL_TOUCHED option to no avail. I suspectthe difference is in calculating the extent, and applying the GeoTransform. Any clues?
- The NoData value set using gdal is neither correctly set or recognised in ArcGIS or QGIS. Instead of a NoData value, the pixels have a value of zero. Could this be to do with the output raster datatype?
NoDataVal = -9999 # Open the data source and read in the extent shpDS = ogr.Open(inPolygonShp) shpLayer = shpDS.GetLayer() # Create the destination data source xMin, xMax, yMin, yMax = shpLayer.GetExtent() xRes = int((xMax - xMin) / inGridSize) yRes = int((yMax - yMin) / inGridSize) rasterDS = gdal.GetDriverByName('GTiff').Create(outputRaster, xRes, yRes, 1, gdal.GDT_Byte) # Define spatial reference rasterDS.SetProjection(shpLayer.GetSpatialRef().ExportToWkt()) rasterDS.SetGeoTransform((xMin, inGridSize, 0, yMax, 0, -inGridSize)) rBand = rasterDS.GetRasterBand(1) rBand.SetNoDataValue(NoDataVal) rBand.Fill(NoDataVal) # Rasterize err = gdal.RasterizeLayer(rasterDS, [1], shpLayer, burn_values=[1], options = ["ALL_TOUCHED=TRUE"]) # http://www.gdal.org/gdal_rasterize.html and using arcpy:
idFld = arcpy.Describe(inPolygonShp).OIDFieldName addMsgAndPrint("OID Field used: "+ idFld) arcpy.PolygonToRaster_conversion(inPolygonShp, idFld, outputRaster, "CELL_CENTER", "NONE", inGridSize)Any advice appreciated.Thanks
أكثر...