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.

Improve this question
4
  • 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 ♦
    Commented Jan 21, 2020 at 3:51
  • Oops yeah. Technically both, I've updated the question so hopefully it makes a bit more sense
    – Nebulous29
    Commented Jan 21, 2020 at 4:01
  • 1
    GDAL_Rasterize -b <band index> 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 gdal.org/programs/gdal_merge.html to 'stack' the bands into a single raster or use driver.Create 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
    Commented Jan 21, 2020 at 4:12
  • 2
    I would rasterize each shp then stack. This post may be helpful for the stacking: gis.stackexchange.com/q/223910/8104
    – Aaron ♦
    Commented Jan 21, 2020 at 4:37

1 Answer 1

Reset to default
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")
Improve this answer
3
  • This looks like a really clean solution, but the pd.concat causes a MemoryError due to the file sizes
    – Nebulous29
    Commented Jan 21, 2020 at 5:42
  • If you have many columns in the shapfile and only need one, you could try this to reduce memory: gis.stackexchange.com/questions/129414/…
    – snowman2
    Commented Jan 21, 2020 at 13:33
  • This got me very close! Thanks for the help :)
    – Nebulous29
    Commented Jan 22, 2020 at 4:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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