Calculate average value in GeoTIFF file inside a shapefile polygon?

Question

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 m³ / hectare for each pixel.

I would like to calculate the average m³ / 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.

Answer 1

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:

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

Answer 2

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   
[snip]
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).