Rasterize a shapefile with Geopandas or fiona - python

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.

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")