How to read geoDataFrame into GDAL function?

Question

I would like to use gdal function in python to convert vector file with polygons to raster (example below). As an input vector data I need to use geoDataFrame (geopandas). Is there a way to read geoDataFrame directly to Gdal function?

gdal.Rasterize(output_raster, vector, xRes=pixel_size,yRes=pixel_size, attribute='waga', outputBounds=[xmin, ymin, xmax, ymax],allTouched=True, outputType=gdal.GDT_Float32)

Answer 1

Sadly there is no direct way to use your GeoDataFrame in gdal's Rasterize function.

In the GDAL/OGR cookbook there's an example on how to perform this by reading a Vector layer from your disk.

Original example from the GDAL/OGR Cookbook

from osgeo import gdal, ogr

# Define pixel_size and NoData value of new raster
pixel_size = 25
NoData_value = -9999

# Filename of input OGR file
vector_fn = 'test.shp'

# Filename of the raster Tiff that will be created
raster_fn = 'test.tif'

# Open the data source and read in the extent
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()

# 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(raster_fn, x_res, y_res, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)

# Rasterize
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[0]))

If you want to rasterize data that you have stored in a GeoPandas GeoDataFrame, you'll have to write it to disk first. Here's a modified version of the code above that assumes you want to rasterize the contents of a GeoDataFrame called source_gdf.

Modified example from the GDAL/OGR Cookbook using GeoPandas GeoDataFrame

# I'm assuming that `source_gdf` has already been populated with geodata.
# The only thing we're doing is exporting the contents of `source_gdf` to disk.
source_gdf_fname = '/path/to/file/output.gpkg'
source_gdf_layername = 'source_gdf'
source_gdf.to_file(source_gdf_fname, layer=source_gdf_layername, driver='GPKG')

from osgeo import gdal, ogr

# Define pixel_size and NoData value of new raster
pixel_size = 25
NoData_value = -9999

# Filename of input OGR file
vector_fn = 'test.shp'

# Filename of the raster Tiff that will be created
raster_fn = 'test.tif'

# Open the data source and read in the extent
# The only change to the code here is that I point to the file 
# and layer we just created at the top of the code. The rest of the
# code snippet is identical to the last. 
source_ds = ogr.Open(source_gdf_fname)
source_layer = source_ds.GetLayer(source_gdf_layername)
x_min, x_max, y_min, y_max = source_layer.GetExtent()

# 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(raster_fn, x_res, y_res, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)

# Rasterize
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[0]))

Answer 2

The standard GDAL python bindings don't interface very easily with geopandas, as you noticed.

There are several alternatives available though. They also use GDAL under the hood, but they have a "more pythonic API" and have better integration with other python libraries like geopandas. One of the more popular ones is rasterio.

This is a code sample that shows how you can rasterize a GeoDataFrame using rasterio.rasterize:

import geopandas as gpd
import rasterio
from rasterio import features
from shapely.geometry import Polygon
import rasterio.plot

# Prepare GeoDataFrame to be rasterized
geometries = [
    Polygon([(0, 5), (5, 5), (5, 0), (0, 5)]),
    Polygon([(10, 10), (10, 15), (15, 10), (10, 10)]),
]
# The waga column will be used as the burn value
data = {"waga": [1.0, 2.0]}
gdf = gpd.GeoDataFrame(data=data, geometry=geometries, crs=31370)

# Prepare some variables
xmin, ymin, xmax, ymax = gdf.total_bounds
pixel_size = 1
width = int((xmax - xmin) // pixel_size)
height = int((ymax - ymin) // pixel_size)
transform = rasterio.transform.from_origin(xmin, ymax, pixel_size, pixel_size)

# Burn geometries
shapes = ((geom, value) for geom, value in zip(gdf.geometry, gdf.waga))
burned = features.rasterize(
    shapes=shapes, out_shape=(width, height), transform=transform, all_touched=True
)

# Write result
output_path = "output.tif"
with rasterio.open(
    output_path,
    mode="w",
    driver="GTiff",
    dtype="float32",
    height=height,
    width=width,
    count=1,
    crs=gdf.crs,
    transform=transform,
    compress="lzw",
) as dest:
    dest.write_band(1, burned)

# Plot result
with rasterio.open(output_path) as src:
    image = src.read(1)
    rasterio.plot.show(image, transform=src.transform)