## R Script for making a raster file of the North Caucasus ####
############# and shapefiles of its country borders ###############
######## 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("script.R") in the command window in R

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

#load several packages
library(rgeos)
library(rgdal)
library(raster)
library(spatial)

####################### raster ################################

##... get SRTM data of Russia (comment out after download for faster compiling)
getData('alt', country='RUS', mask=TRUE)

#... and load raster file from working directory
rus_raster <- raster("RUS1_msk_alt.gri")

#crop rus_raster to Caucasus Ecoregion
e <- extent(37, 50.5, 38.5, 48)
rus_new <- crop(rus_raster, e)


##################### borders #################################
rus_border=getData('GADM', country='RUS', level=1)

rus_id<-rus_border[,c("NAME_1","ID_1")]

#list all region names
print(rus_id$NAME_1)

#read polygon for all regions in North Caucasus:
ady <- rus_border[rus_border$NAME_1 == "Adygey",]
kra <- rus_border[rus_border$NAME_1 == "Krasnodar",]
sta <- rus_border[rus_border$NAME_1 == "Stavropol'",]
kara <- rus_border[rus_border$NAME_1 == "Karachay-Cherkess",]
karb <- rus_border[rus_border$NAME_1 == "Kabardin-Balkar",]
ose <- rus_border[rus_border$NAME_1 == "North Ossetia",]
ing <- rus_border[rus_border$NAME_1 == "Ingush",]
che <- rus_border[rus_border$NAME_1 == "Chechnya",]
dag <- rus_border[rus_border$NAME_1 == "Dagestan",]

#merge borders
north_borders <- ady+kra+sta+kara+karb+ose+ing+che+dag


#crop rus_raster to caucasus
e <- extent(north_borders)
rus_new <- crop(rus_raster, e)

# now use the mask function to limit the rasters to the extend of the border polygons
kra_raster <- mask(rus_new, kra)
ady_raster <- mask(rus_new, ady)
sta_raster <- mask(rus_new, sta)
kara_raster <- mask(rus_new, kara)
karb_raster <- mask(rus_new, karb)
ose_raster <- mask(rus_new, ose)
ing_raster <- mask(rus_new, ing)
che_raster <- mask(rus_new, che)
dag_raster <- mask(rus_new, dag)

#... and make mosaic of these rasters
northcauc_raster <- mosaic(kra_raster, ady_raster,sta_raster,kara_raster,karb_raster,ose_raster,ing_raster,che_raster,dag_raster, fun=mean)

#write raster file of North Caucasus
writeRaster(northcauc_raster, filename="north_caucasus.asc", format="ascii",overwrite=TRUE)
#hdr(northcauc_raster, format="ENVI")

#write shapefile with borders in North Caucasus
wd <- getwd()
writeOGR(north_borders,wd, "north_borders", driver="ESRI Shapefile", overwrite=TRUE)


#uncomment for direct graphic output of map (with corrupted legend, sorry)
#plot(1,2,xlim=c(37,50.5),ylim=c(38.5,47))
#plot(northcauc_raster, add=TRUE)
#plot(north_borders,add=TRUE)

#or save map as png
png("north_caucasus.png",width=15, height=12, units = "cm", res=300)

plot(1,2,xlab="Longitude",ylab="Latitude",xlim=c(37,50.5),ylim=c(38.5,47),main="North Caucasus")

##plot northcauc_raster in map
plot(northcauc_raster,add=TRUE,legend=FALSE)

cauc_raster.range <- c(0, 6000)

############################################## plot legend #################################################

#plot "SRTM ..." below color bar legend
#mgp=c(axis label-axis position, axis label - tick label, axis position - axis line (usually 0))
plot(cauc_raster,add=TRUE,horizontal=TRUE, smallplot=c(.2, .4, .3, .315),legend.only=TRUE,
	axis.args=list(col.axis="white",col="white", mgp = c(3, -0.8, 0)),
	legend.args=list(text='(SRTM 90 m resolution)', side=1, font=1, line=.5, cex=0.3)
	)

#plot color bar legend
#mgp=c(axis label-axis position, axis label - tick label, axis posistion - axis line (usually 0))
#in legend.args=... line=Y bestimmt Y-Versatz des Textes
#tck bestimmt tick mark size
plot(cauc_raster,add=TRUE,horizontal=TRUE, smallplot=c(.2, .4, .3, .315),legend.only=TRUE,
	axis.args=list(
		at=seq(cauc_raster.range[1], cauc_raster.range[2], 1000),
                labels=seq(cauc_raster.range[1], cauc_raster.range[2], 1000), 
                cex.axis=0.4,
		mgp = c(0, -0.3, 0),
		tck=-0.5),
	legend.args=list(text='Elevation (m)', side=1, font=2, line=.2, cex=0.6)
	)

#plot region borders in map
plot(north_borders,add=TRUE)

#a box around the map
box()

#finalize the png file
dev.off()

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


#mark system time again, subtract from t_start and print difference
t_end<-Sys.time()
t_dur=t_end-t_start
print(t_dur)

print("end of file")