I have raster files with precipitation data (4x4 km), monthly, for about 100 years. I have a shapefile for a watershed (which includes several polygons as sub-basins). I need the average precipitation for the watershed for each month. I don't want to do this by hand in ArcMap because it will probably need to be repeated plenty of times (for other watersheds, and when more data becomes available), so I am coding it.

I have written code in R now, using packages raster and rgdal. It works great in the sense that it has a functionality to calculate the weighted average per sub-basin; it checks which fraction of a raster cell falls within the polygon and uses this to calculate a weighted average for the polygon. Because some of the sub-basins are small in comparison to the raster, this is a great feature. And in addition, because I can access the shapefile information such as the area of the sub-basins, I can calculate a weighted average for the whole watershed.

However - this is slow. It also says it in the description of the function extract(). Calculating the weighted value takes 4-5 seconds, so when I ran this script for 30 years worth of monthly data, it took 37 minutes (there are a few more calculations besides the weighted average).

I have a feeling it might be faster to use Python, but I am not sure how to approach this. I did get started with this before the R code (see below) but shifted my attention to R. It seemed a bit daunting since it's been a while that I coded in Python, and I've never used Python for geospatial analyses.

Also, an initial test seems to show that in Python, there is no weighted averaging happening, and that small polygons that do not contain the center of a raster cell, get assigned 0 or None.

So my questions are:

  1. If it is possible in Python, will weighted averaging that way be faster?
  2. How to access other data associated with the shapefile (such as area)?

Below is the very very basic code I have so far.

    from rasterstats import zonal_stats
    import gdal


    # NOTE From Python documentation: "There is no right or wrong way to rasterize a vector.
    # The default strategy is to include all pixels along the line render path (for lines),
    # or cells where the center point is within the polygon (for polygons)

    # get average for every polygon (there are 19)
    data=([f['mean'] for f in stats])

    # take out the "None" values...
    data2=[x for x in data if x is not None]

    # ... so I can calculate the mean for the whole watershed

data looks as follows. Note the None...


My system specs: Windows 7, Python 2.7.9 (using the Enthought Canopy IDE) 32-bit, GDAL-1.11.3, rasterio-0.29.0, Shapely-1.5.13

  • 1
    Welcome to GIS SE! As a new user please take the Tour where you will see that there should be only one question per question. I've collapsed your first two into one but I think that you should edit the remaining one out to be researched/asked separately.
    – PolyGeo
    Nov 6, 2015 at 22:23
  • 1
    If you showed your R code, we might be able to help speed it up. Nov 9, 2015 at 17:11
  • 1
    I would ask the question, why bother? You are adding huge computational overhead to a problem that has a notable scale mismatch. Have you looked at the difference between the overall distribution using the weighted and non-weighted statistic? I would imagine that it matters very little. Certainty, there some proceses and resolutions where this would be a major concern but, I just do not see it in climate aggregrations to watershed basins. Climate data is often quite locally homogenious and weighting the mean by pixel fraction seems like overkill. Sep 6, 2017 at 16:16
  • Did have a look at gis.stackexchange.com/questions/226983/…? This could help, I suppose.
    – yenats
    Oct 26, 2017 at 11:12
  • I wrote an R package for this specific problem (zonal stats with partially covered pixels). It is much faster than the raster package in cases like this.
    – dbaston
    Jul 3, 2018 at 15:26

1 Answer 1


"If it is possible in Python, will weighted averaging that way be faster?"

See answer to next question. Collecting the aggregated stats along with the rasterizing under the hood is what takes the time. You can get your weighted average using the following answer. As far as what is faster, this I cannot comment on...if you're willing to re-code rasterstats, it can be made a fair bit faster but you would lose its nice-to-useness.

"How to access other data associated with the shapefile (such as area)?"

For the data returned by zonal_stats, there is a "count" value which, I believe, is the number of rasterized cells that contribute. If you know the resolution of your raster (and hence area per pixel) then you have area (area = area_per_pixel * 'count').

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