Pixelliser un fichier de formes avec Geopandas ou fiona - python

J'ai besoin de pixelliser un fichier de formes vraiment simple un peu comme ceci http://tinyurl.com/odfbanu . Ce qui est juste un fichier de formes contenant des pays aux États-Unis.J'ai vu cette réponse précédente: GDAL RasterizeLayer ne brûle pas tous les polygones en raster? mais je me demandais s'il y avait un moyen de le faire en utilisant Geopandas ou fiona et peut-être rasterio pour la partie d'écriture tiff.

Mon objectif est donc de le pixelliser et d'attribuer une valeur à chaque polygone partageant une valeur commune, LSAD dans l'exempe.

J'ai donc écrit le début du code inspiré du shongololo dans le fil: Dissoudre des polygones basés sur des attributs avec Python (galbe, 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'

Des trucs vraiment faciles alors maintenant je me demande comment puis-je réellement écrire ces formes sur un tiff. J'ai commencé à travailler avec Geopandas car je pensais que c'était la meilleure option, mais si vous avez une suggestion de Fiona, je suis d'accord.

J'ai trouvé un morceau de code de rasterio qui semble pouvoir prendre une géométrie galbée et la graver dans un nouveau 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))

Vous êtes sur la bonne voie et les géopandas GeoDataFrame sont un bon choix pour la pixellisation sur Fiona. Fiona est un excellent ensemble d'outils, mais je pense que le DataFrame est mieux adapté aux fichiers de formes et aux géométries que les dictionnaires imbriqués.

import geopandas as gpd
import rasterio
from rasterio import features

Configurez vos noms de fichiers

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

Ouvrez le fichier avec GeoPANDAS read_file

counties = gpd.read_file(shp_fn)

Ajoutez la nouvelle colonne (comme dans votre code ci-dessus)

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'

Ouvrez le fichier raster que vous souhaitez utiliser comme modèle pour la gravure d'entités à l'aide de rasterio

rst = rasterio.open(rst_fn)

copier et mettre à jour les métadonnées du raster en entrée pour la sortie

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

Maintenant, gravez les entités dans le raster et écrivez-les

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)

L'idée générale est de créer un itérable contenant des tuples de (géométrie, valeur), où la géométrie est une géométrie galbée et la valeur est ce que vous voulez graver dans le raster à l'emplacement de cette géométrie. Fiona et GeoPANDAS utilisent des géométries galbées, vous avez donc de la chance. Dans cet exemple, un générateur est utilisé pour parcourir les paires (géométrie, valeur) qui ont été extraites du GeoDataFrame et réunies à l'aide de zip ().

Assurez-vous d'ouvrir le out_fnfichier en w+mode, car il devra être utilisé pour la lecture et l'écriture.

geocube est un nouvel outil spécialement conçu pour rasteriser les données géopandas qui encapsule rasterio. Il simplifie le processus et élimine le besoin d'un raster de modèle.

https://github.com/corteva/geocube

Dans le cadre de l'exemple ci-dessus:

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 lettre peut être définie sur la trame de données comme suit:

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

Cependant, seules les valeurs numériques peuvent être pixellisées. Voici un exemple catégorique: https://corteva.github.io/geocube/stable/examples/categorical.html

Donc, au lieu d'utiliser cela, utilisez les nombres au format chaîne et convertissez-les en entiers:

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

Ensuite, pixellisez les données:

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

Enfin, exportez-le vers un raster:

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