Changes between Version 15 and Version 16 of FAQRaster

Timestamp:

Oct 5, 2009, 3:34:55 PM (15 years ago)

Author:

leaffan

Comment:

included a second alternative for blank raster creation

FAQRaster

@@ -28 +28 @@
 One approach 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 
+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 geo-tiff that matches the shapefile extents. It deals with two conceivable alternatives: a) the calculation of the output pixel resolution of the geo-tiff by using the desired dataset dimensions and b) the calculation of the output dimensions based on desired pixel resolutions:
 
 {{{
@@ -46 +45 @@
 shp = 'test.shp'
 tiff = 'test.tif'
+# Alternative A: Set raster dataset dimensions
 tiff_width = 7850
 tiff_height = 3500
+# Alternative B: Set pixel resolution
+x_res = 0.5
+y_res = 0.75 
 
 # Import vector shapefile
@@ -55 +58 @@
 src_extent = src_lyr.GetExtent()
 
+# Alternative A: Calculate pixel resolution from desired dataset dimensions
+x_res = (src_extent[1] - src_extent[0]) / tiff_width
+y_res = (src_extent[3] - src_extent[2]) / tiff_height
+# Alternative B: Calculate raster dataset dimensions from desired pixel resolution
+tiff_width = int(math.ceil(abs(src_extent[1] - src_extent[0]) / x_res))
+tiff_height = int(math.ceil(abs(src_extent[3] - src_extent[2]) / y_res))
+
+# Create raster GeoTransform based on upper left corner and pixel dimensions
+raster_transform = [src_extent[0], x_res, 0.0, src_extent[3], 0.0, -y_res]
+dst_ds.SetGeoTransform(raster_transform)
+
 # 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
-xres = (src_extent[1] - src_extent[0]) / tiff_width
-yres = (src_extent[3] - src_extent[2]) / tiff_height
-raster_transform = [src_extent[0], xres, 0.0, src_extent[3], 0.0, -yres]
-dst_ds.SetGeoTransform( raster_transform )
+dst_ds = raster.Create(tiff, tiff_width, tiff_height, 4, gdal.GDT_Byte)
 
 # Get projection of shapefile and assigned to raster
 srs = osr.SpatialReference()
 srs.ImportFromWkt(src_lyr.GetSpatialRef().__str__())
-dst_ds.SetProjection( srs.ExportToWkt() )
+dst_ds.SetProjection(srs.ExportToWkt())
 
 # Create blank raster with fully opaque alpha band