Extrahieren Sie Rasterwerte innerhalb des Shapefiles mit Pygeoprocessing oder GDAL

11

Ich möchte wissen, wie man alle Rasterwerte innerhalb eines Polygons mit gdal oder pygeoprocessing erhält, ohne das gesamte Raster als Array zu lesen.

pygeoprocessing und gdal können zonale Statistiken erstellen, aber von einer solchen Funktion sind nur die Werte min, max, mean, stdev oder count verfügbar. Wäre es einfach, Werte auf dieselbe Weise zu extrahieren, da Zonenstatistiken auf die Werte zugreifen müssen?

Ich habe hier eine sehr ähnliche Frage gefunden: ( Pixelwert des GDAL-Rasters unter OGR-Punkt ohne NumPy abrufen? ), Aber nur für einen bestimmten "Punkt".

egayer
quelle
Wenn Sie Probleme mit der Verwendung von Rasterio im selben Skript mit gdal haben, habe ich Pygeoprocessing ausprobiert (es wird auch Shapely verwendet) und eine Problemumgehung gefunden.
Xunilk

Antworten:

16

Sie können Rasterio verwenden , um die Rasterwerte innerhalb eines Polygons wie in GIS SE: GDAL Python Cut Geotiff-Bild mit Geojson-Datei zu extrahieren

Ich verwende hier eine Ein-Band-Rasterdatei und GeoPandas für das Shapefile (anstelle von Fiona).

Geben Sie hier die Bildbeschreibung ein

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)

Das Ergebnis out_image ist ein Numpy-maskiertes Array

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

Jetzt verwende ich Wie erhalte ich die Koordinaten einer Zelle in einem Geotif? oder affine Python-Transformationen, um zwischen dem Pixel und den projizierten Koordinaten out_transform als affine Transformation für die Teilmengen-Daten zu transformieren

 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  

Erstellung eines neuen resultierenden GeoDataFrame mit den Werten col, row und height

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

Geben Sie hier die Bildbeschreibung ein

Gen
quelle
Vielen Dank an @gene für diese vollständige Antwort. Ich verstehe jedoch, dass Rasterio mit gdal im selben Skript nicht gut funktioniert, was für mich ein Problem sein kann. Außerdem muss ich Rasterio installieren und es vorher versuchen, um Ihre Antwort zu akzeptieren.
Egayer
Hey @gene, warum musstest du geoms = [mapping(geoms[0])]statt nur verwenden geoms[0]?
Clifgray
1
mapping(geoms[0])= GeoJSON Format des Geometrie
Gen
Hallo Gene, ich habe diesen Fehler "Ungültiger Feldtyp <Typ 'numpy.ndarray'>" in der letzten Zeile (d.to_file)
ilFonta
1
data = out_image.data[0]warf multi-dimensional sub-views are not implementedfür mich, data = out_image[0,:,:]arbeitete aber . Ist dies eine weniger effiziente oder anderweitig problematische Problemumgehung? Irgendeine Idee, warum es wie geschrieben fehlgeschlagen wäre?
Jbaums
2

Wenn Sie Probleme mit der Verwendung von Rasterio im selben Skript mit gdal haben, habe ich Pygeoprocessing ausprobiert (es wird auch Shapely verwendet) und eine Problemumgehung gefunden. Das vollständige Skript (mit Pfaden zu meinen Ebenen) lautet wie folgt:

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

Nachdem ich es ausgeführt habe, habe ich:

Geben Sie hier die Bildbeschreibung ein

Dabei waren die Rasterabtastwerte in jedem Punkt wie erwartet und wurden in die Punktebene integriert.

xunilk
quelle
Hallo Xunilk, was sind die Eingabedateien für dein Skript? Ich möchte nur das Raster und das Shapefile mit den Polygonen verwenden. Vielen Dank
ilFonta