31 Oct 2011

Using IUCN-Data, ArcMap 9.3 and R to Map Species Diversity

I'm overwhelmed by the ever-growing loads of good data that's made available via the web. For instance 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 example I made a map for amphibian species richness without much effort. I used a 5x5 km grid (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, that were produced by 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

9 comments :

  1. Great (and simple) post. Very useful to me.

    Thanks!

    ReplyDelete
  2. Really helpfull,

    Simple and easy to follow.

    Thank you, this is exactly what i was looking for for hours and hours

    Awal,
    Indonesian

    ReplyDelete
  3. I followed this example and it was great, thanks :) However, do you have any links to other country shapefiles which are already gridded, or if not, how would I do this myself manually?

    ReplyDelete
    Replies
    1. Sorry, I don't have other resources at hand - you'll just need to search the web for yourself! When you found suitable species data, you can make your own grid and then follow the example. Making a custom grid is really simple in QGIS, but I guess it also simple in ArcGIS - you'll find plenty of tutorials on the web!

      Delete
  4. Hi
    This post is pretty helpful.Thanks for that.But I was wondering if one was able to do a similar thing if species richness was available in the form of points rather than ranges-polygons.Then how will one count each unique species per grid?

    ReplyDelete
    Replies
    1. I think you meant "species occurrence", rather than "species richness", I guess. I just recently wrote a blog post on how to aggregate species occurrences from point data in QGIS with QspatiaLite: LINK.
      To get you started with QspatiaLite some of my older postings covering QspatiaLite may be of use, supposing you never have dealt with it before!

      Delete
  5. Oh yes species occurrence not richness.Thank you very much.This is extremely helpful

    ReplyDelete