## R Script for area calculations on a raster object of Georgia ###
######## script tested on Ubuntu 12.04 (precise) (32-Bit) #########
############### R version 3.1.0 / JGR version 1.7-16 ##############
############################ May 2014 #############################

#set working directory
#setwd("/home/...")
#type source("area_script.txt") in the command window in R

##mark system time
t_start<-Sys.time()


#load several packages
library(raster)
library(rgdal)
library(geosphere)

######################### calculating the area of a polygon ####################################################

geo_border=getData('GADM', country='GEO', level=0)

#calculate area [m^2] of the polygon
sqm<-areaPolygon(geo_border)

#convert sqm to km^2
sqkm<-sqm/1000000

#print area of Georgia according to shapefile, rounded to one digit
print(paste("Area of Georgia (shapefile):",round(sqkm, digits=1),"km2"))

########################## calculating the area of a raster object ##############################################

##getting SRTM data of Georgia
geo_raster = getData('alt', country='GEO', mask=TRUE)

#get sizes of all cells in raster [km2]
cell_size<-area(geo_raster, na.rm=TRUE, weights=FALSE)

#delete NAs from vector of all raster cells
##NAs lie outside of the rastered region, can thus be omitted
cell_size<-cell_size[!is.na(cell_size)]

#compute area [km2] of all cells in geo_raster
raster_area<-length(cell_size)*median(cell_size)

#print area of Georgia according to raster object
print(paste("Area of Georgia (raster):",round(raster_area, digits=1),"km2"))

##plot
graphics.off()
X11(width=8, height=5)
plot(geo_raster,main="Georgia")

plot(geo_border,add=TRUE)

################################## visualize NAs in raster #################################################

geo_rasterNA <- geo_raster
geo_rasterNA[is.na(geo_rasterNA)] <- 5000

X11(width=8, height=5)

plot(geo_rasterNA,main="Georgia",sub="raster NAs in green")

############################### subset rasters by raster value ################################################

#below-zero areas
geo_raster0 <- geo_raster
geo_raster0[geo_raster0 >-1] <- NA

#calculate area of regions under 0 m asl
#get sizes of all cells under 0 m
cell_size<-area(geo_raster0, na.rm=TRUE, weights=FALSE)

#delete NAs from vector of all raster cells
##NAs lie outside of the rastered region, can thus be omitted
cell_size<-cell_size[!is.na(cell_size)]

#compute area [km2] of all cells in geo_raster
underzero_area<-length(cell_size)*median(cell_size)

#print area below sea level to raster object
print(paste("Area of regions under below sea level:",round(underzero_area, digits=1),"km2"))

X11(width=8, height=5)
plot(geo_raster0,width=15, height=10,main="Georgia areas below sea level",sub=paste("below sea level, area =",round(underzero_area, digits=1),"km2"))
plot(geo_border,add=TRUE)


#######################################################
##lowland zone 0-999 m
geo_raster1000 <- geo_raster

geo_raster1000[geo_raster1000 <=-1] <- NA
geo_raster1000[geo_raster1000 >999] <- NA

#calculate area of regions under 0 m asl
#get sizes of all cells under 0 m
cell_size<-area(geo_raster1000, na.rm=TRUE, weights=FALSE)

#delete NAs from vector of all raster cells
##NAs lie outside of the rastered region, can thus be omitted
cell_size<-cell_size[!is.na(cell_size)]

#compute area [km2] of all cells in geo_raster
lowland_area<-length(cell_size)*median(cell_size)

#print area of Georgia according to raster object
print(paste("Area of lowland regions (0-999 m):",round(lowland_area, digits=1),"km2"))

#plot lowland zone
X11(width=8, height=5)
plot(geo_raster1000,width=15, height=10,main="Georgia lowland areas",sub=paste("0-999 meters asl, area =",round(lowland_area, digits=1),"km2"))
plot(geo_border,add=TRUE)

#######################################################
##montane zone 1000-2400 m
geo_raster2400 <- geo_raster

geo_raster2400[geo_raster2400 <=1000] <- NA
geo_raster2400[geo_raster2400 >2399] <- NA

#calculate area of regions under 0 m asl
#get sizes of all cells under 0 m
cell_size<-area(geo_raster2400, na.rm=TRUE, weights=FALSE)

#delete NAs from vector of all raster cells
##NAs lie outside of the rastered region, can thus be omitted
cell_size<-cell_size[!is.na(cell_size)]

#compute area [km2] of all cells in geo_raster
montane_area<-length(cell_size)*median(cell_size)

#print area of Georgia according to raster object
print(paste("Area of montane regions (1000-2399 m):",round(montane_area, digits=1),"km2"))

#plot montane zone
X11(width=8, height=5)
plot(geo_raster2400,width=15, height=10,main="Georgia montane areas",sub=paste("1000-3999 meters asl, area =",round(montane_area, digits=1),"km2"))
plot(geo_border,add=TRUE)

#######################################################
##alpine zone 2400-3999
geo_raster4000 <- geo_raster

geo_raster4000[geo_raster4000 <=2400] <- NA
geo_raster4000[geo_raster4000 >3999] <- NA

#calculate area of regions under 0 m asl
#get sizes of all cells under 0 m
cell_size<-area(geo_raster4000, na.rm=TRUE, weights=FALSE)

#delete NAs from vector of all raster cells
##NAs lie outside of the rastered region, can thus be omitted
cell_size<-cell_size[!is.na(cell_size)]

#compute area [km2] of all cells in geo_raster
alpine_area<-length(cell_size)*median(cell_size)

#print area of Georgia according to raster object
print(paste("Area of alpine regions (2400-3999 m):",round(alpine_area, digits=1),"km2"))

#plot alpine zone
X11(width=8, height=5)
plot(geo_raster4000,main="Georgia alpine areas",sub=paste("2400-3999 meters asl, area =",round(alpine_area, digits=1),"km2"))
plot(geo_border,add=TRUE)

#######################################################
##nival zone 4000+
geo_raster6000 <- geo_raster

geo_raster6000[geo_raster6000 <=4000] <- NA

#calculate area of regions under 0 m asl
#get sizes of all cells under 0 m
cell_size<-area(geo_raster6000, na.rm=TRUE, weights=FALSE)

#delete NAs from vector of all raster cells
##NAs lie outside of the rastered region, can thus be omitted
cell_size<-cell_size[!is.na(cell_size)]

#compute area [km2] of all cells in geo_raster
nival_area<-length(cell_size)*median(cell_size)

#print area of Georgia according to raster object
print(paste("Area of nival regions (4000+ m):",round(nival_area, digits=1),"km2"))

#plot lowland zone
X11(width=8, height=5)
plot(geo_raster6000,ylim=c(38.5,47),main="Georgia nival areas",sub=paste("4000+ meters asl, area =",round(nival_area, digits=1),"km2"))
plot(geo_border,add=TRUE)

#######################################################

#print area of Georgia as sum of all altitudinal zones
geo_area_comb=underzero_area+lowland_area+montane_area+alpine_area+nival_area
print(paste("Area of Georgia (combined alt. zones):",round(geo_area_comb, digits=1),"km2"))

######################################################################################################################
#Systemzeit nehmen und Berechnungsdauer ausrechnen, drucken
t_end<-Sys.time()
t_dur=t_end-t_start
print(t_dur)

######################################################################################################################################

print("end of file")