I have a polygon and filelist of data I need to convert my polygon to raster and mask my filelist with this mask and campute the statistics in the area in polygon. I'm using GDAL Python.

I can do it using R like this :

maskraster <- rasterize(polygon, raster, mask = TRUE) 
#calcul mean of area in polygons 
Rasternew <- cellStats(maskraster,stat=mean)

How can I do the same using GDAL in Python ?


I have a script which allows me to clip a raster to the polygon of a shapefile, and then compute statistics on it (including, as in this case) the mean. This example includes the modules pillow (https://pypi.org/project/Pillow/) and shapefile (https://pypi.org/project/pyshp/) on Python as well.

from osgeo import gdal, osr
import shapefile
from PIL import Image, ImageDraw

def world_to_image(lat_value,lon_value,dataset):
    This function is required for converting from lat/lon values to the corresponding values in the image coordinates.
    metadata = dataset.GetGeoTransform()
    # metadata is an array with six values: ulx, xres, xskew, uly, yskew, yres

    # The source projection (lat, long coords)
    source = osr.SpatialReference()

    # The target projection (projection of the image)
    target = osr.SpatialReference()

    # Create the transform - this can be used repeatedly
    transform = osr.CoordinateTransformation(source,target)
    new_point = transform.TransformPoint(long_value,lat_value)

    x_pixel = int( round((new_point[0] - metadata[0]) / metadata[1],0)) #matrix[:,*value*]
    y_pixel = int( round((new_point[1] - metadata[3]) / metadata[5],0)) #matrix[*value*,:]
    return y_pixel,x_pixel

dataset = gdal.Open("path/to/file")
band = dataset.GetRasterBand("band number")
band = band.ReadAsArray()
band = band.astype(float) #usually ensures no weirdness later when doing maths steps

sf = shapefile.Reader("path/to/shapefile")
minX, minY, maxX, maxY = sf.bbox #Get corner coordinates of shapefile

ulX,ulY = world_to_image(maxY,minX,dataset) #Convert shapefile corners to image 
lrX,lrY = world_to_image(minY,maxX,dataset)

pxwidth = int(lrX-ulX) #number of pixels for width of image
pxheight = int(lrY-ulY) #number of pixels for height of image
clip = band[ulY:lrY,ulX:lrX] #clipping original image to shapefile bounds

pixels = []
for p in sf.shape(0).points:
    temp = world_to_image(p[1],p[0],dataset)
    pixels.append((temp[1],temp[0])) #convert the points of the polygon to image coordinates 

rasterpoly = Image.new('L',(band.shape),1) #set up new image in the same shape as clipped raster image
rasterize = ImageDraw.Draw(rasterpoly) 
rasterize.polygon(pixels,0) #cuts empty image to exact shape of shapefile
mask = np.array(rasterpoly) #converts it to numpy array

function which cuts the original raster image to the shape of the shapefile, 
leaving nodata points elsewhere.
clip = gdal_array.numpy.choose(mask, (band, np.nan)).astype(gdal_array.numpy.float)

Apologies for the rather lengthy code snippet - it is cut from a longer body of code but I have tried to cut out everything unnecessary. This code should produce a 2D numpy array which you can then perform your statistical functions on. Apologies if this is not exactly the sort of thing you were looking for - it was quite hard to know exactly what you wanted from your post (as someone not familiar with R and therfore your example).

  • that what i need calcul cellstat that is my new question gis.stackexchange.com/questions/308733/… Jan 15 '19 at 18:06
  • as far as I can see that is a completely different question... If my answer above is correct please can you accept the answer :) I will have a look at the other question when I can.
    – JackLidge
    Jan 16 '19 at 10:30
  • that is good thank you and i have some questions cause i don't undersand why you are doing # The target projection (UTM zone WGS 84) , and your exemple is good for one raster with one band but my how i can do that for all 10 bands in my raster ? Jan 16 '19 at 11:36
  • Ahh, that isn't necessarily true that comment, in truth that line of code will convert to whatever coordinate system is used by the input image, I will edit my answer for clarity. I haven't tried with multiple bands, but you could try putting the .GetRasterBand command into a for loop and loading them all into one numpy array. I don't know how this would affect the code below though.
    – JackLidge
    Jan 16 '19 at 11:57
  • in this line clip = gdal_array.numpy.choose(mask, (band, np.nan)).astype(gdal_array.numpy.float) this fonction is for clipping image with mask without nan value ( 0 value ) ? Jan 16 '19 at 15:18

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.