annotate clim_data.R @ 1:008cbcc1e3f6 draft

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