Saya ingin tahu cara mendapatkan semua nilai raster dalam poligon menggunakan gdal atau pygeoprocessing, tanpa membaca seluruh grid sebagai array.
pygeoprocessing dan gdal dapat melakukan statistik zona tetapi hanya min, maks, rata-rata, stdev atau hitungan tersedia dari fungsi tersebut. Karena statistik zona perlu mengakses nilai, apakah akan mudah untuk mengekstraksi nilai dengan cara yang sama?
Saya menemukan pertanyaan yang sangat mirip di sini: ( Mendapatkan nilai piksel raster GDAL di bawah titik OGR tanpa NumPy? ) Tetapi hanya untuk "titik" tertentu.
python
gdal
pygeoprocessing
egayer
sumber
sumber
Jawaban:
Anda dapat menggunakan rasterio untuk mengekstrak nilai raster dalam poligon seperti pada GIS SE: GDAL python memotong gambar geotiff dengan file geojson
Di sini saya menggunakan file raster satu band dan GeoPandas untuk shapefile (bukan Fiona)
Hasil out_image adalah array bertopeng Numpy
Sekarang saya menggunakan Bagaimana cara mendapatkan koordinat sel dalam geotif? atau Python affine mentransformasikan untuk mengubah antara pixel dan koordinat yang diproyeksikan dengan
out_transform
sebagai affine transform untuk data subsetPembuatan GeoDataFrame baru yang dihasilkan dengan nilai col, row, dan elevation
sumber
geoms = [mapping(geoms[0])]
bukan hanyageoms[0]
?mapping(geoms[0])
= Format GeoJSON dari geometridata = out_image.data[0]
melemparkanmulti-dimensional sub-views are not implemented
untuk saya, tetapidata = out_image[0,:,:]
berhasil. Apakah ini solusi yang kurang efisien atau bermasalah? Adakah yang tahu mengapa itu gagal seperti yang tertulis?Jika Anda memiliki masalah dengan menggunakan rasterio dalam skrip yang sama dengan gdal, saya mencoba dengan pygeoprocessing (ini juga menggunakan rupawan) dan saya menemukan solusinya. Skrip lengkap (dengan jalur ke lapisan saya) adalah sebagai berikut:
Setelah menjalankannya, saya mendapat:
di mana nilai-nilai pengambilan sampel raster seperti yang diharapkan di setiap titik dan dimasukkan ke lapisan titik.
sumber