On Mon, Mar 31, 2014 at 9:19 AM, shawn gong <shawngong.g...@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@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 _______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev