Search code examples
pythongdalosgeo

Using Python to Convert a Shapefile to Geotiff

I have a .shp file and I would like to convert it into a GEOTIFF. My shape file consists of a large polygon with many polygons inside. I am using the following code, but the output TIF consists of only the large polygon. 

from osgeo import ogr, gdal
import subprocess

InputVector = Shapefile
OutputImage = OutputfileName

gdalformat = 'GTiff'
datatype = gdal.GDT_Byte
burnVal = 1 
Shapefile = ogr.Open(InputVector)
Shapefile_layer = Shapefile.GetLayer()


Output = gdal.GetDriverByName(gdalformat).Create(OutputImage, RasterXSize, RasterYSize, 1, datatype, options=['COMPRESS=DEFLATE'])
Output.SetProjection(Projection)
Output.SetGeoTransform(GeoTransform) 

Band = Output.GetRasterBand(1)
Band.SetNoDataValue(0)
gdal.RasterizeLayer(Output, [1], Shapefile_layer, burn_values=[burnVal])

subprocess.call("gdaladdo --config COMPRESS_OVERVIEW DEFLATE "+OutputImage+" 2 4 8 16 32 64", shell=True)

I am unsure what I am doing wrong here.

Thanks

Solution

  • My first guess would be to make sure that the larger polygon has "holes"/ NoData where the smaller polygons occur. Ill try and be more clear with an example from ArcMap's erase tool. What you want your large polygon to look like is that with the blue tick (with the erase feature obviously overlaying), what I think your large polygon currently looks like is one with the red tick, which will result in your raster only consisting of the large polygon. My apologies if I am totally wrong, then we'll find another solution.

    enter image description here

    If this is the case (your large polygon is similar to the blue tick), burning all the polygons to 1 might also create confusion. You can create a new field in your shapefile (e.g. UID) and give it an unique numerical ID (it has to be numerical). The you can rasterize based on the new field, for example:

    gdal.RasterizeLayer(Output, [1], Shapefile_layer, options = ['ATTRIBUTE=UID'])