Cortar uma varredura usando rasterio e geopandas

8

Estou cortando um conjunto de fotografias aéreas históricas. Essas fotografias têm grandes áreas em preto nas bordas (valor 0). No entanto, também existem dados válidos com um valor 0. O fluxo de trabalho que estou usando é:

  1. Carregar raster com rasterio
  2. Poligonize a varredura com rasterio.features.shapes ()
  3. Identifique polígonos onde valor = 0 e tamanho> 5000 metros quadrados
  4. Mascarar as imagens originais com polígonos, executar uma máscara invertida

Aqui está o meu código atual para mascarar uma única imagem:

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

Quando executo esse código, recebo o seguinte erro: AttributeError: 'str' object has no attribute 'get'

A documentação para rasterio.mask declara: Os polígonos são ditados do tipo GeoJSON, especificando os limites dos recursos na varredura a serem mantidos. Todos os dados fora dos polígonos especificados serão definidos como nodata.

Estou assumindo que estou dando ao rasterio.mask o tipo errado de "ditado semelhante ao GeoJSON". Eu tentei reformatar o ditado de várias maneiras, sem sucesso. Alguém sabe a maneira correta de converter GeoJSON para um "ditado semelhante a GeoJSON"?

Ou alguém pode fornecer o formato correto de um "ditado semelhante ao GeoJSON"?

Eu sou novo em rasterio e geopandas.

Rosswin
fonte

Respostas:

6

O problema foi resolvido. O problema foi que eu li mal a documentação. Em uma segunda leitura, a documentação do rasterio.mask afirma claramente que os polígonos devem ser uma lista de ditados do tipo GeoJSON. Encontrei o seguinte trecho de código para gerar essas listas a partir desta resposta :

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

Aqui está o 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
fonte