[gdal-dev] Rasterize a shapefile as mask

shawn gong shawngong.gdal at gmail.com
Mon Mar 31 08:19:36 PDT 2014


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])
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20140331/7331af0b/attachment.html>


More information about the gdal-dev mailing list