2

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
  • 7,321
  • 12
  • 29
  • 54
Nebulous29
  • 126
  • 9
  • You mention "I want to convert multiple shapefiles into a single single-band raster" and then later state "I want the raster to have 4 bands, one for each shapefile that equals the shp_V value of one of the shapefiles." Could you please clarify if you want a single band or multiband raster? – Aaron Jan 21 '20 at 03:51
  • Oops yeah. Technically both, I've updated the question so hopefully it makes a bit more sense – Nebulous29 Jan 21 '20 at 04:01
  • 1
    GDAL_Rasterize -b https://gdal.org/programs/gdal_rasterize.html Not used when creating a new raster. You could consider rasterizing to single band rasters then use GDAL_Merge https://gdal.org/programs/gdal_merge.html to 'stack' the bands into a single raster or use driver.Create https://gdal.org/tutorials/raster_api_tut.html (near the bottom) to make a multiband raster and then rasterize to each band in turn. – Michael Stimson Jan 21 '20 at 04:12
  • 2
    I would rasterize each shp then stack. This post may be helpful for the stacking: https://gis.stackexchange.com/q/223910/8104 – Aaron Jan 21 '20 at 04:37

1 Answers1

3

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")
snowman2
  • 7,321
  • 12
  • 29
  • 54
  • This looks like a really clean solution, but the pd.concat causes a MemoryError due to the file sizes – Nebulous29 Jan 21 '20 at 05:42
  • If you have many columns in the shapfile and only need one, you could try this to reduce memory: https://gis.stackexchange.com/questions/129414/only-read-specific-attribute-columns-of-a-shapefile-with-geopandas-fiona#129422 – snowman2 Jan 21 '20 at 13:33
  • This got me very close! Thanks for the help :) – Nebulous29 Jan 22 '20 at 04:08