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".

  • 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.
    – xunilk
    Nov 2 '17 at 17:20

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)

enter image description here

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
print no_data
# 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
     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')

enter image description here

  • Thank you @gene for this complete answer. However, I understand that rasterio is not working well with gdal in the same script which may be a problem for me, plus I need to install rasterio and give a try before to accept your answer.
    – egayer
    Nov 1 '17 at 17:48
  • Hey @gene why did you need to use geoms = [mapping(geoms[0])] instead of just geoms[0]?
    – clifgray
    Dec 30 '18 at 2:12
  • 1
    mapping(geoms[0]) = GeoJSON format of the geometry
    – gene
    Jan 1 '19 at 17:19
  • Hi Gene, I got this error "Invalid field type <type 'numpy.ndarray'>" at the last line (d.to_file)
    – ilFonta
    Jan 31 '19 at 21:17
  • 2
    data = out_image.data[0] threw multi-dimensional sub-views are not implemented for me, but data = out_image[0,:,:] worked. Is this a less efficient or otherwise problematic workaround? Any idea why it would have failed as written?
    – jbaums
    Aug 1 '19 at 6:11

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:

enter image description here

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

  • 1
    Hi Xunilk, what are the input files to your script? I would like to use only the raster and the shapefile with the polygons. Many thanx
    – ilFonta
    Feb 3 '19 at 18:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.