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