|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2
|
|
|
3
|
|
|
4 args <- commandArgs(trailingOnly = TRUE)
|
|
|
5
|
|
|
6 #Rscript clim_data.R 'worldclim' 'var' 'resolution' 'OutputFormat' #'FRA' #'prec1'
|
|
|
7 if (length(args)==0 | args[1]=="-h" | args[1]=="--help"){
|
|
|
8 print ("general script execution : Rscript clim_data.R \'worldclim\' \'var\' resolution \'OutputFormat\' #\'variable-to-plot\' #\'region_code\'")
|
|
|
9 print ("eg : Rscript clim_data.R \'worldclim\' \'prec\' 10 \'raster\' #\'prec1' #\'FRA\'")
|
|
|
10 q('no')
|
|
|
11 }
|
|
|
12
|
|
|
13
|
|
|
14 #Climatic variables dictionaries
|
|
|
15 months<-c("January","February","March","April","May","June","July","August","September","October","November","December")
|
|
|
16 bioclimatic_vars<-c("Annual Mean Temperature", "Mean Diurnal Range (Mean of monthly (max temp - min temp))", "Isothermality (BIO2/BIO7) (x 100)","Temperature Seasonality (standard deviation x100)","Max Temperature of Warmest Month","Min Temperature of Coldest Month","Temperature Annual Range (BIO5-BIO6)","Mean Temperature of Wettest Quarter","Mean Temperature of Driest Quarter","Mean Temperature of Warmest Quarter","Mean Temperature of Coldest Quarter","Annual Precipitation","Precipitation of Wettest Month","Precipitation of Driest Month","Precipitation Seasonality (Coefficient of Variation)","Precipitation of Wettest Quarter","Precipitation of Driest Quarter","Precipitation of Warmest Quarter","Precipitation of Coldest Quarter")
|
|
|
17
|
|
|
18 #Function to create custom plot title
|
|
|
19 get_plot_title<-function(usr_var,usr_var_to_plot){
|
|
|
20 match<-str_extract(usr_var_to_plot,"[0-9]+")
|
|
|
21 if(usr_var %in% c("prec","tmin","tmax")){
|
|
|
22 printable_var<-(months[as.integer(match)])
|
|
|
23 if(usr_var=="prec"){
|
|
|
24 printable_var<-paste(printable_var," precipitations (mm)",sep="")
|
|
|
25 }else if(usr_var=="tmin"){
|
|
|
26 printable_var<-paste(printable_var," minimum temperature (°C *10)",sep="")
|
|
|
27 }else{
|
|
|
28 printable_var<-paste(printable_var," maximum temperature (°C *10)",sep="")
|
|
|
29 }
|
|
|
30 }else if(usr_var=="bio"){
|
|
|
31 printable_var<-(bioclimatic_vars[as.integer(match)])
|
|
|
32 printable_var<-paste("Bioclimatic variable - ",printable_var,sep="")
|
|
|
33 }
|
|
|
34 title<-paste("Worldclim data - ",printable_var,".",sep="")
|
|
|
35 return(title)
|
|
|
36 }
|
|
|
37
|
|
|
38
|
|
|
39 #Call libraries
|
|
|
40 library('raster',quietly=TRUE)
|
|
|
41 library(sp,quietly = TRUE, warn.conflicts = FALSE)
|
|
|
42 library(ncdf4,quietly = TRUE, warn.conflicts = FALSE)
|
|
|
43 library(rgdal,quietly = TRUE, warn.conflicts = FALSE) #To save as geotif
|
|
|
44 library(stringr)
|
|
|
45
|
|
|
46
|
|
|
47 #Get args
|
|
|
48 usr_data=args[1]
|
|
|
49 usr_var=args[2]
|
|
|
50 usr_res=as.numeric(args[3])
|
|
|
51 usr_of=args[4]
|
|
|
52
|
|
|
53
|
|
|
54 # Retrieve 'var' data from WorldClim
|
|
|
55 global.var <- getData(usr_data, download = TRUE, var = usr_var, res = usr_res)
|
|
|
56
|
|
|
57 # Check if we actualy get some
|
|
|
58 if (length(global.var)==0){
|
|
|
59 cat("No data found.")
|
|
|
60 }else{
|
|
|
61 writeRaster(global.var, "output_writeRaster", format=usr_of,overwrite=TRUE)
|
|
|
62 final_msg<-paste("WorldClim data for ", usr_var, " at resolution ", usr_res, " in ", usr_of, " format\n", sep="")
|
|
|
63 cat(final_msg)
|
|
|
64 }
|
|
|
65
|
|
|
66
|
|
|
67
|
|
|
68
|
|
|
69
|
|
|
70
|
|
|
71 #################
|
|
|
72 ##Visualisation##
|
|
|
73 #################
|
|
|
74
|
|
|
75 #Get args
|
|
|
76 if(length(args[5])>=0 && length(args[6])>=0){
|
|
|
77 usr_var_to_plot=args[5]
|
|
|
78 usr_plot_region=args[6]
|
|
|
79 }else{q('no')}
|
|
|
80
|
|
|
81 list_region_mask<-c("FRA","DEU","GBR","ESP","ITA")
|
|
|
82
|
|
|
83 if(usr_plot_region %in% list_region_mask){
|
|
|
84 #Country mask
|
|
|
85 region <- getData("GADM",country=usr_plot_region,level=0)
|
|
|
86 region_mask <- mask(global.var, region)
|
|
|
87 region_var_to_plot_expression<-paste("region_mask$",usr_var_to_plot,sep="")
|
|
|
88 }else{ #All map and resize manualy
|
|
|
89 region_var_to_plot_expression<-paste("global.var$",usr_var_to_plot,sep="")
|
|
|
90 }
|
|
|
91
|
|
|
92
|
|
|
93 region_var_to_plot<-eval(parse(text=region_var_to_plot_expression))
|
|
|
94
|
|
|
95 #PLotmap
|
|
|
96 jpeg(file="worldclim_plot_usr_region.jpeg",bg="white")
|
|
|
97
|
|
|
98 title<-get_plot_title(usr_var,usr_var_to_plot)
|
|
|
99
|
|
|
100
|
|
|
101 if(usr_plot_region=="FRA"){
|
|
|
102 #FRA
|
|
|
103 plot(region_var_to_plot, xlim = c(-7, 12), ylim = c(40, 52), axes=TRUE,xlab="Longitude",ylab="Latitude",main=title)
|
|
|
104 }else if(usr_plot_region=="GBR"){
|
|
|
105 #GBR
|
|
|
106 plot(region_var_to_plot, xlim = c(-10, 5), ylim = c(46, 63), axes=TRUE,xlab="Longitude",ylab="Latitude",main=title)
|
|
|
107 }else if(usr_plot_region=="NA"){
|
|
|
108 #North America :
|
|
|
109 plot(region_ar_to_plot,xlim=c(-180,-50),ylim=c(10,75),xlab="Longitude",ylab="Latitude",main=title)
|
|
|
110 }else if(usr_plot_region=="EU"){
|
|
|
111 #Europe
|
|
|
112 plot(region_var_to_plot,xlim=c(-28,48),ylim=c(34,72),xlab="Longitude",ylab="Latitude",main=title)
|
|
|
113 }else if(usr_plot_region=="DEU"){
|
|
|
114 #DEU
|
|
|
115 plot(region_var_to_plot, xlim = c(5, 15), ylim = c(45, 57),axes=TRUE,xlab="Longitude",ylab="Latitude",main=title)
|
|
|
116 }else if(usr_plot_region=="ESP"){
|
|
|
117 #ESP
|
|
|
118 plot(region_var_to_plot, xlim = c(-10, 6), ylim = c(35, 45), axes=TRUE,xlab="Longitude",ylab="Latitude", main=title)
|
|
|
119 }else if(usr_plot_region=="ITA"){
|
|
|
120 #ITA
|
|
|
121 plot(region_var_to_plot, xlim = c(4, 20), ylim = c(35, 48), axes=TRUE,xlab="Longitude",ylab="Latitude", main=title)
|
|
|
122 }else if(usr_plot_region=="WM"){
|
|
|
123 #Worldmap
|
|
|
124 plot(region_var_to_plot,xlab="Longitude",ylab="Latitude",main=title)
|
|
|
125 }else if(usr_plot_region=="AUS"){
|
|
|
126 #AUS
|
|
|
127 plot(region_var_to_plot,xlim=c(110,155),ylim=c(-45,-10),xlab="Longitude",ylab="Latitude",axes=TRUE,main=title)
|
|
|
128 }else{
|
|
|
129 write("Error with country code.", stderr())
|
|
|
130 q('no')
|
|
|
131 }
|
|
|
132
|
|
|
133 garbage_output<-dev.off
|
|
|
134
|
|
|
135 #Exit
|
|
|
136 q('no')
|
