Changes between Version 10 and Version 11 of FAQRaster

Timestamp:

Jul 31, 2008, 12:07:53 PM (16 years ago)

Author:

springmeyer

Comment:

Add python example for creating a blank GTiff with numpy from shapefile extents

FAQRaster

@@ -20 +20 @@
 
 ''TBD''
+
+== How can I create a blank raster based on a vector files extents for use with gdal_rasterize? ==
+
+Currently the gdal_rasterize utility cannot create a blank raster based a vector datasource although this feature has been requested (ticket 1599).
+Until that feature is add you'll need to find a way to generate a raster that matches the extent of your vector layer that you plan to use to burn features into the raster.
+
+One way to do this is to use gdal_translate to clip out and flatten an existing raster file that covers your area of interest. See this [http://lists.osgeo.org/pipermail/gdal-dev/2008-February/016061.html example of syntax].
+
+Another approach is to use one of the language bindings to GDAL to create a blank raster from scratch. The tricky part is understanding the !GeoTransform syntax.
+This python snippet shows how to read in a shapefile and output a tiff that matches the shapefile extents 
+ * Note: you'll need to modify the `px` variable based on your desired resolution and raster dimensions
+
+{{{
+#!/usr/bin/env python
+
+from osgeo import gdal
+from osgeo import osr
+from osgeo import ogr
+import numpy
+
+shp = 'test.shp'
+tiff = 'test.tif'
+px = .001
+tiff_width = 7850
+tiff_height = 3500
+
+# Import vector shapefile
+vector = ogr.GetDriverByName('ESRI Shapefile')
+src_ds = vector.Open(shp)
+src_lyr = src_ds.GetLayerByIndex(index=0)
+src_extent = src_lyr.GetExtent()
+
+# Create new raster layer with 4 bands
+raster = gdal.GetDriverByName('GTiff')
+dst_ds = raster.Create( tiff, tiff_width, tiff_height, 4, gdal.GDT_Byte)
+
+# Create raster GeoTransform based on upper left corner and pixel resolution
+raster_transform = [src_extent[0], px, 0.0, src_extent[3], 0.0, -px]
+dst_ds.SetGeoTransform( raster_transform )
+
+# Get projection of shapefile and assigned to raster
+srs = osr.SpatialReference()
+srs.ImportFromWkt(src_lyr.GetSpatialRef().__str__())
+dst_ds.SetProjection( srs.ExportToWkt() )
+
+# Create blank raster with fully opaque alpha band
+zeros = numpy.zeros( (tiff_height, tiff_width), numpy.uint8 )
+dst_ds.GetRasterBand(1).WriteArray( zeros )
+dst_ds.GetRasterBand(2).WriteArray( zeros )
+dst_ds.GetRasterBand(3).WriteArray( zeros )
+opaque = numpy.ones((tiff_height,tiff_width), numpy.uint8 )*255
+dst_ds.GetRasterBand(4).WriteArray( opaque )
+
+}}}
+
 
 == Can I use gdal_rasterize to generate non-solid polygons? ==