4

I'm writing some code that's supposed to take multiple multipolygon layers and merge them into 1 raster. I can load in the datasets just fine and render them to an image (although I use this just for debugging), here's what I'm doing right now:

noiselayers["noise"] = [[1, 0.1], [2, 0.8]]

    #create a source in memory
    ds = ogr.GetDriverByName("Memory").CreateDataSource("wrk")
    source = ds.CreateLayer("", None, ogr.wkbMultiPolygon )
    source.CreateField(ogr.FieldDefn("VAL", ogr.OFSTFloat32))

    for layerdata in noiselayers[selected]:
        data = requests.get("http://localhost/api/layer/export/%s" % str(layerdata[0]))
        geometry = json.loads(data.text)

        for g in geometry:
            geom = ogr.CreateGeometryFromWkt(g['the_geom'])

            structure = ogr.Feature(source.GetLayerDefn())
            structure.SetGeometry(geom)
            structure.SetField("VAL", layerdata[1])

            #save geometry to layer
            source.CreateFeature(structure)


    pixelsize = 1000
    #novalue = -9999

    #get the layer bounds dynamically from the layer        
    #x_min, x_max, y_min, y_max = layer.GetExtent()

    #extracted values from ns_play_area_final_EPSG_3035
    x_min = 2777236.78057
    x_max = 5065039.36824
    y_min = 2815620.0976
    y_max = 4485652.62185

    x_res = int((x_max - x_min) / pixelsize)
    y_res = int((y_max - y_min) / pixelsize)

    target = gdal.GetDriverByName("MEM").Create("", x_res, y_res, 1, gdal.GDT_Byte)
    target.SetGeoTransform((x_min, pixelsize, 0, y_max, 0, -pixelsize))            

    gdal.RasterizeLayer(target, [1], source)

This grabs a bunch of data from our API and rasterizes it just fine. However, I need to merge multiple layers with different weights (in the noiselayers variable, the 2nd value of each list is the weight). When it's a higher value it should increase the value more, for now just setting the value of the square to that value is enough. I've tried using the options (as seen here: https://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py), but they never seem to do anything, here's some of the things I tried to use the options:

gdal.RasterizeLayer(target, [1], source, None, None, [1000], ["ATTRIBUTE=VAL"])

gdal.RasterizeLayer(target, [1], source, options = ["ATTRIBUTE=VAL"])
  • any combination of the 2, the options never seem to do anything other than break the rasterizing. I've also tried inversing (using INVERSE=TRUE, or any combination of capitalizing/not capitalizing that) the raster, but that also never seemed to work.

Does anyone have a working example of how to use attributes and/or other hints on how to solve my issue?

2

Those options never worked for me either. I ended up using a work around by setting attribute filters on the layer object. How about this:

for index, weight in noiselayers[selected]:
    source.SetAttributeFilter('VAL={}'.format(weight))
    gdal.RasterizeLayer(target, [1], source, burn_values=[weight])

# reset the filter
source.SetAttributeFilter('')

you can find the attribute filter documentation here.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.