####################################
# Costa Rica Land Cover 		   #
# GlobCover 2008				   #
#  joewheatley.net Nov 2009        #
####################################

#library(sp)
library(rgdal)

# costa rica provinces

con <- url("http://gadm.org/data/rda/CRI_adm1.RData")
print(load(con))
close(con)

#shell("gdalwarp GLOBCOVER_200412_200606_V2.2_CentralAmerica_Reg.tif -te -86 8 -82 12 costaRica.tif")
costa <- readGDAL("costaRica.tif")


# puntarenas is slot 6. polygons 27 & 25 are mainland. islands e.g. cocos island are omitted

temp <- Polygons(list(Polygon(gadm@polygons[[6]]@Polygons[[27]]@coords),Polygon(gadm@polygons[[6]]@Polygons[[25]]@coords)),"puntarenas")
temp <- SpatialPolygons(list(temp),proj4string=CRS(proj4string(gadm)))
punt.sp <- SpatialPolygonsDataFrame(temp, data.frame(cbind(2,2), row.names=c("puntarenas")))  # puntarenas


puntarenas <- costa
puntarenas.overlay <- overlay(costa,punt.sp)  # 1 in interior of puntarenas polygons, 0 outside
puntarenas$band1 <- costa$band1*puntarenas.overlay

levs <- sort(unique(puntarenas$band1));  # land cover units present

# Globcover Regional Legend is in .csv file CentralAmericaLegend.csv

puntlegend <- read.csv("CentralAmerica_Legend.csv",header=T)
lev_values <- puntlegend$Value
lev_classes <- puntlegend$Label

present <- which(lev_values %in% levs)
lev.present <- lev_values[present]  # lev_values which are present in puntarenas
classes.present <- lev_classes[present]

# Globcover colors  - convert RGB to R colors according to GlobCover prescription

lev_cols <- c(rep(0,length(puntlegend$Red)))
for (i in 1: length(puntlegend$Red)){
lev_cols[i] <- rgb(puntlegend$Red[i],puntlegend$Green[i],puntlegend$Blue[i],maxColorValue=255)
}

cols.present <- lev_cols[present]

#legend

windows(40,40)
plot.new()
legend("left",as.vector(classes.present),fill=cols.present,cex=0.8,bty="n",box.col="white")

#map in traditional plot system

windows(60,40)
nf <- layout(matrix(c(1,2),1,2,byrow=T),c(4,1),c(4),T)
#layout.show(nf)

bks <- levs-0.5;
bks <- c(bks,max(levs)+5)  # breaks for image

op <- par(mar=c(3, 3, 3, 1))
plot(punt.sp,axes=T)
image(puntarenas,breaks=bks ,col=cols.present, add=T)
title(main="Puntarenas Province, Costa Rica")
plot(punt.sp,add=T,axes=T)
text(-84,8,"Data from GlobCover http://ionia1.esrin.esa.int/")
par(mar=c(3, 3, 3, 1))
barplot(table(puntarenas$band1*0.09),names.arg="",col=cols.present,las=2,horiz=T,axes=F,main="Area (km2)")
axis(1, at = c(0,20000,40000), labels=c("0k","20k","40k"))
par(op)