From: Frank W. <war...@gd...> - 2002-01-18 05:55:25
|
On Fri, Jan 18, 2002 at 11:15:13AM +1000, Denham Robert wrote: > Openev discussants, > > I was wondering about ways I could approach converting vector coverages, > in particular polygon coverages into raster layers. Is this possible in > Openev or in other free systems ? > > Regards > Robert Denham Robert, The polygon to raster support isn't currently accessable from the OpenEV GUI, but it can be done from a script. I wrote the following script to try and demonstrate how the polygon to raster process works. This is actually a freestanding script. It uses OpenEV methods, but doesn't need to be run from within a gui session. I did find I needed to add a get_extents() method on the GvShapes (contain object for holding a files worth shapes such as polygons). If you can't easily update OpenEV from CVS, you will have to comment out the "extents = polygons.get_extents()" call, and replace it with something like "extents = (xmin, ymin, xsize, ysize)" for your layer. Let me know if you have any questions. Best regards, -- ---------------------------------------+-------------------------------------- I set the clouds in motion - turn up | Frank Warmerdam, war...@po... light and sound - activate the windows | http://pobox.com/~warmerdam and watch the world go round - Rush | Geospatial Programmer for Rent #!/usr/bin/env python import sys import gview import gvalg import gdal import Numeric # ---------------------------------------------------------------------------- # Get input/output arguments. # ---------------------------------------------------------------------------- if len(sys.argv) > 2: raster_filename = sys.argv[2] else: raster_filename = 'raster_out.tif' if len(sys.argv) > 1: poly_filename = sys.argv[1] else: poly_filename = '/u/data/esri/shape/eg_data/polygon.shp' # ---------------------------------------------------------------------------- # Read polygons. # ---------------------------------------------------------------------------- polygons = gview.GvShapes( shapefilename = poly_filename ) # ---------------------------------------------------------------------------- # Get extents of polygons, and establish appropriate bounds for raster file. # # (NOTE: the get_extents() method on the GvShapes contain is a new binding # in Python I added Jan 17, 2001) # ---------------------------------------------------------------------------- extents = polygons.get_extents() if extents[2] > extents[3]: xsize = 1024 ysize = int((1024 * extents[3])/extents[2]) else: ysize = 1024 xsize = int((1024 * extents[2])/extents[3]) pixel_size = extents[2] / xsize geotransform = ( extents[0], pixel_size, 0.0, extents[1] + pixel_size * ysize, 0.0, -pixel_size ) # ---------------------------------------------------------------------------- # Create raster file, and set georeferencing transform. # ---------------------------------------------------------------------------- rfile = gdal.GetDriverByName('GTiff').Create(raster_filename,xsize,ysize,1) rfile.SetGeoTransform( geotransform ) raster = gview.GvRaster( dataset = rfile ) rfile = None # ---------------------------------------------------------------------------- # Rasterize # ---------------------------------------------------------------------------- gvalg.rasterize_shapes( raster, polygons, 255 ) raster = None |