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.
You can use the rasterstats package for this.
For example (assuming you have a geopandas.GeoDataFrame called gdf):
from rasterstats import zonal_stats
with rasterio.open("/path/to/raster.tif") as src:
affine = src.transform
array = src.read(1)
df_zonal_stats = pd.DataFrame(zonal_stats(gdf, array, affine=affine))
# adding statistics back to original GeoDataFrame
gdf2 = pd.concat([gdf, df_zonal_stats], axis=1)
Best Answer
You can do zonal statistics from a GeoDataFrame directly on a GeoTiff using rasterstats.
There are some good examples of rasterstats integration on the wiki