Extraia valores raster dentro do shapefile com pygeoprocessing ou gdal

11

Gostaria de saber como obter todos os valores de varredura em um polígono usando gdal ou pygeoprocessing, sem ler toda a grade como uma matriz.

pygeoprocessing e gdal podem fazer estatísticas zonais, mas apenas o min, max, mean, stdev ou count estão disponíveis nessa função. Como as estatísticas zonais precisam acessar os valores, seria fácil extrair valores da mesma maneira?

Eu encontrei uma pergunta muito semelhante aqui: ( Obtendo valor de pixel da varredura GDAL sob o ponto OGR sem NumPy? ), Mas apenas para um "ponto" específico.

egayer
fonte
Se você tiver problemas usando o rasterio no mesmo script que o gdal, eu estava testando o pygeoprocessing (também usa formas bem torneadas) e descobri uma solução alternativa.
Xunilk

Respostas:

16

Você pode usar o rasterio para extrair os valores do raster dentro de um polígono, como no GIS SE: GDAL, imagem de geotiff de corte python com arquivo geojson

Eu uso aqui um arquivo raster de uma banda e GeoPandas para o shapefile (em vez de Fiona)

insira a descrição da imagem aqui

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)

O resultado out_image é um array mascarado 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)

Agora eu uso Como obter as coordenadas de uma célula em um geotif? ou Python affine transforma para transformar entre o pixel e as coordenadas projetadas com out_transform a transformação affine para os dados do subconjunto

 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  

Criação de um novo GeoDataFrame resultante com os valores de col, linha e elevação

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

insira a descrição da imagem aqui

gene
fonte
Obrigado @gene por esta resposta completa. No entanto, entendo que o rasterio não está funcionando bem com o gdal no mesmo script, o que pode ser um problema para mim, além disso, preciso instalar o rasterio e tentar antes para aceitar sua resposta.
exemplo
Hey @gene, por que você precisou usar em geoms = [mapping(geoms[0])]vez de apenas geoms[0]?
Clifgray
1
mapping(geoms[0])= Formato GeoJSON da geometria
gene
Oi Gene, recebi este erro "Tipo de campo inválido <tipo 'numpy.ndarray'>" na última linha (arquivo.d.to)
ilFonta
1
data = out_image.data[0]jogou multi-dimensional sub-views are not implementedpara mim, mas data = out_image[0,:,:]funcionou. Esta é uma solução alternativa menos eficiente ou problemática? Alguma idéia de por que teria falhado como está escrito?
jbaums
2

Se você tiver problemas usando o rasterio no mesmo script que o gdal, eu estava testando o pygeoprocessing (também usa formas bem torneadas) e descobri uma solução alternativa. O script completo (com caminhos para minhas camadas) é o seguinte:

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

Depois de executá-lo, obtive:

insira a descrição da imagem aqui

onde os valores de amostragem raster foram os esperados em cada ponto e incorporados à camada de pontos.

xunilk
fonte
Olá Xunilk, quais são os arquivos de entrada para o seu script? Eu gostaria de usar apenas o raster e o shapefile com os polígonos. Muitos thanx
ilFonta 03/02/19