Ekstrak nilai raster di dalam shapefile dengan pygeoprocessing atau gdal

12

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.

egayer
sumber
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.
xunilk

Jawaban:

16

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)

masukkan deskripsi gambar di sini

import rasterio
from rasterio.mask import mask
import geopandas as gpd
shapefile = gpd.read_file("extraction.shp")
# extract the geometry in GeoJSON format
geoms = shapefile.geometry.values # list of shapely geometries
geometry = geoms[0] # shapely geometry
# transform to GeJSON format
from shapely.geometry import mapping
geoms = [mapping(geoms[0])]
# extract the raster values values within the polygon 
with rasterio.open("raster.tif") as src:
     out_image, out_transform = mask(src, geoms, crop=True)

Hasil out_image adalah array bertopeng Numpy

# no data values of the original raster
no_data=src.nodata
print no_data
-9999.0
# extract the values of the masked array
data = out_image.data[0]
# extract the row, columns of the valid values
import numpy as np
row, col = np.where(data != no_data) 
elev = np.extract(data != no_data, data)

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 subset

 from rasterio import Affine # or from affine import Affine
 T1 = out_transform * Affine.translation(0.5, 0.5) # reference the pixel centre
 rc2xy = lambda r, c: (c, r) * T1  

Pembuatan GeoDataFrame baru yang dihasilkan dengan nilai col, row, dan elevation

d = gpd.GeoDataFrame({'col':col,'row':row,'elev':elev})
# coordinate transformation
d['x'] = d.apply(lambda row: rc2xy(row.row,row.col)[0], axis=1)
d['y'] = d.apply(lambda row: rc2xy(row.row,row.col)[1], axis=1)
# geometry
from shapely.geometry import Point
d['geometry'] =d.apply(lambda row: Point(row['x'], row['y']), axis=1)
# first 2 points
d.head(2)
     row  col   elev       x          y            geometry  
 0    1    2  201.7!  203590.58  89773.50  POINT (203590.58 89773.50)  
 1    1    3  200.17  203625.97  89773.50  POINT (203625.97 89773.50)

# save to a shapefile
d.to_file('result.shp', driver='ESRI Shapefile')

masukkan deskripsi gambar di sini

gen
sumber
Terima kasih @ gen atas jawaban lengkap ini. Namun, saya mengerti bahwa rasterio tidak bekerja dengan baik dengan gdal dalam skrip yang sama yang mungkin menjadi masalah bagi saya, ditambah saya perlu menginstal rasterio dan mencoba sebelumnya untuk menerima jawaban Anda.
egayer
Hey @gene, mengapa Anda harus menggunakan, geoms = [mapping(geoms[0])]bukan hanya geoms[0]?
clifgray
1
mapping(geoms[0])= Format GeoJSON dari geometri
gen
Hai Gene, saya mendapatkan kesalahan ini "Jenis bidang tidak valid <ketik 'numpy.ndarray'>" di baris terakhir (d.to_file)
ilFonta
1
data = out_image.data[0]melemparkan multi-dimensional sub-views are not implementeduntuk saya, tetapi data = out_image[0,:,:]berhasil. Apakah ini solusi yang kurang efisien atau bermasalah? Adakah yang tahu mengapa itu gagal seperti yang tertulis?
jbaums
2

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:

import pygeoprocessing.geoprocessing as geop
from shapely.geometry import shape, mapping, Point
from osgeo import gdal
import numpy as np 
import fiona

path = '/home/zeito/pyqgis_data/'

uri1 = path + 'aleatorio.tif'

info_raster2 = geop.get_raster_info(uri1)

geop.create_raster_from_vector_extents(base_vector_path = path + 'cut_polygon3.shp',
                                       target_raster_path = path + 'raster_from_vector_extension.tif',
                                       target_pixel_size = info_raster2['pixel_size'],
                                       target_pixel_type = info_raster2['datatype'],
                                       target_nodata = -999,
                                       fill_value = 1)

uri2 = path + 'raster_from_vector_extension.tif'

info_raster = geop.get_raster_info(uri2)

cols = info_raster['raster_size'][0]
rows = info_raster['raster_size'][1]

geotransform = info_raster['geotransform']

xsize =  geotransform[1]
ysize = -geotransform[5]

xmin = geotransform[0]
ymin = geotransform[3]

# create one-dimensional arrays for x and y
x = np.linspace(xmin + xsize/2, xmin + xsize/2 + (cols-1)*xsize, cols)
y = np.linspace(ymin - ysize/2, ymin - ysize/2 - (rows-1)*ysize, rows)

# create the mesh based on these arrays
X, Y = np.meshgrid(x, y)

X = X.reshape((np.prod(X.shape),))
Y = Y.reshape((np.prod(Y.shape),))

coords = zip(X, Y)

shapely_points = [ Point(point[0], point[1]) for point in coords ]

polygon = fiona.open(path + 'cut_polygon3.shp')
crs = polygon.crs
geom_polygon = [ feat["geometry"] for feat in polygon ]

shapely_geom_polygon = [ shape(geom) for geom in geom_polygon ]

within_points = [ (pt.x, pt.y) for pt in shapely_points if pt.within(shapely_geom_polygon[0]) ]

src_ds = gdal.Open(uri1)
rb = src_ds.GetRasterBand(1)

gt = info_raster2['geotransform']

values = [ rb.ReadAsArray(int((point[0] - gt[0]) / gt[1]), #x pixel
                          int((point[1] - gt[3]) / gt[5]), #y pixel
                          1, 1)[0][0] 
           for point in within_points ]

#creation of the resulting shapefile
schema = {'geometry': 'Point','properties': {'id': 'int', 'value':'int'},}

with fiona.open('/home/zeito/pyqgis_data/points_proc.shp', 'w', 'ESRI Shapefile', schema, crs)  as output:

    for i, point in enumerate(within_points):
        output.write({'geometry':mapping(Point(point)),'properties': {'id':i, 'value':str(values[i])}})

Setelah menjalankannya, saya mendapat:

masukkan deskripsi gambar di sini

di mana nilai-nilai pengambilan sampel raster seperti yang diharapkan di setiap titik dan dimasukkan ke lapisan titik.

xunilk
sumber
Hai Xunilk, apa file input untuk skrip Anda? Saya hanya ingin menggunakan raster dan shapefile dengan poligon. Many thanx
ilFonta