17

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))
2
  • The answer is about GDALrasterize, I am precisely asking if someone has an idea about doing the same thing using Geopandas and rasterio. Not a duplicate Jun 18 '15 at 15:32
  • Found a piece of Code that might help, post edited Jun 19 '15 at 22:54
33

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.

1
  • What does compress='lzw' do?
    – A Merii
    Nov 20 '20 at 18:10
9

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

enter image description here

Lastly, export it to a raster:

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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