[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