I'm a forest owner and have information about my forest areas in ESRI shapefile format, in ETRS-TM35FIN coordinate system. I managed to create a text dump about the forest area in the following style:

Shape:10 (Polygon)  nVertices=831, nParts=1
  Bounds:(365190.963,6785661.028, 0)
      to (372638.972,6788929.792, 0)
     (366426.668,6787577.830, 0) Ring
     (366426.948,6787576.866, 0)
     (366427.680,6787575.980, 0)
     (366428.396,6787574.551, 0)
     (366426.668,6787577.830, 0)

(the data has been anonymized, by the way, so if it doesn't look like valid data in ETRS-TM35FIN system, I subtracted a random offset from x/y values, same for each x, same for each y)

I also downloaded a GeoTIFF file which contains the amount of trees in m3 / hectare for each pixel.

I would like to calculate the average m3 / hectare value inside my forest area, discarding all information about neighboring forest areas.

How to do this? Is there some existing software that is able to load GeoTIFF data and calculate averages within ESRI shapefile polygons? Do I have to write the software on my own? I'm a very capable software developer, so given pointers to the right direction (i.e. suitable libraries), I'm probably able to write the software required to calculate the average.

I'd like to restrict answers to free and very cheap software, i.e. I'm not willing to pay more than $100 for the required software, and I suspect many GIS software systems are far more expensive than that.


2 Answers 2


The process of summarizing raster data within the bounds of some 2d vectors is commonly called zonal statistics and can be accomplished using open source software in a number of ways. For example:

  1. QGIS 3.x has the zonal statistics tool built in: enter image description here

  2. In python you can use the package "rasterstats" https://pypi.org/project/rasterstats/

  3. In R you can use the packages 'maptools', 'rgdal' and 'raster' to perform the extraction. Use 'foreach' and 'doparalell' to process the operations on multiple cores. https://rpubs.com/rural_gis/254726


I already managed to write a Python script so I'll describe how I did it:

Use numpy as np, imageio, shapely.geometry.Point and shapely.geometry.polygon.Polygon.

All of the following assumes the shapefile and GeoTIFF are in the same coordinate system (now they're in ETRS-TM35FIN both so the requirement is met). If they're not in the same coordinate system, you need to project the shapefile coordinates to the GeoTIFF coordinate system.

Do the following once (this is slow for large polygons):

Use the free shpdump tool to dump the shapefile to the text format. Create a Python Polygon from the text data. (Some manual editing required, but can be done with regular expressions. The Python Polygon doesn't require last and first coordinate to be the same, so you can comment out the last coordinate identical to the first.)

Use the listgeo tool to show the bounds of the image:

      ModelPixelScaleTag (1,3):
         16                16                0   
Upper Left    (  500000.000, 7338000.000)
Lower Left    (  500000.000, 7242000.000)
Upper Right   (  692000.000, 7338000.000)
Lower Right   (  692000.000, 7242000.000)
Center        (  596000.000, 7290000.000)

Now we know that fileleft = 500000 and fileupper = 7338000 and scaling = 16.0.

Create a np.zeros mask array in the shape of the GeoTIFF image, using np.bool data type.

Calculate the minimum and maximum coordinates of the polygon within the image. Add/subtract few pixels for safety.

For each x,y pixel within the minimum/maximum coordinates, create a Point in the ETRS-TM35FIN coordinate system describing the center of a pixel. Check whether it's in the Polygon, if so, set mask value to 1. Note the midpoint of a pixel is fileleft + (x+0.5)*scaling and fileupper - (y+0.5)*scaling, so do remember the 0.5.

Store the numpy mask to a file, preferably gzipped (it gzips well).

Do the following once in 2 years when a new GeoTIFF arrives describing new forest growth:

Load the numpy mask array to mask.

Load the GeoTIFF image to im using im = imageio.imread(fname).

Calculate np.sum(im*mask)/np.sum(mask).

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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