Wednesday, 11 July 2018

rasterization - How to convert file .shp to .tif using ogr or python or gdal


I need to convert files. I'm trying to use gdal_rasterize



gdal_rasterize -b 1 -burn 0 -l file file.shp work.tif

and the output is the following:



ERROR 10: Pointer 'hDS' is NULL in 'GDALGetProjectionRef'. Warning 1: The input vector layer has a SRS, but the output raster dataset SRS is unknown. Ensure output raster dataset has the same SRS, otherwise results might be incorrect. Segmentation fault (core dumped)




Answer



To solve your projection error, you could obtain the projection information from an existing ("reference") image prior to rasterisation. This ensures that the output raster has the same projection, extent and pixel resolution as other images you are currently working with.


Here is an example of how this can be achieved in Python/GDAL. I have tested it on point shapefiles only, but it should work for polygons as well.


# A script to rasterise a shapefile to the same projection & pixel resolution as a reference image.

from osgeo import ogr, gdal
import subprocess

InputVector = 'VectorName.shp'
OutputImage = 'Result.tif'

RefImage = 'Image_Name.tif'

gdalformat = 'GTiff'
datatype = gdal.GDT_Byte

burnVal = 1 #value for the output image pixels
##########################################################
# Get projection info from reference image
Image = gdal.Open(RefImage, gdal.GA_ReadOnly)

# Open Shapefile
Shapefile = ogr.Open(InputVector)
Shapefile_layer = Shapefile.GetLayer()

# Rasterise

print("Rasterising shapefile...")
Output = gdal.GetDriverByName(gdalformat).Create(OutputImage, Image.RasterXSize, Image.RasterYSize, 1, datatype, options=['COMPRESS=DEFLATE'])
Output.SetProjection(Image.GetProjectionRef())
Output.SetGeoTransform(Image.GetGeoTransform())

# Write data to band 1
Band = Output.GetRasterBand(1)
Band.SetNoDataValue(0)
gdal.RasterizeLayer(Output, [1], Shapefile_layer, burn_values=[burnVal])


# Close datasets
Band = None
Output = None
Image = None
Shapefile = None

# Build image overviews
subprocess.call("gdaladdo --config COMPRESS_OVERVIEW DEFLATE "+OutputImage+" 2 4 8 16 32 64", shell=True)
print("Done.")

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...