3

I would like to convert the individual records (polygons/features) of a shapefile to raster using GDAL/OGR. Here is the example shape and the mask raster. Rasterizing the full layer with all records (two, in this example) at once works:

from osgeo import gdal
from osgeo import ogr
from osgeo import gdalconst

raster_file = r"c:\Temp\mask.tif"
shape_file = r"c:\Temp\test2.shp"

# open raster and obtain extent and properties
data = gdal.Open(raster_file, gdalconst.GA_ReadOnly)
geo_transform = data.GetGeoTransform()
x_min = geo_transform[0]
y_max = geo_transform[3]
x_max = x_min + geo_transform[1] * data.RasterXSize
y_min = y_max + geo_transform[5] * data.RasterYSize
x_res = data.RasterXSize
y_res = data.RasterYSize
pixel_width = geo_transform[1]

# open shapefile and create layer
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shape_file, 0)
layer = dataSource.GetLayer()

# convert the layer with field 'presence' to raster
target_ds = gdal.GetDriverByName('GTiff').Create(r"C:\Temp\test2_rasterized.tif", x_res, y_res, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width))
gdal.RasterizeLayer(target_ds, [1], layer, options=["ATTRIBUTE=presence"])

When I loop over the individual features, I can't find a way to export the feature to raster:

from osgeo import gdal
from osgeo import ogr
from osgeo import gdalconst
import os

raster_file = r"c:\Temp\mask.tif"
shape_file = r"c:\Temp\test2.shp"

# open raster and obtain extent and properties
data = gdal.Open(raster_file, gdalconst.GA_ReadOnly)
geo_transform = data.GetGeoTransform()
x_min = geo_transform[0]
y_max = geo_transform[3]
x_max = x_min + geo_transform[1] * data.RasterXSize
y_min = y_max + geo_transform[5] * data.RasterYSize
x_res = data.RasterXSize
y_res = data.RasterYSize
pixel_width = geo_transform[1]

# open shapefile and create layer
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shape_file, 0)
layer = dataSource.GetLayer()
nr_features = layer.GetFeatureCount()

for feat_id in range(nr_features):
    feature = layer.GetFeature(feat_id)
    output = os.path.join(r"C:\Temp", "rasterized"+str(feat_id)+".tif")
    target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Byte)
    target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width))

    # this does not work since I do not pass a layer to the command
    gdal.RasterizeLayer(target_ds, [1], feature, options=["ATTRIBUTE=presence"])

    # e.g. is it possible to achieve it with Rasterize?
    gdal.Rasterize(...) 

Is what I am after at all possible using GDAL/OGR? Or do I need a workaround like converting each feature to an individual layer and then rasterize each layer?

  • Complete workaround for converting each feature to an individual layer and then rasterize each one it can be observed in my answer. – xunilk Aug 2 '19 at 12:15
2

Complete workaround for converting each feature to an individual layer and then rasterize each one it can be observed in following code (where my own paths to layers were used):

from osgeo import gdal
from osgeo import ogr
from osgeo import gdalconst
import os

raster_file = "/home/zeito/pyqgis_data/mask.tif"
shape_file = "/home/zeito/pyqgis_data/test2.shp"

# open raster and obtain extent and properties
data = gdal.Open(raster_file, gdalconst.GA_ReadOnly)
geo_transform = data.GetGeoTransform()
x_min = geo_transform[0]
y_max = geo_transform[3]
x_max = x_min + geo_transform[1] * data.RasterXSize
y_min = y_max + geo_transform[5] * data.RasterYSize
x_res = data.RasterXSize
y_res = data.RasterYSize
pixel_width = geo_transform[1]

# open shapefile and create layer
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shape_file, 0)
layer = dataSource.GetLayer()
nr_features = layer.GetFeatureCount()

for feat_id in range(nr_features):
    feature = layer.GetFeature(feat_id)
    id_no = feature.GetField("id_no")
    presence = feature.GetField("presence")
    outShapefile = os.path.join("/home/zeito/pyqgis_data/", "feature" + str(feat_id)+".shp")
    output = os.path.join("/home/zeito/pyqgis_data/", "rasterized" + str(feat_id)+".tif")

    outDriver = ogr.GetDriverByName("ESRI Shapefile")

    # Remove output shapefile if it already exists
    if os.path.exists(outShapefile):
        outDriver.DeleteDataSource(outShapefile)

    # Create the output shapefile
    outDataSource = outDriver.CreateDataSource(outShapefile)
    outLayer = outDataSource.CreateLayer("feature" + str(feat_id), geom_type=ogr.wkbPolygon)

    # Add an ID field
    idField = ogr.FieldDefn("id_no", ogr.OFTInteger)
    outLayer.CreateField(idField)

    # Add a presence field
    presField = ogr.FieldDefn("presence", ogr.OFTInteger)
    outLayer.CreateField(presField)

    # Create the feature and set values
    featureDefn = outLayer.GetLayerDefn()
    feat = ogr.Feature(featureDefn)
    geom = feature.GetGeometryRef()

    feat.SetGeometry(geom)
    feat.SetField("id_no", id_no)
    feat.SetField("presence", presence)
    outLayer.CreateFeature(feat)
    feat = None

    # Save and close DataSource
    outDataSource = None

    target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Byte)
    target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width))

    new_dataSource = driver.Open(outShapefile, 0)
    new_layer = new_dataSource.GetLayer()

    # this works now 
    gdal.RasterizeLayer(target_ds, [1], new_layer, options=["ATTRIBUTE=presence"])
    target_ds = None

    new_dataSource = None

dataSource = None

After running above code at Python Console of QGIS 3.4, vector and raster layers were obtained as expected.

Initial situation:

enter image description here

Rasterized feature 0:

enter image description here

Rasterized feature 1:

enter image description here

| improve this answer | |
  • Thanks @xunilk - this works! So it seems I indeed need to go the route via the individual layers. At the end I have thousands of records and I would need to delete the shapefiles again (which is ok to do) but do you know if it possible to keep the shapes/layers in-memory? – openwater Aug 2 '19 at 22:05
  • 1
    Yes, it's possible. See section "Load data to memory" in pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html . – xunilk Aug 2 '19 at 23:48
  • Thanks - I adapted your code so that the layer is created in memory. I think this is more suitable for my task since I am only interested in the rasterized output. May I edit your code? – openwater Aug 3 '19 at 16:23
0

Not sure if anyone is struggling with this but I spent a day on this and the solution was simple. I was itereating through features in a shapefile to export images of a specific feature with an empty mask outside each features geometry. Easier than I thought but it took me a while.

    feature = layer.GetFeature(i)
    geom = feature.GetGeometryRef()
    minX, maxX, minY, maxY = geom.GetEnvelope() ###gets extents to create image within

    json_coords = geom.ExportToJson() ###convert geometry to json coordinate

  ###use the json coordinates as culineDSName and cropToCutline as True to create the mask with GdalWarp

     output_options = gdal.WarpOptions(outputBounds=[minX, minY, maxX, maxY], dstAlpha=True,  multithread=True, cutlineDSName = json_coords, cropToCutline=True)
    gdal.Warp(output_file, Raster, options=output_options)
| improve this answer | |
  • What is your variable or import "Raster"? – openwater Nov 22 '19 at 9:35
  • Your "Raster" variable seems to be my initial variable "raster_file". However, when I run your code, the resulting tiffs contain random pixels within the original polygons only and not a full representation of the polygon. The tiffs have two bands where pixels contain either value "0" or "no Data" for band 1 and either value "0" or "255" for band 2. Can you add a screenshot of your result and paste the whole code that you use? – openwater Nov 22 '19 at 10:08
  • My Raster Variable is a large Raster file from which I'm clipping from. As it iterates through the shapefile it will make a small tiff file for each polygon feature with an empty mask outside of each features geometry. Within each feautres geometry it will show whatever was in the original Raster file that I'm clipping form. My final output is a folder full of smaller images of what is within each features geometry. Sorry I can't post my results because I can't share specifically what I'm working on. – Mike Walters Nov 25 '19 at 16:31

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.