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? Commented Jan 21, 2020 at 3:51
  • Oops yeah. Technically both, I've updated the question so hopefully it makes a bit more sense 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. 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 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 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/… Commented Jan 21, 2020 at 13:33
  • This got me very close! Thanks for the help :) Commented Jan 22, 2020 at 4:08

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.