First get the data
library(vegan)
## Loading required package: permute
## Attaching package: 'permute'
## The following object(s) are masked from 'package:stats':
##
## nobs
## Warning: replacing previous import 'nobs' when loading 'permute'
## This is vegan 2.0-0
library(rgdal)
## Loading required package: sp
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.8.0, released 2011/01/12 Path to GDAL shared
## files: /usr/share/gdal/1.8 Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September
## 2009, [PJ_VERSION: 470] Path to PROJ.4 shared files: (autodetected)
library(foreign)
setwd("/home/duncan/Dropbox/Public/Analyses/REDDEAM")
d <- readOGR("shapes", "Trees_cacti_agave")
## OGR data source with driver: ESRI Shapefile
## Source: "shapes", layer: "Trees_cacti_agave"
## with 161707 features and 6 fields
## Feature type: wkbPoint with 2 dimensions
eco <- readOGR("shapes", "ecoregions")
## OGR data source with driver: ESRI Shapefile
## Source: "shapes", layer: "ecoregions"
## with 694 features and 17 fields
## Feature type: wkbPolygon with 2 dimensions
Run an overlay to set up a gid
olay <- over(d, eco)
dd <- data.frame(olay, d)
tb <- table(dd$name)
selected <- names(tb[tb > 3])
dd <- subset(dd, name %in% selected)
Now produce a matrix
tb <- table(dd$cveecon4, dd$name)
tb <- tb[apply(tb, 1, sum) > 3, ]
tb <- tb > 0
dismat <- betadiver(tb, "sim")
clus <- hclust(dismat, method = "average")
plot(clus)
l1a <- cutree(clus, 8)
l1a <- data.frame(cveecon4 = names(l1a), l1a)
writeOGR(eco, "shapes", "testeco", "ESRI Shapefile")
## Error: layer exists, use a new layer name
a <- read.dbf("shapes/testeco.dbf")
d2 <- merge(l1a, a, all.y = T, sort = FALSE)
d2 <- d2[order(d2$cat), ]
write.dbf(d2, "shapes/testeco.dbf")
