[gdal-dev] Rasterize a shapefile as mask

[Kyle Shannon]

On Mon, Mar 31, 2014 at 9:19 AM, shawn gong <shawngong.gdal at gmail.com> wrote:
> Hi list,
>
> I have a RSAT-2 image and a shoreline vector (consists of 2 polygons). I
> want to generate a mask with the same # of lines and pixels as the RSAT-2
> image, and make inside polygons as 0, outside as 1.
>
> The example that I could find is UTM based, which I tried and sort of
> worked. But the mask image is not oblique as the RSAT-2 image. When I used
> the lat/lon projection, the gdal.RasterizeLayer() doesn't work. Would
> someone tell me why?
>
> thanks,
> Shawn
>
>
> UTM version that worked:
>     from osgeo import gdal, ogr, gdal_array
>
>     # Filename of input OGR file
>     vector_file = 'shoreline_utm.shp'
>
>     # Open the data source and read in the extent
>     vector_ds = ogr.Open(vector_file)
>     vector_layer = vector_ds.GetLayer()
>     #vector_srs = vector_layer.GetSpatialRef()
>     x_min, x_max, y_min, y_max = vector_layer.GetExtent()
>
>     # raster Tiff that will be created
>     new_raster = 'shoreline_utm.tif'
>
>     pixel_spacing = 9.77
>     line_spacing = 5.5
>     xsize = int((x_max - x_min) / pixel_spacing)
>     ysize = int((y_max - y_min) / line_spacing)
>
>     new_raster_ds = gdal.GetDriverByName('GTiff').Create(new_raster, xsize,
> ysize, 1, gdal.GDT_Byte)
>     new_raster_ds.SetGeoTransform((x_min, pixel_spacing, 0, y_max, 0,
> -line_spacing))
>
>     band = new_raster_ds.GetRasterBand(1)
>     band.Fill(1)
>     gdal.RasterizeLayer(new_raster_ds, [1], vector_layer, burn_values=[0])
>
>
> lat/lon version:
>     from osgeo import gdal, ogr, gdal_array
>
>     # Filename of input OGR file
>     vector_file = 'shoreline.shp'
>
>     # Open the data source and read in the extent
>     vector_ds = ogr.Open(vector_file)
>     vector_layer = vector_ds.GetLayer()
>
>     # raster Tiff that will be created
>     new_raster = 'shorelineTest.tif'
>
>     # master Tiff whose size will be copied
>     master = 'product.xml'
>     master_ds = gdal.Open(master)
>     xsize, ysize = master_ds.RasterXSize, master_ds.RasterYSize
>
>     new_raster_ds = gdal.GetDriverByName('GTiff').Create(new_raster, xsize,
> ysize, 1, gdal.GDT_Byte)
>
>     # use georeferencing from the main product.xml dataset
>     gdal_array.CopyDatasetInfo(master_ds, new_raster_ds)
>
>     band = new_raster_ds.GetRasterBand(1)
>     band.Fill(1)
>
>     gdal.RasterizeLayer(new_raster_ds, [1], vector_layer, burn_values=[0])
>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev

Shawn,
According to the docs, the layer has to be in the same srs as the
raster.  Are you doing this?  If it's not, you can provide a
transformer, or reproject the geometries yourself.  See:

http://gdal.org/gdal__alg_8h.html#adfe5e5d287d6c184aab03acbfa567cb1

notably the pfnTransformer arg.

-- 
Kyle