Download and read ASCII csv of targers and shapefile of planning units

 

    require(rgdal)
    require(sp)
    
    pu.data <- read.csv("http://s3.amazonaws.com/nrcs-sgi/EvansModels/HaitiAllTargets.csv")  
      row.names(pu.data) <- pu.data$ID
    download.file("http://s3.amazonaws.com/nrcs-sgi/EvansModels/PU_POLYS.shp",
                  destfile="PU_POLYS.shp", mode="wb")
    download.file("http://s3.amazonaws.com/nrcs-sgi/EvansModels/PU_POLYS.shx",
                  destfile="PU_POLYS.shx", mode="wb")
    download.file("http://s3.amazonaws.com/nrcs-sgi/EvansModels/PU_POLYS.dbf",
                  destfile="PU_POLYS.dbf", mode="wb")
    download.file("http://s3.amazonaws.com/nrcs-sgi/EvansModels/PU_POLYS.prj",
                  destfile="PU_POLYS.prj")
    rarity.sp <- readOGR(getwd(), "PU_POLYS")
    
Function to calculate H' rarity measure. The expected of the metric are:
      H' > 1  Over dispersed
      H' = 1 In equilibrium 
      H' < 1 Under represented (rare)

    H <- function (x, t=1, MARGIN=2) {
       H <- as.data.frame(array(0, dim=c( dim(x)[1], 0 )))
        rownames(H) <- rownames(x)
         tcounts <- apply(x, MARGIN=1, function(x){ length(x[x > 0]) } )
         total <- apply(x, MARGIN=2, sum) # COLUMN TOTALS
            for (n in 1:dim(x)[2] ) {
                s = total[n]
                  p <- (x[,n] / s)
                  a <- length(x[x[,n] > 0 ,]) / length(x[,n])
                r <- sqrt(ifelse(p >= 0, p, 0)) / sqrt(a)   
              H <- cbind(H, r)
            }
          names(H) <- names(x)
            hcounts <- apply(H, MARGIN=1, function(x){ length(x[x >= t]) } )  
            H <- cbind(ID=as.numeric(as.character(rownames(H))), TCOUNT=tcounts, H,  
                       H=(rowSums(H)), maxH=apply(H, MARGIN=1, max), RARITY=hcounts )
         return ( cbind(H, sH=(H$H / max(H$H)), sHmax=(H$maxH / max(H$maxH)) ) )
      }
       

Calculate rarity of conservation targets in protofolio

 

    rarity <- H(pu.data[,2:ncol(pu.data)])
    
Join rarity to planning units shapefile

    rarity.sp@data=data.frame(rarity.sp@data, rarity[match(rarity.sp@data[,"UNIT_ID"],
                              rarity[,"ID"]),])  
      rarity.sp@data <- rarity.sp@data[,-c(1,3)]
        rarity.sp@data <- na.omit(rarity.sp@data)    
        
Plot distribution and counts of H' <= 1   

    plot(density(rarity[rarity$H <= 1 ,]$H), main="DISTRIBUTION OF H <= 1")
      paste("H <= 1", dim(rarity[rarity$H <= 1 ,])[1], sep=" - N=")
   

Plot rarity of planning units

 

    require(RColorBrewer)
    require(classInt)   
     plotvar <- rarity.sp@data[,"H"]
       nclr <- 6
         plotclr <- rev(brewer.pal(nclr,"Spectral"))
           cuts <- classIntervals(plotvar, n=nclr, style="kmeans")  
             colcode <- findColours(cuts, plotclr)
        plot(rarity.sp, col=colcode, border=NA)
          legend("topleft", legend=rev(c("Very Rare","Rare","Moderatly Rare",
                 "Somewhat Common","Common","Over Dispersed")),
                 fill=attr(colcode, "palette"), cex=0.6, bty="n")