Online Trick

This blogpost shows how to find the dominant species and species counts within sampling polygons. The Species-layer that I'll use here is comprised of overlapping polygons which represent the distribution of several species. The Regions-layer represents areas of interest over which we would like to calculate some measures like species count, dominant species and area occupied by the dominant species.

Since QGIS now makes import/export and querying of spatial data easy, we can use the spatiaLite engine to join the intersection of both layers to the region table and then aggregate this intersections by applying max- and count-function on each region. We'll also keep the identity and the area-value of the species with the largest intersecting area.

For the presented example I'll use

Regions, which is a polygon layer with a areas of interest

Species, which is a polygon layer with overlapping features, representing species

Do the calculation in 2 easy steps:

Import the layers to a spatiaLite DB with the Import function of the plugin (example data: HERE)

Run the query. For later use you can load this table to QGIS or export with the plugin's export button.

SELECT   
  t.region AS region,
  t.species AS sp_dom,
  count(*) AS sp_number,
  max(t.sp_area) / 10000 AS sp_dom_area
  FROM ( 
      SELECT
          g.region AS region, s.species AS species,
          area(intersection(g.Geometry, s.Geometry)) AS sp_area
          FROM Regions AS g JOIN Sp_Distribution AS s 
          ON INTERSECTS(g.Geometry, s.Geometry)  
  ) AS t
GROUP BY t.region
ORDER BY t.region

Addendum:
If you wish to calculate any other diversity measures, like Diversity- or Heterogenity-Indices, you might just run the below query (which actually is the subquery from above) and feed the resulting table to any statistic-software!

The output table will contain region's IDs, each intersecting species and the intersection area.
The intersection area, which is the species' area per polygon, is the metric that would be used for the calculation of diversity / heterogenity measures, etc. of regions!

SELECT 
  g.region AS regID, 
  s.species AS sp,
  AREA(INTERSECTION(g.geometry, s.geometry)) AS sp_area
FROM Regions AS g JOIN Sp_Distribution AS s 
ON INTERSECTS(g.Geometry,s.Geometry)
ORDER BY regID, sp_area ASC

I tested this on

QGIS 2.6 Brighton

with the QspatiaLite Plugin installed

Following up my previous blogpost on using SpatiaLite for the calculation of diversity metrics from spatial data, I'll add this SQL-query which counts unique species-names from the intersection of species polygons and a circle-buffer around centroids of an input grid. The species number within the bufferarea are joined to a newly created grid. I use a subquery which grabs only those cells from the rectangular input grid, for which the condition that the buffer-area around the grid-cell's centroid covers the species unioned polygons at least to 80%.

Example data is HERE. You can use the shipped qml-stylefile for the newly generated grid. It labels three grid-cells with the species counts for illustration.

Import grid- and Sp_distr-layers with QspatiaLite Plugin

Run query and choose option "Create spatial table and load in QGIS", mind to set "geom" as geometry column

select 
    g1.PKUID as gID,
    count (distinct s.species) as sp_num_inbu, 
    g1.Geometry AS geom
from (
 select g.*
 from(select Gunion(geometry) as geom
           from Sp_distr) as u, grid as g
 where area(intersection(buffer(centroid(g.geometry), 500), u.geom)) > pow(500, 2)*pi()*0.8
) as g1 join Sp_distr as s on intersects(buffer(centroid( g1.Geometry), 500), s.Geometry)
group by gID

..I'm overwhelmed by the ever-growing loads of data that's made available via the web. I.e., IUCN collects and hosts spatial species data which is free for download. I'm itching to play with all this data... And, in the end there may arise some valuable outcome:

In the below examples I made a map for amphibian species richness - without much effort. I used a 5x5 km grid (constant areas, provided by Statistik Austria) and amphibian species ranges and intersected this data to yield species numbers per grid-cell. For aggregation of the tables from the polygons produced by the intersection I needed to switch from ArcMap to R because ESRI's "Summary Statistics" only cover MIN, MAX, MEAN, FIRST/LAST... and there is no easy way to apply custom functions. So, I just accessed the dbf that I wanted to aggregate with R and did the calculations there.. then reloaded the layer in ArcMap and was done..



Workflow-documentation:

(1) add downloaded data:

first add "L005000_LB.shp" (5000 m grid, zip-folder: l005000_lb.zip), this will make arcmap to choose this shp's GCS as custom,
which is "MGI".

then add "SPECIES1109_DBO_VIEW_AMPHIBIANS" (distribution ranges as polygons for each species, multipart, in case a species' areal is fragmented, zipped gdb-folder: ALL_AMPHIBIANS_Oct2010.zip).

when the warning about the GCS pops up, choose transformation from "GCS_WGS_1984" to "GCS_MGI"
using "MGI_to_WGS1984_3".


(2) intersect the two:

input featuters are "L005000_LB.shp" & "SPECIES1109_DBO_VIEW_AMPHIBIANS",
the resulting shp will be called "Amph_Sp_R5000_inters.shp".

choose the grid first, because this is used to defines the GCS for the output,
which we want to be "MGI"!

the grid quadrats are intersect with all species' polygons at a time, producing several "pieces"
of on grid cell, often more than one for a given species..
examine this by going to the attribute table and sorting by FID_L00500.

this is because ArcMap obviously intersects the to features at once, rather than one species
polygon after another - this would yield only one intersection per gridcell and species..

so, to eventually get species numbers we need to calculate unique species values per gridcell.
ArcMap's Summary Statistics provide a limited set of functions: SUM, MEAN, MAX, MIN, COUNT, FIRST, LAST
but we need no. of unique values -

we'll switch to R, and R will easily do this:

# for dbf import you'll need:
require(foreign)

# read data, mind to change to your download-path:
dat <- read.dbf("C:\\Gis_Daten\\Amphibians\\Amph_Sp_R5000_inters.dbf")

# check data:
str(dat)

# "FID_L00500" holds the identifier for each gridcell
# "BINOMIAL" holds the names of species
# the intersection produced polygons for each gridcell and species..
# so, aggregating  

# set up custom function that calculates no. of unique values:
x <- c(1,1,2,3,3,3,4)
n.unique <- function(x) {length(unique(x))}
n.unique(x)

# the following bit will give the species no. per gridcell:
grid.no.sp <- aggregate(BINOMIAL ~ FID_L00500, FUN = n.unique, data = dat)

# change column name "BINOMIAL" to "Sp.No.":
names(grid.no.sp)[2] <- "Sp.No."

# add ID column and re-order columns:
grid.no.sp$OID <- row.names(grid.no.sp)
grid.no.sp <- grid.no.sp[, c(3,1,2)]

# write file:
write.csv(grid.no.sp, file = "C:\\Gis_Daten\\Amphibians\\grid.no.sp.CSV",
          row.names = FALSE)

(3) join this CSV to your grid shp-file and your done! (!) i noticed that ArcMap often makes trouble when trying to use txt or csv files for joins & relates. if you experience the same: skip the write.csv() and alternatively add the below lines to your R script: (you'll need to remove the grid-layer from your ArcMap Session, otherwise it may not be writeable..)

# doing the join in R:
grid <- read.dbf("C:\\Gis_Daten\\Raster_Österr\\l005000_lb\\l005000_LB.dbf")
str(grid)

# bring into same order before joining:

grid <- grid[order(as.character(grid$NAME)), ]
grid.no.sp <- grid.no.sp[order(as.character(grid.no.sp$NAME)), ]

# join & check if they are properly aligned:
print(join <- data.frame(grid, grid.no.sp))
str(join)

write.dbf(join, file = "C:\\Gis_Daten\\Raster_Österr\\l005000_lb\\l005000_LB.dbf")

(4) then add the grid-file again and check the attribute table - sp.no. should be there. apply symbology and you're done!!

Note: Data is courtesy of STATISTIK AUSTRIA and IUCN

Here's a short follow up on some previous posting about the use of QspatiaLite for the aggregation of species distribution data. In this case the species data comes as a point layer. For each cell of a 1000 x 1000 m grid (1) the number of individuals per species, (2) the total number of individuals and (3) the number of different species should be calculated.


There is a layer with 10 different species (variabel name is "sp") across the whole extent with names "1", "2", "3", .. , "10" and a layer with the grid cells (variable name = "id") numbered consecutively, from 1 to 150.

In the attribute table of the below screenshot you see that I selected the grid cell with id=1 and the points (=species) within this cell. There are 8 individuals - "4", "5", "6" and "10" occure once, whereas "2" and "8" occure twice.

The query table in the screenshot is the result for (3).


For (1) you have to query for points/species within grid cells and group over grid cells and species and take the count from that aggregation

SELECT
 t.gID AS gID,
 t.sp AS Sp,
 count(*) AS NrIndSp
FROM (SELECT
  g.id AS gID, 
  s.sp AS sp
FROM grid AS g JOIN Sp_distr AS s 
ON within(s.Geometry, g.Geometry)
) as t GROUP BY t.gId, t.sp

For (2) you simple need to query for points/species within grid cells and aggregate over grid cells:

Select 
 t.gID,
 count(*) as NrInd
From (SELECT
  g.id AS gID, 
  s.sp AS sp
FROM grid AS g JOIN Sp_distr AS s 
ON within(s.Geometry, g.Geometry)
) as t 
GROUP BY t.gID 
ORDER BY t.gID

For (3) you'll first need to aggregate over grid cells and points/species, and then again aggregate over this query table by grid cells which will finally give you the distinct species!

SELECT 
 v.gID,
 count(*) AS SpNr
FROM (SELECT 
 t.gID,
 t.sp
FROM (SELECT
  g.id AS gID, 
  s.sp AS sp
FROM grid AS g JOIN Sp_distr AS s 
ON within(s.Geometry, g.Geometry)
) as t GROUP BY t.gId, t.sp
) as v GROUP BY v.gId