Pangkas raster menggunakan rasterio dan geopanda

8

Saya memangkas satu set foto udara bersejarah. Foto-foto ini memiliki area hitam besar di tepinya (nilai 0). Namun ada juga data yang valid dengan nilai 0. Alur kerja yang saya gunakan adalah:

  1. Memuat raster dengan rasterio
  2. Polygonize raster dengan rasterio.features.shapes ()
  3. Identifikasi poligon di mana nilai = 0 dan ukuran> 5000 meter persegi
  4. Menyamarkan gambar asli dengan poligon, melakukan topeng terbalik

Ini kode saya saat ini untuk menutupi satu gambar:

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

Ketika saya menjalankan kode ini saya menerima kesalahan berikut: AttributeError: 'str' object has no attribute 'get'

The dokumentasi untuk rasterio.mask negara: Poligon yang GeoJSON seperti dicts menentukan batas-batas fitur dalam raster yang akan disimpan. Semua data di luar poligon yang ditentukan akan ditetapkan ke nodata.

Saya berasumsi bahwa saya memberikan rasterio.mask jenis yang salah dari "dict GeoJSON-like". Saya telah mencoba untuk memformat ulang dikt beberapa cara tanpa hasil. Adakah yang tahu cara yang benar untuk mengonversi GeoJSON menjadi "dict seperti GeoJSON"?

Atau adakah yang bisa memberikan format "geoJSON-like dict" yang benar?

Saya baru mengenal rasterio dan geopanda.

rosswin
sumber

Jawaban:

6

Masalah telah teratasi. Masalahnya adalah saya salah membaca dokumentasi. Pada bacaan kedua, dokumentasi rasterio.mask dengan jelas menyatakan bahwa poligon harus menjadi daftar dict seperti GeoJSON. Saya menemukan potongan kode berikut untuk menghasilkan daftar dari jawaban ini :

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

Berikut adalah kode lengkap yang berfungsi:

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
sumber