Rasterize a shapefile with Geopandas or fiona - python

Question

I need to rasterize a really simple shapefile a bit like this http://tinyurl.com/odfbanu. Which is just a shapefile containing counties in the US. I have seen this previous answer: GDAL RasterizeLayer does not Burn All Polygons to Raster? but I was wondering if there is a way to do it using Geopandas or fiona and maybe rasterio for the tiff writing part.

So my goal is to rasterize it and assign a value to every polygon sharing a common value, LSAD in the example.

So I wrote the beginning of the code inspired by shongololo in the thread: Dissolving polygons based on attributes with Python (shapely, fiona)?.

from geopandas import GeoDataFrame

name_in = 'cb_2013_us_county_20m.shp'

#Open the file with geopandas
counties = GeoDataFrame.from_file(name_in)

#Add a column to the Geodataframe containing the new value
for i in range (len(counties)):
    LSAD = counties.at[i,'LSAD']
    if LSAD == 00 :
        counties['LSAD_NUM'] == 'A'
    elif LSAD == 03 :
        counties['LSAD_NUM'] == 'B'
    elif LSAD == 04 :
        counties['LSAD_NUM'] == 'C'
    elif LSAD == 05 :
        counties['LSAD_NUM'] == 'D'
    elif LSAD == 06 :
        counties['LSAD_NUM'] == 'E'
    elif LSAD == 13 :
        counties['LSAD_NUM'] == 'F'
    elif LSAD == 15 :
        counties['LSAD_NUM'] == 'G'  
    elif LSAD == 25 :
        counties['LSAD_NUM'] == 'I'          
    else :
        counties['LSAD_NUM'] == 'NA'

Really easy stuff, so now I am wondering how I can actually write those shapes to a tiff. I began working with Geopandas as I believed that was the best option but if you have a fiona suggestion I am up for it too.

I found a piece of code from rasterio which seems to be able to take a shapely geometry and burn it into a new raster http://tinyurl.com/op49uek

# I guess that my goal should be to load my list of geometries under geometry to be able to pass it to rasterio later on
geometry = {'type':'Polygon','coordinates':[[(2,2),(2,4.25),(4.25,4.25),(4.25,2),(2,2)]]}
    
with rasterio.drivers():
    result = rasterize([geometry], out_shape=(rows, cols))
    with rasterio.open(
            "test.tif",
            'w',
            driver='GTiff',
            width=cols,
            height=rows,
            count=1,
            dtype=numpy.uint8,
            nodata=0,
            transform=IDENTITY,
            crs={'init': "EPSG:4326"}) as out:
                 out.write_band(1, result.astype(numpy.uint8))

Answer 1

You are on the right track and the geopandas GeoDataFrame is a good choice for rasterization over Fiona. Fiona is a great toolset, but I think that the DataFrame is better suited to shapefiles and geometries than nested dictionaries.

import geopandas as gpd
import rasterio
from rasterio import features

Set up your filenames

shp_fn = 'cb_2013_us_county_20m.shp'
rst_fn = 'template_raster.tif'
out_fn = './rasterized.tif'

Open the file with GeoPANDAS read_file

counties = gpd.read_file(shp_fn)

Add the new column (as in your above code)

for i in range (len(counties)):
    LSAD = counties.at[i,'LSAD']
    if LSAD == 00 :
        counties['LSAD_NUM'] == 'A'
    elif LSAD == 03 :
        counties['LSAD_NUM'] == 'B'
    elif LSAD == 04 :
        counties['LSAD_NUM'] == 'C'
    elif LSAD == 05 :
        counties['LSAD_NUM'] == 'D'
    elif LSAD == 06 :
        counties['LSAD_NUM'] == 'E'
    elif LSAD == 13 :
        counties['LSAD_NUM'] == 'F'
    elif LSAD == 15 :
        counties['LSAD_NUM'] == 'G'  
    elif LSAD == 25 :
        counties['LSAD_NUM'] == 'I'          
    else :
        counties['LSAD_NUM'] == 'NA'

Open the raster file you want to use as a template for feature burning using rasterio

rst = rasterio.open(rst_fn)

copy and update the metadata from the input raster for the output

meta = rst.meta.copy()
meta.update(compress='lzw')

Now burn the features into the raster and write it out

with rasterio.open(out_fn, 'w+', **meta) as out:
    out_arr = out.read(1)

    # this is where we create a generator of geom, value pairs to use in rasterizing
    shapes = ((geom,value) for geom, value in zip(counties.geometry, counties.LSAD_NUM))

    burned = features.rasterize(shapes=shapes, fill=0, out=out_arr, transform=out.transform)
    out.write_band(1, burned)

The overall idea is to create an iterable containing tuples of (geometry, value), where the geometry is a shapely geometry and the value is what you want to burn into the raster at that geometry's location. Both Fiona and GeoPANDAS use shapely geometries so you are in luck there. In this example a generator is used to iterate through the (geometry,value) pairs which were extracted from the GeoDataFrame and joined together using zip().

Make sure you open the out_fn file in w+ mode, because it will have to be used for reading and writing.

Answer 2

geocube is a new tool specifically designed for rasterizing geopandas data that wraps rasterio. It simplifies the process and eliminates the need for a template raster.

https://github.com/corteva/geocube

In the context of the example above:

from geocube.api.core import make_geocube
import geopandas

counties = geopandas.read_file("zip://cb_2013_us_county_20m.zip/cb_2013_us_county_20m.shp")

The letter can be set on the dataframe like so:

counties["LSAD_LETTER"] = 'NA'
lsad_letter = counties.LSAD_LETTER.copy()
lsad_letter[counties.LSAD=='00'] = 'A'
lsad_letter[counties.LSAD=='03'] = 'B'
lsad_letter[counties.LSAD=='04'] = 'C'
lsad_letter[counties.LSAD=='05'] = 'D'
lsad_letter[counties.LSAD=='06'] = 'E'
lsad_letter[counties.LSAD=='13'] = 'F'
lsad_letter[counties.LSAD=='15'] = 'G'
lsad_letter[counties.LSAD=='25'] = 'I'
counties["LSAD_LETTER"] = lsad_letter

However, only numeric values can be rasterized. Here is a categorical example: https://corteva.github.io/geocube/stable/examples/categorical.html

So, instead of using that, use the numbers in string format and convert to integers:

counties["LSAD_NUM"] = counties.LSAD.astype(int)

Then, rasterize the data:

cube = make_geocube(
    counties,
    measurements=["LSAD_NUM"],
    resolution=(1, -1),
)

Lastly, export it to a raster:

cube.LSAD_NUM.rio.to_raster("lsad_num.tif")

Answer 3

Just to complement the answers already provided ...

If you do not have a template raster you can generate the metadata from your shapefile.

import geopandas as gpd
from rasterio import transform

shp = gpd.read_file(shp_name)
bbox = shp.total_bounds
xmin, ymin, xmax, ymax = bbox
res = 100 # desired resolution
w = (xmax - xmin) // res 
h = (ymax - ymin) // res

out_meta = {
    "driver": "GTiff",
    "dtype": "uint8",
    "height": h,
    "width": w,
    "count": 1,
    "crs": shp.crs,
    "transform": transform.from_bounds(xmin, ymin, xmax, ymax, w, h),
    "compress": 'lzw'
}