[gdal-dev] Burn 3d-lines from python

[Even Rouault]

Le mardi 11 août 2015 11:36:32, Simon Lyngby Kokkendorff a écrit :
> Hello list,
> 
>   I'm trying to rasterize some 3d-lines programtically from python. Using
> gdal 1.11.2 from the osgeo4w64 environment. Seems to work fine, except that
> I get a z-offset of 255 (one byte) in my output raster. Works fine if I use
> gdal_rasterize with the -3d option.
> 
>  The following python code should reproduce the problem:
> 
> from osgeo import ogr,gdal
> print gdal.__version__
> georef=[0,1.0,0,10.5,0,-1.0]
> m_drv_ogr=ogr.GetDriverByName("Memory")
> line_ds = m_drv_ogr.CreateDataSource( "dummy")
> layer = line_ds.CreateLayer( "lines", None, ogr.wkbLineString25D)
> layerdefn=layer.GetLayerDefn()
> line=ogr.Geometry(ogr.wkbLineString25D)
> line.AddPoint(0,5.5,0)
> line.AddPoint(10,5.5,10)
> feature=ogr.Feature(layerdefn)
> feature.SetGeometry(line)
> res=layer.CreateFeature(feature)
> layer.ResetReading()
> m_drv_gdal=gdal.GetDriverByName("MEM")
> raster_ds=m_drv_gdal.Create("dummy",10,10,1,gdal.GDT_Float64)
> raster_ds.SetGeoTransform(georef)
> ok=gdal.RasterizeLayer(raster_ds,[1], layer, options=["BURN_VALUE_FROM=Z"])
> A=raster_ds.ReadAsArray()
> print A
> 
> Will output:
> 
> [[   0.    0.    0.    0.    0.    0.    0.    0.    0.    0.]
>  [   0.    0.    0.    0.    0.    0.    0.    0.    0.    0.]
>  [   0.    0.    0.    0.    0.    0.    0.    0.    0.    0.]
>  [   0.    0.    0.    0.    0.    0.    0.    0.    0.    0.]
>  [   0.    0.    0.    0.    0.    0.    0.    0.    0.    0.]
>  [ 255.  256.  257.  258.  259.  260.  261.  262.  263.  264.]
>  [   0.    0.    0.    0.    0.    0.    0.    0.    0.    0.]
>  [   0.    0.    0.    0.    0.    0.    0.    0.    0.    0.]
>  [   0.    0.    0.    0.    0.    0.    0.    0.    0.    0.]
>  [   0.    0.    0.    0.    0.    0.    0.    0.    0.    0.]]

Simon,

looking closer, the issue is that the Python bindings of gdal.RasterizeLayer() 
have a default constant burn value of 255 if not explicitly specified. This 
constant burn value comes in addition to the other source of burning, such as 
Z.
So the fix is to add burn_values = [0] as a keyword of gdal.RasterizeLayer()

I guess the rationale of that was that without any option, the rasterizing 
would burn 255, which is a reasonable default. But it is true that it is 
rather unexpected if you specify BURN_VALUE_FROM=Z

Even


> 
> Cheers,
> Simon Kokkendorff, Dansish Geodata Agency

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com