Recadrer un raster à l'aide de rasterio et de géopandas

8

Je recadre un ensemble de photographies aériennes historiques. Ces photographies ont de grandes zones noires sur les bords (valeur 0). Cependant, il existe également des données valides avec une valeur 0. Le flux de travail que j'utilise est:

  1. Charger un raster avec rasterio
  2. Polygoniser le raster avec rasterio.features.shapes ()
  3. Identifier les polygones dont la valeur = 0 et la taille> 5000 mètres carrés
  4. Masquez les images originales avec des polygones, effectuez un masque inversé

Voici mon code actuel pour masquer une seule image:

import rasterio
from rasterio import features
from rasterio import mask
import json
import geopandas as gpd

results = []
final_results = []

with rasterio.open(r"C:\1927_oahu\tif\_Line1_6to8_0.tif") as src:
    src_meta = src.meta
    src_affine = src_meta.get("transform")

    band = src.read(1)

    for geometry, raster_value in features.shapes(band, transform=src_affine):
        if raster_value == 0:
        result = {'properties': {'raster_value': raster_value}, 'geometry': geometry}
        results.append(result)

gpd_results = gpd.GeoDataFrame.from_features(results)

gpd_results["area"] = gpd_results["geometry"].area

gpd_results_filtered = gpd_results[gpd_results["area"] > 5000]

gpd_filtered_json_str = gpd_results_filtered.to_json()

gpd_filtered_json_dict = json.loads(gpd_filtered_json_str)

for k, v in gpd_filtered_json_dict.iteritems():
    if k == "features":
        for items in v:
            #final_results = {"coordinates": (items.get("geometry").get("coordinates"))}
            final_results = {"geometry": (items.get("geometry").get("coordinates"))}

masked_band, masked_transform = mask.mask(src, final_results, invert=True)

src_meta.update(dtype=rasterio.uint8, nodata=0)
with rasterio.open(os.path.join(r"C:\1927_oahu_output", "out.tif"), 'w', **src_meta) as dst:
    dst.write_band(1, masked_band.astype(rasterio.uint8))

Lorsque j'exécute ce code, je reçois l'erreur suivante: AttributeError: 'str' object has no attribute 'get'

La documentation de rasterio.mask indique: Les polygones sont des modèles de type GeoJSON spécifiant les limites des entités dans le raster à conserver. Toutes les données en dehors des polygones spécifiés seront définies sur nodata.

Je suppose que je donne à rasterio.mask le mauvais type de "dict de type GeoJSON". J'ai essayé de reformater le dict de plusieurs façons sans succès. Quelqu'un connaît-il la bonne façon de convertir GeoJSON en un "dict de type GeoJSON"?

Ou quelqu'un peut-il fournir le format correct d'un "dict de type GeoJSON"?

Je suis nouveau dans la rasterio et les géopandas.

Rosswin
la source

Réponses:

6

Le problème est résolu. Le problème était que j'avais mal lu la documentation. Lors d'une deuxième lecture, la documentation de rasterio.mask indique clairement que les polygones doivent être une liste de dict de type GeoJSON. J'ai trouvé l'extrait de code suivant pour générer ces listes à partir de cette réponse :

geoms = [feature["geometry"] for feature in shapefile]

Voici le code complet qui fonctionne:

import rasterio
from rasterio import features
from rasterio import mask
import json
import geopandas as gpd
import os

results = []
final_results = []

with rasterio.open(r"C:\1927_oahu\tif\_Line1_6to8_0.tif") as src:
    src_meta = src.meta.copy()
    src_affine = src_meta.get("transform")

    band = src.read(1)

    for geometry, raster_value in features.shapes(band, transform=src_affine):
        if raster_value == 1:
            result = {'properties': {'raster_value': raster_value}, 'geometry': geometry}
            results.append(result)

        gpd_results = gpd.GeoDataFrame.from_features(results)

        gpd_results["area"] = gpd_results["geometry"].area

        gpd_results_filtered = gpd_results[gpd_results["area"] > 5000]

        gpd_filtered_json_str = gpd_results_filtered.to_json()

        gpd_filtered_json_dict = json.loads(gpd_filtered_json_str)

        final_results = [feature["geometry"] for feature in gpd_filtered_json_dict["features"]]

        masked_band, masked_transform = mask.mask(src, final_results, invert=True)

        masked_band[masked_band > 255] = 0

        src_meta.update(dtype=rasterio.uint8, height=int(masked_band.shape[1]), width=int(masked_band.shape[2]), nodata=0, transform=masked_transform, compress='lzw')

        with rasterio.open(r"C:\1927_oahu\output\_Line1_6to8_0.tif"), 'w', **src_meta) as dst:
                dst.write(masked_band.astype(rasterio.uint8))   
Rosswin
la source