Extraire des valeurs raster dans un fichier de formes avec pygeoprocessing ou gdal

12

Je voudrais savoir comment obtenir toutes les valeurs raster dans un polygone en utilisant gdal ou pygeoprocessing, sans lire la grille entière comme un tableau.

pygeoprocessing et gdal peuvent faire des statistiques zonales mais seuls les min, max, moyenne, stdev ou count sont disponibles à partir d'une telle fonction. Étant donné que les statistiques zonales doivent accéder aux valeurs, serait-il facile d'extraire les valeurs de la même manière?

J'ai trouvé une question très similaire ici: ( Obtenir la valeur en pixels du raster GDAL sous le point OGR sans NumPy? ) Mais uniquement pour un "point" particulier.

égayer
la source
Si vous rencontrez des problèmes en utilisant rasterio dans le même script avec gdal, j'essayais avec le pygeoprocessing (il utilise également galbé) et j'ai trouvé une solution de contournement.
xunilk

Réponses:

16

Vous pouvez utiliser rasterio pour extraire les valeurs raster dans un polygone comme dans GIS SE: image géotiff coupée en python GDAL avec fichier geojson

J'utilise ici un fichier raster à une bande et GeoPandas pour le shapefile (au lieu de Fiona)

entrez la description de l'image ici

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)

Le résultat out_image est un tableau masqué 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)

Maintenant, j'utilise Comment obtenir les coordonnées d'une cellule dans un géotif? ou Python affine transforme pour se transformer entre le pixel et les coordonnées projetées avec out_transform comme transformation affine pour les données du sous-ensemble

 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  

Création d'une nouvelle GeoDataFrame résultante avec les valeurs de col, de ligne et d'élévation

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

entrez la description de l'image ici

gène
la source
Merci @gene pour cette réponse complète. Cependant, je comprends que rasterio ne fonctionne pas bien avec gdal dans le même script, ce qui peut être un problème pour moi, et je dois installer rasterio et essayer avant d'accepter votre réponse.
egayer
Hé @gene, pourquoi avez-vous dû utiliser geoms = [mapping(geoms[0])]plutôt que juste geoms[0]?
clifgray
1
mapping(geoms[0])= Format GeoJSON de la géométrie
gène
Salut Gene, j'ai eu cette erreur "Type de champ non valide <type 'numpy.ndarray'>" à la dernière ligne (d.to_file)
ilFonta
1
data = out_image.data[0]jeté multi-dimensional sub-views are not implementedpour moi, mais a data = out_image[0,:,:]fonctionné. Est-ce une solution de contournement moins efficace ou autrement problématique? Une idée de pourquoi cela aurait échoué tel qu'écrit?
jbaums
2

Si vous rencontrez des problèmes en utilisant rasterio dans le même script avec gdal, j'essayais avec le pygeoprocessing (il utilise également galbé ) et j'ai trouvé une solution de contournement. Le script complet (avec les chemins vers mes couches) est le suivant:

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

Après l'avoir exécuté, j'ai eu:

entrez la description de l'image ici

où les valeurs d'échantillonnage raster étaient comme prévu en chaque point et incorporées à la couche de points.

xunilk
la source
Salut Xunilk, quels sont les fichiers d'entrée de votre script? Je souhaite utiliser uniquement le raster et le fichier de formes avec les polygones. Many
thanx