Mercurial > repos > mnhn65mo > get_worldclim_data

comparison clim_data.R @ 0:bdf7d9f2ddae draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author mnhn65mo
date Thu, 09 Aug 2018 05:01:04 -0400
parents
children 008cbcc1e3f6
comparison
equal deleted inserted replaced
-1:000000000000 0:bdf7d9f2ddae
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')