Recorta un ráster con rasterio y geopandas

8

Estoy recortando un conjunto de fotografías aéreas históricas. Estas fotografías tienen grandes áreas negras en los bordes (valor 0). Sin embargo, también hay datos válidos con un valor 0. El flujo de trabajo que estoy usando es:

  1. Cargar ráster con rasterio
  2. Poligonalice el ráster con rasterio.features.shapes ()
  3. Identifique polígonos donde valor = 0 y tamaño> 5000 metros cuadrados
  4. Enmascara las imágenes originales con polígonos, realiza una máscara invertida

Aquí está mi código actual para enmascarar una sola imagen:

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

Cuando ejecuto este código, recibo el siguiente error: AttributeError: 'str' object has no attribute 'get'

La documentación de los estados de rasterio.mask : los polígonos son dictados similares a GeoJSON que especifican los límites de las características en el ráster que se mantendrán. Todos los datos fuera de los polígonos especificados se establecerán en nodata.

Supongo que le estoy dando a rasterio.mask el tipo incorrecto de "Dict similar a GeoJSON". He intentado reformatear el dict de varias maneras sin éxito. ¿Alguien sabe la forma correcta de convertir GeoJSON a un "dict similar a GeoJSON"?

¿O alguien puede proporcionar el formato correcto de un "dict similar a GeoJSON"?

Soy nuevo en rasterio y geopandas.

Rosswin
fuente

Respuestas:

6

El problema está resuelto. El problema fue que leí mal la documentación. En una segunda lectura, la documentación de rasterio.mask establece claramente que los polígonos deberían ser una lista de dictados similares a GeoJSON. Encontré el siguiente fragmento de código para generar esas listas a partir de esta respuesta :

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

Aquí está el código completo que está funcionando:

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
fuente