LecoS update – Landcover Polygon overlay

Hi folks,

it’s already late at night and after nearly 4 hours of confused crap coding i finally managed to implement a new amazing feature in my QGIS plugin LecoS. Up to now it was only possible to calculate the landcover for a single raster. However most people usually want to have the values for their respective plots, points of interests, sites, … . For this reason i coded the new Polygon overlay tool, which could extract raster values and save them directly to the vector layers attribute table.

Here is an example of the new features with some of the data mentioned in my previous post about the horrific Herakles farm palm oil project directly in the vicinity of some of the last and oldest protected rainforests belonging to a biodiversity hotspot, the Guinean Forests of West Africa.

The green color is classified as mature forest with >= 60% tree cover

The green color is classified as mature forest with >= 60% tree cover

Quite often during the campaign emerged the question how much of the protected area is still covered by mature forest trees. I want to resolve this question (in a admittedly very imprecise way) and downloaded a classified raster from the CARPE project, which displays the forest cover in the region. I furthermore added the shapes of the 4 largest protected areas in the vicinity. Now i open the new Batch polygon overlay tool in the LecoS menu (Menubar -> raster) and choose the options as displayed below. It now calculates the total proportion of the first class (mature tree cover >= 60%) for each vector feature (national park). Results are added directly to the shapes attribute table.lecosBatchClass

Although my plugin was originally designed just to analyse classified raster images, i furthermore implemented some basic methods to extract underlying raster values to polygons. Here is a simple example with the previous shapes:

SRTM DEM for the vicinity of Korup NP.

SRTM DEM for the vicinity of Korup NP.

I’ve downloaded a SRTM Digital Elevation Model for the vicinity of Korup National park. Like before i added the geometries of all major protected areas in the area. Now i simply want to extract the mean Height from the SRTM DEM for each protected Area.

This could be achieved by starting the “Landcover polygon overlay tool”, select the raster and vector layer (beforehand make sure that they have the same spatial projection) and then choose NO as answer in the central combobox. Choose “Mean of Raster values” and click on OK.

The Overlay tool will directly add the calculated values to the attribute table for each vector feature. However QGis doesn’t update the attribute table for some reason and therefore you have to load in your shapefile to the QGis layer browser again. Don’t worry, the values are most likely being saved in the attribute table. The column name is always a combination of the ras_ + the used method, ergo LCprop for Landcover Proportion (to keep it short).

For the simple example above the results look like this:


As you can see in nearly all PA is a very high percentage of over 60% mature tree cover

On the GIS-Stackexchange forum are dozens of questions ([1],[2],[3],[4]) regarding this topic and I hope that this update will come in handy to calculate something like the sum or the mean of all the values from a underlying raster shape per polygon feature. If you have a point-layer you should use the very good Point-Sampling Plugin instead, which is also available in QGis plugin downloader.

Good night 🙂


Tags: , , , , ,

About Martin Jung

PhD researcher at the University of Sussex. Interested in nature conservation, ecology and biodiversity as well as statistics, GIS and 'big data'

12 responses to “LecoS update – Landcover Polygon overlay”

  1. Maxim Dubinin says :

    Regarding second part, did you reinvent Zonal Statistics? There are two plugins/flavors of it already in QGIS 😉

    • Curlew says :

      You’re right and i remember trying it once. I just tried it again to accomplish the second part with zonal statistics and it results in a QGis crash (1.9 Master). Zonal stats simply aborts at 66% and returns this non helpful error ( http://i.imgur.com/dUVnj.png ).
      I implemented the function for non-classified rasters because it wasn’t much work to add it and it could work in many operations where other tools fail, because it uses the powerful numpy and scipy libraries as backbone instead of QGIS functions. Nevertheless i believe its always good to have multiple options (i heavily depend on R and SAGA for instance, because there are often some (known) bugs in the software), while using open-source GIS software.

  2. Pedro says :


    Great work!

    I’ve been doing some tests and the results for the mean height from a DEM, with LecoS and Zonal statistics are slightly different. It’s normal?

    It would be interesting to offer additional statistics such as minimum, maximum, standard deviation, quartiles, etc., of Raster values.

    Thank you very much!

    Best regars!

    • Curlew says :

      What scale of difference are we talking about? Most likely those differences come up sometimes, because the method isn’t always 100% exact. The whole batch overlay tool works by clipping the raster shapes in a loop for each polygon feature. Sometimes it results in slightly bigger or smaller raster subsets and this could be the source of the differences.
      Anyway, thanks for the suggestions. Maybe i’ll implement them when i’ve more time to code.

      Nja, damn it. There goes my quiet Sunday. I uploaded a new version, which adds the option to calculate min/max/med/std/lowerQ/upperQ
      Have fun 🙂

      • Pedro says :


        Sorry to have contributed to your little rested Sunday! 🙂
        Thank you very much for this new update!!

        I leave just a few suggestions, but not for this Sunday! 🙂

        Instead of using radio buttons, you could consider using check boxes, to calculate various statistics simultaneously.

        The magnitude of the difference is not too large. On 27 polygons, it varies from 0.008 to 2.49. The mean difference is ~0.6 and the median is ~0.41. These differences are not directly proportional to the area of the polygons, that is, the biggest differences do not arise associated with larger polygons.

        Thanks again for you great work!!


  3. Janis Kloks says :


    Firstly thank you for the great Tools!

    But i have an error here that I don’t understand what it means:

    TypeError: coordinate list must contain at least 2 coordinates

    In my Attribute Table there are already 2 coordinates (X;Y) and the shapefile contains 1231 polygons each with an ID (some of ID’s are the same)

    and with another shapefile it says:

    shape mismatch: objects cannot be broadcast to a single shape

    Any helpful thoughts? Or should i ask this somewhere else?

    • Curlew says :


      i would really like to help you, but i can’t do that so well without having a look at the data myself. All i can do is interpret the errors from remote (you can imagine that this is quite hard) and make suggestions.

      In general: All of your files need to be in the same spatial projection. If they are not, re-project them (rightclick -> save as or use gdaltools). This is a crucial point and many errors are resolved after that. Your shape mismatch Error indicates that (1) either some of your polygons are outside the raster file, which is not allowed or (2) they are not in the same projection or (3) there are some weird features that overlay each other, incorrect polygon node placements and other gis-crap (i kinda got a similar error from a sample file before and didn’t found a solution for that). Try to simplify it somehow.

      The coordinates don’t have to be inside the attribute table (by the way, i wonder how you do that with polygons. LecoS accepts only polygons and not point-locations. Use the qgis Point-Sampling plugin instead).

      Try merging some of your features together or make a test-run with a random quickly drawn polygon on your raster.

      Did you also tried what you want to do with other tools or software (SAGA zonal statistics or the QGis-plugin Zonal statistics?). The “Grid Statistics for polygons” tool coming with SAGA in my Sextante toolbox is what i usually use for this task.
      Maybe there is some issues with your shapefile in general (double id´s for different features are usually not very good 😉 )


      • Janis Kloks says :

        thank’s for the answer! Yes I understand, but the file size is quite big, I can tell what data it is – Raster is CLC2006 land cover data and vector layer is like Admin 1 – States, Provinces from NaturalEarth.

        Regarding the errors, my projections are set right, because I did make a test of a smaller polygon already, and it worked perfectly. So it’s probably the administrative vector layer that has areas outside of the Raster layer as you said.

        To get the coordinates into Attribute Table from a vector layer you can define two columns in excel as an example, where they contain lat/long coordinates(I don’t remember exactly how it was done, but you have to use some simple excel formulas power 🙂 , save it as .csv and then in Qgis use – Layer – Add Delimited Text Layer

        What I want to do is get Raster statistics for each class in every polygon, so that I can say how many km2 of different land cover type i have in every administrative unit.

        I corrected that with ID’s, you’re right!

      • Janis Kloks says :


        it seems to work now, just made the Vector layer into separate polygon layer – SAGA – Polygon parts to separate polygons, but it takes too long to calculate 1 Class.

        However I had success exporting statistics for all the classes just for a single Raster with no polygon overlay and it was very fast, but i do not know how to automate the process for many Raster files at once, but then I would have to split my Raster layer into many separate layers (admin divisions) at first place

        I know that there are GME Tools that can do it, but it’s not yet opensource, it uses ArcMap..


        • Curlew says :

          It always takes a while if you’re using grids and shapes with large dimensions. Try the SAGA function “Grid statistics for Polygons”. Also consider using R or python for this task (if you want a scripted solution). Maybe in the future i will code Sextante support for LecoS, than you make a batch file computation with the function or use the sextante modeler.

  4. Janis Kloks says :

    Thanks, the Saga function solution worked, I’m separating now the Raster, otherwise it sums up all the pixels, eliminating the classes. Unfortunately I dont have enough knowledge for scripting solution, but the sextante support sounds good!

Sussex Research Hive

Supporting the research community at the University of Sussex

Small Pond Science

Research, teaching, and mentorship in the sciences

Landscape Ecology 2.0

intersecting landscape ecology, open science, and R


The Research Blog of IIASA

Amy Whitehead's Research

the ecological musings of a conservation biologist

Michael McCarthy's Research

School of BioSciences, The University of Melbourne

The Rostrum

science, statistics, policy and more


Environmental Change - Understand, Predict, Adapt

Dynamic Ecology

Multa novit vulpes


METeorological Visualisation Utilities using R for Science and Teaching

A Birder´s Blog

"Everybody loves what they know"

Trust Me, I'm a Geographer

Using Technology to Explore Our World

Duncan Golicher's weblog

Research, scripts and life in Chiapas

Daniel J. Hocking

Ecology, conservation biology, & statistical modeling for a changing world

%d bloggers like this: