Rasterising multiple shapefiles into a single raster

Question

I'm not even sure you can do this, but I want to convert multiple shapefiles into a single single-band raster.

I'm pretty sure its going to be easier to first create a multi-band raster and then use each combination of band values and a lookup table of some sort to assign a pixel a single value, producing a single-band raster. But I'm stuck on creating the initial multi-band raster.

Lets say I have 4 shapefiles (shp1, shp2, shp3, shp4).
The attribute table in each shapefile has a shp_V column that has an integer ranging from 5-8.
I want the raster to have 4 bands, one for each shapefile that equals the shp_V value of one of the shapefiles.

All examples I can find with rasterio or gdal only rasterises one shapefile, which I know how to do.

Currently I have this, but it obviously only converts one shapefile to one single-band raster:

attribute = 'shp_V'
shp1 = "Shapefiles/shp1.shp"
shp2 = "Shapefiles/shp2.shp"
shp3 = "Shapefiles/shp3.shp"
shp4 = "Shapefiles/shp4.shp"
out_tif = "Rasters/raster1234.tif"

cmd = 'gdal_rasterize -a {} -co COMPRESS=DEFLATE 
       -co BIGTIFF=YES -tr 0.3 0.3 {} {}'.format(attribute, shp1, out_tif)
print(cmd)
os.system(cmd)

The shapefiles and resulting rasters can get quite large, so any solutions which are a bit easier on processing power would also be a benefit.

snowman2 · Answer

I would recommend a combination of geopandas and geocube.

Here is some untested set of code that should get you pretty close to what you want to do.

Step 1: Combine the shapefiles

import pandas
import geopandas

gpd1 = geopandas.read_file("Shapefiles/shp1.shp")
gpd2 = geopandas.read_file("Shapefiles/shp2.shp")
gpd3 = geopandas.read_file("Shapefiles/shp3.shp")
gpd4 = geopandas.read_file("Shapefiles/shp4.shp")

merged_gpd = pandas.concat([gpd1, gpd2, gpd3, gpd4]).set_geometry("geometry")

Step 3: Rasterize the data

from geocube.api.core import make_geocube

grid = make_geocube(
    vector_data=merged_gpd,
    measurements=["shp_V"],
    resolution=(-0.3, 0.3),
)

Step 4: Export to raster

grid.shp_V.rio.to_raster("Rasters/raster1234.tif", compress="DEFLATE", bigtiff="YES")

Rasterising multiple shapefiles into a single raster

One Answer

Add your own answers!

Ask a Question