Extrayendo valores ráster dentro de shapefile con pygeoprocessing, gdal, shapely y fiona

El módulo "third party" pygeoprocessing de python facilita algunas funciones que pueden realizarse con GDAL e internamente emplea fiona y shapely por lo que si logran instalarlo sin problemas mediante pip tendrán acceso también a estos módulos.

Lo que si he podido detectar es que no tengo acceso a muchos de los métodos descritos en el manual por lo que sospecho que ha sido pensado más bien para Python 3 o debo forzar la instalación de la última versión con pip en Python 2.7.

El código completo desarrollado puede observarse a continuación:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73

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

Después de ejecutado en la Python Console de QGIS y cargadas todas las capas involucradas (excepto la de estensión del shapefile), tal como se esperaba, se obtuvo lo siguiente:

Mediante el plugin Value Tool de QGIS se corroboró que el muestreo ráster fue ejecutado de la manera correcta.