Rasterizing a GDAL layer
EDIT: I guess I'd use qGIS python bindings: http://www.qgis.org/wiki/Python_Bindings
That's the easiest way I can think of. I remember hand rolling something before, but it's ugly. qGIS would be easier, even if you had to make a separate Windows installation (to get python to work with it) then set up an XML-RPC server to run it in a separate python process.
I you can get GDAL to rasterize properly that's great too.
I haven't used gdal for a while, but here's my guess:
burn_values
is for false color if you don't use Z-values. Everything inside your polygon is [255,0,0]
(red) if you use burn=[1,2,3],burn_values=[255,0,0]
. I'm not sure what happens to points - they might not plot.
Use gdal.RasterizeLayer(ds,bands,layer,burn_values, options = ["BURN_VALUE_FROM=Z"])
if you want to use the Z values.
I'm just pulling this from the tests you were looking at: http://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py
Another approach - pull the polygon objects out, and draw them using shapely, which may not be attractive. Or look into geodjango (I think it uses openlayers to plot into browsers using JavaScript).
Also, do you need to rasterize? A pdf export might be better, if you really want precision.
Actually, I think I found using Matplotlib (after extracting and projecting the features) was easier than rasterization, and I could get a lot more control.
EDIT:
A lower level approach is here:
http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/gdal2grd.py\
Finally, you can iterate over the polygons (after transforming them into a local projection), and plot them directly. But you better not have complex polygons, or you will have a bit of grief. If you have complex polygons ... you are probably best off using shapely and r-tree from http://trac.gispython.org/lab if you want to roll your own plotter.
Geodjango might be a good place to ask .. they will know a lot more than me. Do they have a mailing list? There's also lots of python mapping experts around, but none of them seem to worry about this. I guess they just plot it in qGIS or GRASS or something.
Seriously, I hope that somebody who knows what they are doing can reply.
Yogesh Gandhi
Updated on November 21, 2020Comments
-
Yogesh Gandhi over 3 years
Edit
Here is the proper way to do it, and the documentation:
import random from osgeo import gdal, ogr RASTERIZE_COLOR_FIELD = "__color__" def rasterize(pixel_size=25): # Open the data source orig_data_source = ogr.Open("test.shp") # Make a copy of the layer's data source because we'll need to # modify its attributes table source_ds = ogr.GetDriverByName("Memory").CopyDataSource( orig_data_source, "") source_layer = source_ds.GetLayer(0) source_srs = source_layer.GetSpatialRef() x_min, x_max, y_min, y_max = source_layer.GetExtent() # Create a field in the source layer to hold the features colors field_def = ogr.FieldDefn(RASTERIZE_COLOR_FIELD, ogr.OFTReal) source_layer.CreateField(field_def) source_layer_def = source_layer.GetLayerDefn() field_index = source_layer_def.GetFieldIndex(RASTERIZE_COLOR_FIELD) # Generate random values for the color field (it's here that the value # of the attribute should be used, but you get the idea) for feature in source_layer: feature.SetField(field_index, random.randint(0, 255)) source_layer.SetFeature(feature) # Create the destination data source x_res = int((x_max - x_min) / pixel_size) y_res = int((y_max - y_min) / pixel_size) target_ds = gdal.GetDriverByName('GTiff').Create('test.tif', x_res, y_res, 3, gdal.GDT_Byte) target_ds.SetGeoTransform(( x_min, pixel_size, 0, y_max, 0, -pixel_size, )) if source_srs: # Make the target raster have the same projection as the source target_ds.SetProjection(source_srs.ExportToWkt()) else: # Source has no projection (needs GDAL >= 1.7.0 to work) target_ds.SetProjection('LOCAL_CS["arbitrary"]') # Rasterize err = gdal.RasterizeLayer(target_ds, (3, 2, 1), source_layer, burn_values=(0, 0, 0), options=["ATTRIBUTE=%s" % RASTERIZE_COLOR_FIELD]) if err != 0: raise Exception("error rasterizing layer: %s" % err)
Original question
I'm looking for information on how to use
osgeo.gdal.RasterizeLayer()
(the docstring is very succinct, and I can't find it in the C or C++ API docs. I only found a doc for the java bindings).I adapted a unit test and tried it on a .shp made of polygons:
import os import sys from osgeo import gdal, gdalconst, ogr, osr def rasterize(): # Create a raster to rasterize into. target_ds = gdal.GetDriverByName('GTiff').Create('test.tif', 1280, 1024, 3, gdal.GDT_Byte) # Create a layer to rasterize from. cutline_ds = ogr.Open("data.shp") # Run the algorithm. err = gdal.RasterizeLayer(target_ds, [3,2,1], cutline_ds.GetLayer(0), burn_values=[200,220,240]) if err != 0: print("error:", err) if __name__ == '__main__': rasterize()
It runs fine, but all I obtain is a black .tif.
What's the
burn_values
parameter for ? CanRasterizeLayer()
be used to rasterize a layer with features colored differently based on the value of an attribute ?If it can't, what should I use ? Is AGG suitable for rendering geographic data (I want no antialiasing and a very robust renderer, able to draw very large and very small features correctly, possibly from "dirty data" (degenerate polygons, etc...), and sometimes specified in large coordinates) ?
Here, the polygons are differentiated by the value of an attribute (the colors don't matter, I just want to have a different one for each value of the attribute).
-
Yogesh Gandhi over 14 yearsYes, I need to rasterize (I edited the question to show what I want to do). Maybe there there is an option like "BURN_VALUE_FROM_Z" that could use an attribute ?
-
Yogesh Gandhi over 14 yearsAlso, I don't understand why I end up with a black image in my test.
-
Yogesh Gandhi over 14 yearsI have been able to use the RasterizeLayer with the help of the folks on #gdal. The problem was in the transformation/extent shifting, you have to make the source and destination extents match or everything is clipped out. This is actually explained in the docs, I don't know how I missed it :D Thanks anyway for your research, I'll accept your answer and edit my question with the fixed code.