Rasterizar un shapefile con Geopandas o fiona - python

Question

Necesito rasterizar un shapefile muy simple, un poco como esto http://tinyurl.com/odfbanu . Que no es más que un shapefile que contiene condados de EE.UU.. He visto esta respuesta anterior: ¿GDAL RasterizeLayer no graba todos los polígonos en raster? pero me preguntaba si hay una manera de hacerlo usando Geopandas o fiona y tal vez rasterio para la parte de escritura tiff.

Así que mi objetivo es rasterizarlo y asignar un valor a cada polígono que comparta un valor común, LSAD en el ejemplo.

Así que escribí el principio del código inspirado por shongololo en el hilo: ¿Disolver polígonos basándose en atributos con 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'

Realmente fácil, así que ahora me pregunto cómo puedo escribir esas formas en un tiff. Empecé a trabajar con Geopandas como yo creía que era la mejor opción, pero si usted tiene una sugerencia fiona estoy dispuesto a ello también.

He encontrado un trozo de código de rasterio que parece ser capaz de tomar una geometría con forma y grabarlo en un nuevo 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

Vas por buen camino y el GeoDataFrame de geopandas es una buena opción para la rasterización sobre Fiona. Fiona es un gran conjunto de herramientas, pero creo que el DataFrame es más adecuado para shapefiles y geometrías que para diccionarios anidados.

import geopandas as gpd
import rasterio
from rasterio import features

Configure sus nombres de archivo

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

Abrir el fichero con GeoPANDAS read_file

counties = gpd.read_file(shp_fn)

Añade la nueva columna (como en el código anterior)

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'

Abra el archivo raster que desea utilizar como plantilla para la grabación de características mediante rasterio

rst = rasterio.open(rst_fn)

copiar y actualizar los metadatos de la trama de entrada para la salida

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

Ahora grabar las características en el raster y escribirlo

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)

La idea general es crear un iterable que contenga tuplas de (geometría, valor), donde la geometría es una geometría con forma y el valor es lo que quieres grabar en el raster en la ubicación de esa geometría. Tanto Fiona como GeoPANDAS utilizan geometrías conformadas, así que estás de suerte. En este ejemplo se utiliza un generador para iterar a través de los pares (geometría,valor) que fueron extraídos del GeoDataFrame y unidos mediante zip().

Asegúrese de abrir el out_fn archivo en w+ porque tendrá que utilizarse para leer y escribir.

Answer 2

Geocube es una nueva herramienta diseñada específicamente para rasterizar datos geopandas que envuelve rasterio. Simplifica el proceso y elimina la necesidad de una plantilla raster.

https://github.com/corteva/geocube

En el contexto del ejemplo anterior:

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

La letra puede establecerse en el marco de datos de este modo:

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

Sin embargo, sólo pueden rasterizarse los valores numéricos. He aquí un ejemplo categórico: https://corteva.github.io/geocube/stable/examples/categorical.html

Así que, en lugar de usar eso, usa los números en formato cadena y conviértelos a enteros:

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

A continuación, rasterice los datos:

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

Por último, expórtalo a una trama:

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