Last year Hansen et al. published Global Forest Change map as set of raster data available in TIFF format. Basically those are 30-meters resolution satellite images with pixel values ranging from 0 (no forest) to 100 (full forest).

I would like to programmatically calculate tree cover area over certain country that has its boundaries described by a shapefile. Since I just started with GIS, would you be so kind and advise me where should I begin? What are the available libraries that could help me?

I guess my problem is not new at all, but with my limited knowledge I was unable to find proper resources (or recognise them as such). My preference would be C++ or Python, given that the latter has some high-performance libraries.

For starter, I have found this post that describes how to handle it in ArcGIS. However, since I want to process lot of shapefiles (500MB) and raster data is huge in size, using such an approach is not an option.

  • Just as a comment: Hansen et al. already published country-wise forest tree-cover extent for the baseline 2000 and 2012. See the Supplementary Information sciencemag.org/content/342/6160/850/suppl/DC1 – Curlew Jun 23 '14 at 11:32
  • Thanks! I am using countries since then I can validate my results against Hansen's and hence check whether I am doing the right thing. In the end I will use e.g. administrative boundaries or other polygons. – Lukasz Tracewski Jun 23 '14 at 11:52

What you try to do, GIS people call "Zonal Statistics". There are many ways to do it, either in a Desktop Gis or "scriptwise".

In QGIS there is a Plug In called Zonal Statistics https://docs.qgis.org/2.6/en/docs/user_manual/plugins/plugins_zonal_statistics.html which you find in Rastermethods.

As you mentioned that you either want to do it in C++ or Python, I recommend the latter. Python has some very powerful Geo libraries. The one you are looking for is called rasterstats. https://github.com/perrygeo/python-rasterstats and/or http://pythonhosted.org/rasterstats/. On the github page (as linked above) you can see, that what you are trying to do, is very easy.

from rasterstats import zonal_stats
stats = zonal_stats("country_shapefile.shp", "forestcover.tif")
> ['count', 'min', 'max', 'mean']
[f['mean'] for f in stats]
> [756.6057470703125, 114.660084635416666]
| improve this answer | |
  • Thanks for your answer! Although the proposed solution certainly can nicely summarise rasters, it does not scale well beyond relatively small areas. I am upvoting this answer as a good one, but the right one turned out to be use of Google Earth Engine (which, btw, provides Python API). With the latter I managed to process planetary-scale datasets. – Lukasz Tracewski Dec 17 '16 at 7:16

If you want to use a programming language then you can use GDAL which is available for both Python and C++.

| improve this answer | |
  • 1
    Can you provide something a little more concrete, e.g. an example or reference to a specific API call that calculates such area? I have seen GDAL before, but from API reference nor tutorial I could not tell which function should do the trick. Once I know the name, I can work my way from there. – Lukasz Tracewski Jun 23 '14 at 11:46
  • Here's an example using a QGIS plugin, you can dig through the code if you prefer to not use the plugin: conservationecology.wordpress.com/qgis-plugins-and-scripts/… – camdenl Jun 23 '14 at 12:22

Having tried a few options, the best by far turned out to be using Google Earth Engine, a cloud computing platform for geospatial computations. It's currently beta, with no SLA (service-level agreement), but also freely available. That's exactly where Global Forest Change map was produced.

Great advantage of this solution is that it is scalable. The proposed Python solution works for small scales, but the moment you need to process area of size of countries and continents, it becomes a bottleneck.

In our research we calculated forest area for distribution maps of almost 12000 species, covering together roughly 1000 million square kilometres. Here's the complete code together with link to the article we published.

| improve this answer | |

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.