### R Script for making a raster file of the Caucasus Ecoregion ### ############# 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("script3.txt") in the command window in R ##mark system time t_start<-Sys.time() #load several packages ##mark system time t_start<-Sys.time() #load several packages library(rgeos) library(rgdal) library(raster) library(spatial) ####################### raster preparation ################################ #... load raster file of North Caucasus from working directory ## for creating it see http://caucasus-spiders.info/r-spatial/raster-basics-1-the-north-caucasus/ north_raster <- raster("north_caucasus.asc") ##getting SRTM data of Georgia geo_raster = getData('alt', country='GEO', mask=TRUE) ##getting SRTM data of Armenia arm_raster = getData('alt', country='ARM', mask=TRUE) ##getting SRTM data of Azerbaijan aze_raster = getData('alt', country='AZE', mask=TRUE) #combine all rasters as one for the whole Caucasus Ecoregion cauc_raster <- mosaic(north_raster,geo_raster,aze_raster,arm_raster, fun=mean) ##################### shapefile preparation ################################# #read shapefile "north_borders" from working directory ## requires filesnorth_borders.shp, north_borders.dbf, north_borders.shx and north_borders.prj #these are not plotted in the map, its just for showing how to load the shapefile north_borders=readOGR(".", "north_borders") #save map as png png("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="Caucasus Ecoregion") ##plot cauc_raster in map plot(cauc_raster,add=TRUE,legend=FALSE) ############################################## plot legend ################################################# #set range for plotting the legend cauc_raster.range <- c(0, 6000) #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,2,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) #get all other borders for the region rus_border=getData('GADM', country='RUS', level=1) geo_border=getData('GADM', country='GEO', level=0) arm_border=getData('GADM', country='ARM', level=0) aze_border=getData('GADM', country='AZE', level=0) irn_border=getData('GADM', country='IRN', level=0) tur_border=getData('GADM', country='TUR', level=0) ukr_border=getData('GADM', country='UKR', level=0) kaz_border=getData('GADM', country='KAZ', level=0) ukr_border=getData('GADM', country='UKR', level=0) tkm_border=getData('GADM', country='TKM', level=0) # ... and plot them plot(rus_border,add=TRUE,fill=FALSE,) plot(geo_border,add=TRUE,fill=FALSE,) plot(aze_border,add=TRUE,fill=FALSE,) plot(arm_border,add=TRUE,fill=FALSE) plot(irn_border,add=TRUE) plot(tur_border,add=TRUE) plot(ukr_border,add=TRUE) plot(kaz_border,add=TRUE) plot(tkm_border,add=TRUE) #a box around the map box() #finalize the png file dev.off() ########################### plot map in R window ######################### #end all open R graphics graphics.off() #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(cauc_raster, add=TRUE) #plot(north_borders,add=TRUE) plot(rus_border,add=TRUE,fill=FALSE,) plot(geo_border,add=TRUE,fill=FALSE,) plot(aze_border,add=TRUE,fill=FALSE,) plot(arm_border,add=TRUE,fill=FALSE) plot(irn_border,add=TRUE) plot(tur_border,add=TRUE) plot(ukr_border,add=TRUE) plot(kaz_border,add=TRUE) plot(tkm_border,add=TRUE) ############################################################################ #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")