Extract raster values within shapefile with pygeoprocessing or gdal

Question

I would like to know how to get all the raster values within a polygon using gdal or pygeoprocessing, without reading the entire grid as an array.

pygeoprocessing and gdal can do zonal statistics but only the min, max, mean, stdev or count are available from such a function. Since zonal statistics need to access the values, would it be easy to extract values the same way ?

I found a very similar question here : (Getting pixel value of GDAL raster under OGR point without NumPy?) but only for a particular "point".

Answer 1

You can use rasterio to extract the raster values within a polygon as in GIS SE: GDAL python cut geotiff image with geojson file

I use here a one band raster file and GeoPandas for the shapefile ( instead of Fiona)

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)

The out_image result is a Numpy masked 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)

Now I use How to I get the coordinates of a cell in a geotif? or Python affine transforms to transform between the pixel and projected coordinates with out_transform as the affine transform for the subset data

 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  

Creation of a new resulting GeoDataFrame with the col, row and elevation values

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

Answer 2

If you have problems by using rasterio in the same script with gdal, I was trying out with pygeoprocessing (it also uses shapely) and I found out a workaround. Complete script (with paths to my layers) is as follows:

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

After running it, I got:

where raster sampling values were as expected in each point and incorporated to point layer.