#################################### # 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)