1 |
30f84063
|
Benoit Parmentier
|
##################################### METHODS COMPARISON part 5 ##########################################
|
2 |
|
|
#################################### Spatial Analysis ############################################
|
3 |
|
|
#This script utilizes the R ojbects created during the interpolation phase. #
|
4 |
e9998467
|
Benoit Parmentier
|
#This scripts focuses on a detailed study of differences in the predictions of CAI_kr and FUsion_Kr
|
5 |
|
|
#Differences are examined through:
|
6 |
|
|
#1) per land cover classes
|
7 |
|
|
#2) per elevation classes
|
8 |
|
|
#3) through spiatial transects
|
9 |
|
|
#
|
10 |
|
|
#Note this code is for exploratory analyses so some sections are not succinct and
|
11 |
|
|
#can be improve for repeatability and clarity.
|
12 |
30f84063
|
Benoit Parmentier
|
#AUTHOR: Benoit Parmentier #
|
13 |
e9998467
|
Benoit Parmentier
|
#DATE: 12/04/2012 #
|
14 |
30f84063
|
Benoit Parmentier
|
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#491 -- #
|
15 |
|
|
###################################################################################################
|
16 |
|
|
|
17 |
|
|
###Loading R library and packages
|
18 |
|
|
library(gtools) # loading some useful tools such as mixedsort
|
19 |
|
|
library(mgcv) # GAM package by Wood 2006 (version 2012)
|
20 |
|
|
library(sp) # Spatial pacakge with class definition by Bivand et al. 2008
|
21 |
|
|
library(spdep) # Spatial package with methods and spatial stat. by Bivand et al. 2012
|
22 |
|
|
library(rgdal) # GDAL wrapper for R, spatial utilities (Keitt et al. 2012)
|
23 |
|
|
library(gstat) # Kriging and co-kriging by Pebesma et al. 2004
|
24 |
|
|
library(automap) # Automated Kriging based on gstat module by Hiemstra et al. 2008
|
25 |
|
|
library(spgwr)
|
26 |
|
|
library(gpclib)
|
27 |
|
|
library(maptools)
|
28 |
|
|
library(graphics)
|
29 |
|
|
library(parallel) # Urbanek S. and Ripley B., package for multi cores & parralel processing
|
30 |
|
|
library(raster)
|
31 |
|
|
library(rasterVis)
|
32 |
|
|
library(plotrix) #Draw circle on graph
|
33 |
|
|
library(reshape)
|
34 |
|
|
######### Functions used in the script
|
35 |
|
|
#loading R objects that might have similar names
|
36 |
|
|
load_obj <- function(f)
|
37 |
|
|
{
|
38 |
|
|
env <- new.env()
|
39 |
|
|
nm <- load(f, env)[1]
|
40 |
|
|
env[[nm]]
|
41 |
|
|
}
|
42 |
|
|
|
43 |
|
|
plot_transect<-function(list_trans,r_stack,title_plot,disp=TRUE){
|
44 |
|
|
#This function creates plot of transects for stack of raster images.
|
45 |
|
|
#The parameters are:
|
46 |
|
|
#list_trans: list of files containing the transects lines in shapefile format
|
47 |
|
|
#r_stack: raster stack of files
|
48 |
|
|
#title_plot: plot title
|
49 |
|
|
#disp: dispaly and save from X11 if TRUE
|
50 |
|
|
nb<-length(list_trans)
|
51 |
|
|
t_col<-rainbow(nb)
|
52 |
|
|
list_trans_data<-vector("list",nb)
|
53 |
|
|
for (i in 1:nb){
|
54 |
|
|
trans_file<-list_trans[[i]][1]
|
55 |
|
|
filename<-sub(".shp","",trans_file) #Removing the extension from file.
|
56 |
|
|
transect<-readOGR(".", filename) #reading shapefile
|
57 |
|
|
trans_data<-extract(r_stack, transect)
|
58 |
|
|
if (disp==FALSE){
|
59 |
|
|
png(file=paste(list_trans[[i]]),".png",sep="")
|
60 |
|
|
}
|
61 |
|
|
for (k in 1:ncol(trans_data[[1]])){
|
62 |
|
|
y<-trans_data[[1]][,k]
|
63 |
|
|
x<-1:length(y)
|
64 |
|
|
if (k!=1){
|
65 |
|
|
lines(x,y,col=t_col[k])
|
66 |
|
|
}
|
67 |
|
|
if (k==1){
|
68 |
|
|
plot(x,y,type="l",xlab="Position index", ylab="temperature",col=rainbow(k))
|
69 |
|
|
}
|
70 |
|
|
}
|
71 |
|
|
title(title_plot[i])
|
72 |
|
|
legend("topright",legend=layerNames(r_stack),
|
73 |
|
|
cex=1.2, col=t_col,
|
74 |
|
|
lty=1)
|
75 |
|
|
|
76 |
|
|
if (disp==TRUE){
|
77 |
|
|
savePlot(file=paste(list_trans[[i]][2],".png",sep=""),type="png")
|
78 |
|
|
}
|
79 |
|
|
if (disp==FALSE){
|
80 |
|
|
dev.off()
|
81 |
|
|
}
|
82 |
|
|
list_trans_data[[i]]<-trans_data
|
83 |
|
|
}
|
84 |
|
|
names(list_trans_data)<-names(list_trans)
|
85 |
|
|
return(list_trans_data)
|
86 |
|
|
}
|
87 |
|
|
|
88 |
|
|
plot_transect_m<-function(list_trans,r_stack,title_plot,disp=TRUE,m_layers){
|
89 |
|
|
#This function creates plot of transects for stack of raster images.
|
90 |
|
|
#Arguments:
|
91 |
|
|
#list_trans: list of files containing the transects lines in shapefile format
|
92 |
|
|
#r_stack: raster stack containing the information to extect
|
93 |
|
|
#title_plot: plot title
|
94 |
|
|
#disp: display and save from X11 if TRUE or plot to png file if FALSE
|
95 |
|
|
#m_layers: index for layerers containing alternate units to be drawned on a differnt scale
|
96 |
|
|
#RETURN:
|
97 |
|
|
#list containing transect information
|
98 |
|
|
|
99 |
|
|
nb<-length(list_trans)
|
100 |
|
|
t_col<-rainbow(nb)
|
101 |
|
|
list_trans_data<-vector("list",nb)
|
102 |
|
|
|
103 |
|
|
#For scale 1
|
104 |
|
|
for (i in 1:nb){
|
105 |
|
|
trans_file<-list_trans[[i]][1]
|
106 |
|
|
filename<-sub(".shp","",trans_file) #Removing the extension from file.
|
107 |
|
|
transect<-readOGR(".", filename) #reading shapefile
|
108 |
|
|
trans_data<-extract(r_stack, transect)
|
109 |
|
|
if (disp==FALSE){
|
110 |
|
|
png(file=paste(list_trans[[i]]),".png",sep="")
|
111 |
|
|
}
|
112 |
|
|
#Plot layer values for specific transect
|
113 |
|
|
for (k in 1:ncol(trans_data[[1]])){
|
114 |
|
|
y<-trans_data[[1]][,k]
|
115 |
|
|
x<-1:length(y)
|
116 |
|
|
m<-match(k,m_layers)
|
117 |
|
|
|
118 |
|
|
if (k==1 & is.na(m)){
|
119 |
|
|
plot(x,y,type="l",xlab="Position index", ylab="temperature",col=t_col[k])
|
120 |
|
|
axis(2,xlab="",ylab="tmax (in degree C)")
|
121 |
|
|
}
|
122 |
|
|
if (k==1 & !is.na(m)){
|
123 |
|
|
plot(x,y,type="l",col=t_col[k],axes=F) #plotting fusion profile
|
124 |
|
|
axis(4,xlab="",ylab="tmax (in degree C)")
|
125 |
|
|
|
126 |
|
|
}
|
127 |
|
|
if (k!=1 & is.na(m)){
|
128 |
|
|
#par(new=TRUE) # new plot without erasing old
|
129 |
|
|
lines(x,y,type="l",col=t_col[k],axes=F) #plotting fusion profile
|
130 |
|
|
#axis(2,xlab="",ylab="tmax (in degree C)")
|
131 |
|
|
}
|
132 |
|
|
if (k!=1 & !is.na(m)){
|
133 |
|
|
par(new=TRUE) # key: ask for new plot without erasing old
|
134 |
|
|
plot(x,y,type="l",col=t_col[k],axes=F) #plotting fusion profile
|
135 |
|
|
#axis(4,xlab="",ylab="tmax (in degree C)")
|
136 |
|
|
}
|
137 |
|
|
|
138 |
|
|
}
|
139 |
|
|
title(title_plot[i])
|
140 |
|
|
legend("topright",legend=layerNames(r_stack),
|
141 |
|
|
cex=1.2, col=t_col,
|
142 |
|
|
lty=1)
|
143 |
|
|
|
144 |
|
|
if (disp==TRUE){
|
145 |
|
|
savePlot(file=paste(list_trans[[i]][2],".png",sep=""),type="png")
|
146 |
|
|
}
|
147 |
|
|
if (disp==FALSE){
|
148 |
|
|
dev.off()
|
149 |
|
|
}
|
150 |
|
|
list_trans_data[[i]]<-trans_data
|
151 |
|
|
}
|
152 |
|
|
names(list_trans_data)<-names(list_trans)
|
153 |
|
|
return(list_trans_data)
|
154 |
|
|
}
|
155 |
|
|
|
156 |
|
|
transect_from_spdf<-function (spdf,selected_features){
|
157 |
|
|
#This function produces a transect from a set of selected points in a point layer
|
158 |
|
|
# Arguments:
|
159 |
|
|
# spdf: SpatialPointDataFrame
|
160 |
|
|
# selected_features: index of ssubset points used in the transect line
|
161 |
|
|
# Return: SpatialLinesDataframe object corresponding to the transect
|
162 |
|
|
# Author: Benoit Parmentier
|
163 |
|
|
# Date: 11-29-2012
|
164 |
|
|
|
165 |
|
|
dat_id<-spdf[selected_features,] #creating new subset from spdf
|
166 |
|
|
spdf_proj<-proj4string(dat_id)
|
167 |
|
|
matrix_point_coords<-coordinates(dat_id)
|
168 |
|
|
#Add possibility of keeping attributes?
|
169 |
|
|
#Transform a sequence of points with coords into Spatial Lines
|
170 |
|
|
#Note that X is the ID, modify for dataframe?
|
171 |
|
|
trans4<-SpatialLines(list(Lines(list(Line(coordinates(matrix_point_coords))),"X")))
|
172 |
|
|
tmp<-as.data.frame(dat_id[1,])
|
173 |
|
|
row.names(tmp)<-rep("X",1)
|
174 |
|
|
trans4<-SpatialLinesDataFrame(trans4,data=tmp)
|
175 |
|
|
proj4string(trans4)<-spdf_proj
|
176 |
|
|
return(trans4)
|
177 |
|
|
}
|
178 |
|
|
|
179 |
|
|
###Parameters and arguments
|
180 |
|
|
|
181 |
|
|
infile1<- "ghcn_or_tmax_covariates_06262012_OR83M.shp" #GHCN shapefile containing variables for modeling 2010
|
182 |
|
|
#infile2<-"list_10_dates_04212012.txt" #List of 10 dates for the regression
|
183 |
|
|
infile2<-"list_365_dates_04212012.txt" #list of dates
|
184 |
|
|
infile3<-"LST_dates_var_names.txt" #LST dates name
|
185 |
|
|
infile4<-"models_interpolation_05142012.txt" #Interpolation model names
|
186 |
|
|
infile5<-"mean_day244_rescaled.rst" #mean LST for day 244
|
187 |
|
|
inlistf<-"list_files_05032012.txt" #list of raster images containing the Covariates
|
188 |
|
|
infile6<-"OR83M_state_outline.shp"
|
189 |
|
|
#stat_loc<-read.table(paste(path,"/","location_study_area_OR_0602012.txt",sep=""),sep=",", header=TRUE)
|
190 |
|
|
|
191 |
e9998467
|
Benoit Parmentier
|
out_prefix<-"methods_comp5_12042012_"
|
192 |
30f84063
|
Benoit Parmentier
|
nb_transect<-4
|
193 |
|
|
##### LOAD USEFUL DATA
|
194 |
|
|
|
195 |
|
|
#obj_list<-"list_obj_08262012.txt" #Results of fusion from the run on ATLAS
|
196 |
e9998467
|
Benoit Parmentier
|
path_wd<-"/home/parmentier/Data/IPLANT_project/methods_interpolation_comparison_10242012" #Jupiter LOCATION on Atlas for kriging #Jupiter Location on XANDERS
|
197 |
30f84063
|
Benoit Parmentier
|
#path<-"/Users/benoitparmentier/Dropbox/Data/NCEAS/Oregon_covariates/" #Local dropbox folder on Benoit's laptop
|
198 |
e9998467
|
Benoit Parmentier
|
setwd(path_wd)
|
199 |
|
|
path_data_cai<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_CAI" #Change to constant
|
200 |
|
|
path_data_fus<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_GAM"
|
201 |
|
|
#list files that contain model objects and ratingin-testing information for CAI and Fusion
|
202 |
|
|
obj_mod_fus_name<-"results_mod_obj__365d_GAM_fusion_const_all_lstd_11022012.RData"
|
203 |
|
|
obj_mod_cai_name<-"results_mod_obj__365d_GAM_CAI2_const_all_10312012.RData"
|
204 |
|
|
|
205 |
|
|
gam_fus<-load_obj(file.path(path_data_fus,obj_mod_fus_name))
|
206 |
|
|
gam_cai<-load_obj(file.path(path_data_cai,obj_mod_cai_name)) #This contains all the info
|
207 |
|
|
sampling_date_list<-gam_fus$sampling_obj$sampling_dat$date
|
208 |
|
|
|
209 |
|
|
### Projection for the current region
|
210 |
30f84063
|
Benoit Parmentier
|
proj_str="+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs";
|
211 |
|
|
#User defined output prefix
|
212 |
|
|
|
213 |
e9998467
|
Benoit Parmentier
|
### MAKE THIS A FUNCTION TO LOAD STACK AND DEFINE VALID RANGE...
|
214 |
30f84063
|
Benoit Parmentier
|
#CRS<-proj4string(ghcn) #Storing projection information (ellipsoid, datum,etc.)
|
215 |
|
|
lines<-read.table(paste(path,"/",inlistf,sep=""), sep="") #Column 1 contains the names of raster files
|
216 |
|
|
inlistvar<-lines[,1]
|
217 |
|
|
inlistvar<-paste(path,"/",as.character(inlistvar),sep="")
|
218 |
|
|
covar_names<-as.character(lines[,2]) #Column two contains short names for covaraites
|
219 |
|
|
|
220 |
|
|
s_raster<- stack(inlistvar) #Creating a stack of raster images from the list of variables.
|
221 |
|
|
layerNames(s_raster)<-covar_names #Assigning names to the raster layers
|
222 |
|
|
projection(s_raster)<-proj_str
|
223 |
|
|
|
224 |
|
|
#Create mask using land cover data
|
225 |
|
|
pos<-match("LC10",layerNames(s_raster)) #Find the layer which contains water bodies
|
226 |
|
|
LC10<-subset(s_raster,pos)
|
227 |
|
|
LC10[is.na(LC10)]<-0 #Since NA values are 0, we assign all zero to NA
|
228 |
|
|
mask_land<-LC10<100 #All values below 100% water are assigned the value 1, value 0 is "water"
|
229 |
|
|
mask_land_NA<-mask_land
|
230 |
|
|
mask_land_NA[mask_land_NA==0]<-NA #Water bodies are assigned value 1
|
231 |
|
|
|
232 |
|
|
data_name<-"mask_land_OR"
|
233 |
|
|
raster_name<-paste(data_name,".rst", sep="")
|
234 |
|
|
writeRaster(mask_land, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI)
|
235 |
|
|
#writeRaster(r2, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI)
|
236 |
|
|
|
237 |
|
|
pos<-match("ELEV_SRTM",layerNames(s_raster)) #Find column with name "ELEV_SRTM"
|
238 |
|
|
ELEV_SRTM<-raster(s_raster,layer=pos) #Select layer from stack on 10/30
|
239 |
|
|
s_raster<-dropLayer(s_raster,pos)
|
240 |
|
|
ELEV_SRTM[ELEV_SRTM <0]<-NA
|
241 |
|
|
mask_ELEV_SRTM<-ELEV_SRTM>0
|
242 |
|
|
|
243 |
|
|
#Change this a in loop...
|
244 |
|
|
pos<-match("LC1",layerNames(s_raster)) #Find column with name "value"
|
245 |
|
|
LC1<-raster(s_raster,layer=pos) #Select layer from stack
|
246 |
|
|
s_raster<-dropLayer(s_raster,pos)
|
247 |
|
|
LC1[is.na(LC1)]<-0
|
248 |
|
|
pos<-match("LC2",layerNames(s_raster)) #Find column with name "value"
|
249 |
|
|
LC2<-raster(s_raster,layer=pos) #Select layer from stack
|
250 |
|
|
s_raster<-dropLayer(s_raster,pos)
|
251 |
|
|
LC2[is.na(LC2)]<-0
|
252 |
|
|
pos<-match("LC3",layerNames(s_raster)) #Find column with name "value"
|
253 |
|
|
LC3<-raster(s_raster,layer=pos) #Select layer from stack
|
254 |
|
|
s_raster<-dropLayer(s_raster,pos)
|
255 |
|
|
LC3[is.na(LC3)]<-0
|
256 |
|
|
pos<-match("LC4",layerNames(s_raster)) #Find column with name "value"
|
257 |
|
|
LC4<-raster(s_raster,layer=pos) #Select layer from stack
|
258 |
|
|
s_raster<-dropLayer(s_raster,pos)
|
259 |
|
|
LC4[is.na(LC4)]<-0
|
260 |
|
|
pos<-match("LC6",layerNames(s_raster)) #Find column with name "value"
|
261 |
|
|
LC6<-raster(s_raster,layer=pos) #Select layer from stack
|
262 |
|
|
s_raster<-dropLayer(s_raster,pos)
|
263 |
|
|
LC6[is.na(LC6)]<-0
|
264 |
|
|
pos<-match("LC7",layerNames(s_raster)) #Find column with name "value"
|
265 |
|
|
LC7<-raster(s_raster,layer=pos) #Select layer from stack
|
266 |
|
|
s_raster<-dropLayer(s_raster,pos)
|
267 |
|
|
LC7[is.na(LC7)]<-0
|
268 |
|
|
pos<-match("LC9",layerNames(s_raster)) #Find column with name "LC9", this is wetland...
|
269 |
|
|
LC9<-raster(s_raster,layer=pos) #Select layer from stack
|
270 |
|
|
s_raster<-dropLayer(s_raster,pos)
|
271 |
|
|
LC9[is.na(LC9)]<-0
|
272 |
|
|
|
273 |
|
|
LC_s<-stack(LC1,LC2,LC3,LC4,LC6,LC7)
|
274 |
|
|
layerNames(LC_s)<-c("LC1_forest","LC2_shrub","LC3_grass","LC4_crop","LC6_urban","LC7_barren")
|
275 |
|
|
LC_s <-mask(LC_s,mask_ELEV_SRTM)
|
276 |
|
|
plot(LC_s)
|
277 |
|
|
|
278 |
|
|
s_raster<-addLayer(s_raster, LC_s)
|
279 |
|
|
|
280 |
|
|
#mention this is the last... files
|
281 |
|
|
|
282 |
|
|
#Read region outline...
|
283 |
|
|
filename<-sub(".shp","",infile6) #Removing the extension from file.
|
284 |
|
|
reg_outline<-readOGR(".", filename) #reading shapefile
|
285 |
|
|
|
286 |
|
|
############ PART 4: RESIDUALS ANALYSIS: ranking, plots, focus regions ##################
|
287 |
|
|
############## EXAMINING STATION RESIDUALS ###########
|
288 |
|
|
########### CONSTANT OVER 365 AND SAMPLING OVER 365
|
289 |
|
|
#Plot daily_deltaclim_rast, bias_rast,add data_s and data_v
|
290 |
|
|
|
291 |
|
|
# RANK STATION by average or median RMSE
|
292 |
|
|
# Count the number of times a station is in the extremum group of outliers...
|
293 |
|
|
# LOOK at specific date...
|
294 |
|
|
|
295 |
|
|
#Examine residuals for a spciefic date...Jan, 1 using run of const_all i.e. same training over 365 dates
|
296 |
|
|
path_data_cai<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_CAI" #Change to constant
|
297 |
|
|
path_data_fus<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_GAM"
|
298 |
|
|
|
299 |
|
|
date_selected<-"20100103"
|
300 |
|
|
|
301 |
|
|
oldpath<-getwd()
|
302 |
|
|
setwd(path_data_cai)
|
303 |
|
|
|
304 |
|
|
################ VISUALIZATION !!!!!!!! ############
|
305 |
|
|
#updated the analysis
|
306 |
|
|
|
307 |
|
|
dates<-c("20100103","20100901")
|
308 |
|
|
i=2
|
309 |
|
|
|
310 |
|
|
for(i in 1:length(dates)){
|
311 |
|
|
|
312 |
|
|
date_selected<-dates[i]
|
313 |
|
|
oldpath<-getwd()
|
314 |
|
|
setwd(path_data_cai)
|
315 |
|
|
file_pat<-glob2rx(paste("*tmax_predicted*",date_selected,"*_365d_GAM_CAI2_const_all_10312012.rst",sep="")) #Search for files in relation to fusion
|
316 |
|
|
lf_cai2c<-list.files(pattern=file_pat) #Search for files in relation to fusion
|
317 |
|
|
rast_cai2c<-stack(lf_cai2c) #lf_cai2c CAI results with constant sampling over 365 dates
|
318 |
|
|
rast_cai2c<-mask(rast_cai2c,mask_ELEV_SRTM)
|
319 |
|
|
|
320 |
|
|
oldpath<-getwd()
|
321 |
|
|
setwd(path_data_fus)
|
322 |
|
|
file_pat<-glob2rx(paste("*tmax_predicted*",date_selected,"*_365d_GAM_fusion_const_all_lstd_11022012.rst",sep="")) #Search for files in relation to fusion
|
323 |
|
|
lf_fus1c<-list.files(pattern=file_pat) #Search for files in relation to fusion
|
324 |
|
|
rast_fus1c<-stack(lf_fus1c)
|
325 |
|
|
rast_fus1c<-mask(rast_fus1c,mask_ELEV_SRTM)
|
326 |
|
|
|
327 |
|
|
#PLOT ALL MODELS
|
328 |
|
|
#Prepare for plotting
|
329 |
|
|
|
330 |
|
|
setwd(path) #set path to the output path
|
331 |
|
|
|
332 |
|
|
s_range<-c(minValue(rast_fus1c),maxValue(rast_fus1c)) #stack min and max
|
333 |
|
|
s_range<-c(min(s_range),max(s_range))
|
334 |
|
|
col_breaks <- pretty(s_range, n=50)
|
335 |
|
|
lab_breaks <- pretty(s_range, n=5)
|
336 |
|
|
temp_colors <- colorRampPalette(c('blue', 'white', 'red'))
|
337 |
|
|
X11(width=18,height=12)
|
338 |
|
|
par(mfrow=c(3,3))
|
339 |
|
|
for (k in 1:length(lf_fus1c)){
|
340 |
|
|
fus1c_r<-raster(rast_fus1c,k)
|
341 |
|
|
plot(fus1c_r, breaks=col_breaks, col=temp_colors(length(col_breaks)-1),
|
342 |
|
|
axis=list(at=lab_breaks, labels=lab_breaks))
|
343 |
|
|
}
|
344 |
|
|
plot(rast_fus1c,col=temp_colors(49))
|
345 |
|
|
savePlot(paste("fig1_diff_models_fusion_",date_selected,out_prefix,".png", sep=""), type="png")
|
346 |
|
|
dev.off()
|
347 |
|
|
s_range<-c(minValue(rast_cai2c),maxValue(rast_cai2c)) #stack min and max
|
348 |
|
|
s_range<-c(min(s_range),max(s_range))
|
349 |
|
|
col_breaks <- pretty(s_range, n=50)
|
350 |
|
|
lab_breaks <- pretty(s_range, n=5)
|
351 |
|
|
temp_colors <- colorRampPalette(c('blue', 'white', 'red'))
|
352 |
|
|
X11(width=18,height=12)
|
353 |
|
|
par(mfrow=c(3,3))
|
354 |
|
|
for (k in 1:length(lf_fus1c)){
|
355 |
|
|
cai2c_r<-raster(rast_cai2c,k)
|
356 |
|
|
plot(cai2c_r, breaks=col_breaks, col=temp_colors(length(col_breaks)-1),
|
357 |
|
|
axis=list(at=lab_breaks, labels=lab_breaks))
|
358 |
|
|
}
|
359 |
|
|
plot(rast_cai2c,col=temp_colors(49))
|
360 |
|
|
savePlot(paste("fig2_diff_models_cai_",date_selected,out_prefix,".png", sep=""), type="png")
|
361 |
|
|
dev.off()
|
362 |
|
|
#PLOT CAI_Kr and Fusion_Kr
|
363 |
|
|
|
364 |
|
|
rast_fus_pred<-raster(rast_fus1c,1) # Select the first model from the stack i.e fusion with kriging for both steps
|
365 |
|
|
rast_cai_pred<-raster(rast_cai2c,1)
|
366 |
|
|
layerNames(rast_cai_pred)<-paste("cai",date_selected,sep="_")
|
367 |
|
|
layerNames(rast_fus_pred)<-paste("fus",date_selected,sep="_")
|
368 |
|
|
#Plot side by side
|
369 |
|
|
X11(width=16,height=9)
|
370 |
|
|
rast_pred<-stack(rast_cai_pred,rast_fus_pred)
|
371 |
|
|
layerNames(rast_pred)<-c(paste('CAI_kr',date_selected,sep=" "),paste('Fusion_kr',date_selected,sep=" "))
|
372 |
|
|
s.range <- c(min(minValue(rast_pred)), max(maxValue(rast_pred)))
|
373 |
|
|
col.breaks <- pretty(s.range, n=50)
|
374 |
|
|
lab.breaks <- pretty(s.range, n=5)
|
375 |
|
|
temp.colors <- colorRampPalette(c('blue', 'white', 'red'))
|
376 |
|
|
plot(rast_pred, breaks=col.breaks, col=temp.colors(length(col.breaks)-1),
|
377 |
|
|
axis=list(at=lab.breaks, labels=lab.breaks))
|
378 |
|
|
savePlot(paste("fig3_diff_CAI_fusion_",date_selected,out_prefix,".png", sep=""), type="png")
|
379 |
|
|
|
380 |
|
|
#Scatter plot of fus vs cai
|
381 |
|
|
plot(values(rast_fus_pred),values(rast_cai_pred),ylab="CAI",xlab="Fusion",axis=FALSE)
|
382 |
|
|
title(paste("CAI and fusion scatterplot on ",date_selected,sep=""))
|
383 |
|
|
savePlot(paste("fig4_diff_image_scatterplot_CAI_fusion_",date_selected,out_prefix,".png", sep=""), type="png")
|
384 |
|
|
dev.off()
|
385 |
|
|
|
386 |
|
|
## Start difference analysis
|
387 |
|
|
#Calculate difference image for the date selected
|
388 |
|
|
rast_diff<-rast_fus_pred-rast_cai_pred
|
389 |
|
|
layerNames(rast_diff)<-paste("diff",date_selected,sep="_")
|
390 |
|
|
mean_val<-cellStats(rast_diff,mean)
|
391 |
|
|
sd_val<-cellStats(rast_diff,sd)
|
392 |
|
|
|
393 |
|
|
#View classified diff and outliers...
|
394 |
|
|
diff_n_outlier<-rast_diff< (2*-sd_val) #Create negative and positive outliers...
|
395 |
|
|
diff_p_outlier<-rast_diff> (2*sd_val) #Create negative and positive outliers...
|
396 |
|
|
diff_outlier<-stack(diff_n_outlier,diff_p_outlier)
|
397 |
|
|
layerNames(diff_outlier)<-c("Negative_diff_outliers","Positive_diff_outliers")
|
398 |
|
|
bool_ramp<-colorRampPalette(c("black","red"))
|
399 |
|
|
X11()
|
400 |
|
|
plot(diff_outlier,col=bool_ramp(2))
|
401 |
|
|
savePlot(paste("fig5_diff_image_outliers_CAI_fusion_",date_selected,out_prefix,".png", sep=""), type="png")
|
402 |
|
|
dev.off()
|
403 |
|
|
tmp<-overlay(diff_n_outlier,ELEV_SRTM,fun=function(x,y){return(x*y)})
|
404 |
|
|
#could use mask
|
405 |
|
|
tmp[tmp==0]<-NA
|
406 |
|
|
mean_Elev_n_outliers<-cellStats(tmp,mean)
|
407 |
|
|
mean_Elev<-cellStats(ELEV_SRTM,mean)
|
408 |
|
|
print(c(mean_Elev_n_outliers,mean_Elev),digits=7) #This shows that outliers are in higher areas
|
409 |
|
|
# on average: 1691m compared to 1044m
|
410 |
|
|
##postive outliers and land cover
|
411 |
|
|
#LC2 (shrub), LC1(forest),LC3(grass),LC4(crop)
|
412 |
|
|
tmp<-overlay(diff_p_outlier,LC2,fun=function(x,y){return(x*y)})
|
413 |
|
|
tmp[tmp==0]<-NA
|
414 |
|
|
mean_LC2_p_outliers<-cellStats(tmp,mean) #There is more shrub (44.84% than on average 22.32)
|
415 |
|
|
mean_LC2<-cellStats(LC2,mean)
|
416 |
|
|
print(c(mean_LC2_p_outliers,mean_LC2),digits=7) #This shows that outliers have in higher
|
417 |
|
|
#proportion of shurb (44% against 25%)
|
418 |
|
|
tmp<-overlay(diff_p_outlier,LC3,fun=function(x,y){return(x*y)})
|
419 |
|
|
tmp[tmp==0]<-NA
|
420 |
|
|
mean_LC3_p_outliers<-cellStats(tmp,mean) #There is more grass (42.73% than on average 14.47)
|
421 |
|
|
mean_LC3<-cellStats(LC3,mean)
|
422 |
|
|
print(c(mean_LC3_p_outliers,mean_LC3),digits=7) #This shows that outliers have in higher
|
423 |
|
|
#proportion of shurb (44% against 25%)
|
424 |
|
|
tmp<-overlay(diff_p_outlier,LC4,fun=function(x,y){return(x*y)})
|
425 |
|
|
tmp[tmp==0]<-NA
|
426 |
|
|
mean_LC4_p_outliers<-cellStats(tmp,mean) #There is more grass (42.73% than on average 14.47)
|
427 |
|
|
mean_LC4<-cellStats(LC4,mean)
|
428 |
|
|
print(c(mean_LC4_p_outliers,mean_LC4),digits=7) #This shows that outliers have in higher
|
429 |
|
|
|
430 |
|
|
#CREATE A TABLE
|
431 |
|
|
|
432 |
|
|
####
|
433 |
|
|
#View histogram
|
434 |
|
|
hist(rast_diff)
|
435 |
|
|
|
436 |
|
|
### More Land cover analysis related to references...
|
437 |
|
|
|
438 |
|
|
LC2<-mask(LC2,mask_ELEV_SRTM)
|
439 |
|
|
cellStats(LC2,"countNA") #Check that NA have been assigned to water and areas below 0 m
|
440 |
|
|
|
441 |
|
|
LC2_50_m<- LC2>50
|
442 |
|
|
|
443 |
|
|
LC2_50<-LC2_50_m*LC2
|
444 |
|
|
diff_LC2_50<-LC2_50_m*rast_diff
|
445 |
|
|
cellStats(diff_LC2_50,"mean")
|
446 |
|
|
plot(LC2)
|
447 |
|
|
plot(LC2_50)
|
448 |
|
|
freq(LC2_50)
|
449 |
|
|
|
450 |
|
|
#Forest NOW
|
451 |
|
|
LC1<-mask(LC1,mask_ELEV_SRTM)
|
452 |
|
|
cellStats(LC1,"countNA") #Check that NA have been assigned to water and areas below 0 m
|
453 |
|
|
|
454 |
|
|
LC1_50_m<- LC1>50
|
455 |
|
|
LC1_100_m<- LC1>=100
|
456 |
|
|
LC1_50_m[LC1_50_m==0]<-NA
|
457 |
|
|
LC1_100_m[LC1_100_m==0]<-NA
|
458 |
|
|
LC1_50<-LC1_50_m*LC1
|
459 |
|
|
LC1_100<-LC1_100_m*LC1
|
460 |
|
|
plot(LC1)
|
461 |
|
|
plot(LC1_50_m)
|
462 |
|
|
freq(LC1_50_m)
|
463 |
|
|
diff_LC1_50<-LC1_50_m*rast_diff
|
464 |
|
|
diff_LC1_100<-LC1_100_m*rast_diff
|
465 |
|
|
|
466 |
|
|
plot(diff_LC1_50)
|
467 |
|
|
cellStats(diff_LC1_50,"mean")
|
468 |
|
|
cellStats(diff_LC1_100,"mean")
|
469 |
|
|
plot(values(diff_LC1_50),values(LC1_50))
|
470 |
|
|
plot(values(diff_LC1_100),values(LC1_100))
|
471 |
|
|
x<-brick(LC1,rast_diff)
|
472 |
|
|
|
473 |
|
|
#Summarize results using plot
|
474 |
|
|
#LC1 and LC3 and LC4
|
475 |
|
|
avl<-c(0,10,1,10,20,2,20,30,3,30,40,4,40,50,5,50,60,6,60,70,7,70,80,8,80,90,9,90,100,10)#Note that category 1 does not include 0!!
|
476 |
|
|
rclmat<-matrix(avl,ncol=3,byrow=TRUE)
|
477 |
|
|
LC1_rec<-reclass(LC1,rclmat) #Loss of layer names when using reclass
|
478 |
|
|
LC2_rec<-reclass(LC2,rclmat) #Loss of layer names when using reclass
|
479 |
|
|
LC3_rec<-reclass(LC3,rclmat) #Loss of layer names when using reclass
|
480 |
|
|
LC4_rec<-reclass(LC4,rclmat) #Loss of layer names when using reclass
|
481 |
|
|
LC6_rec<-reclass(LC6,rclmat) #Loss of layer names when using reclass
|
482 |
|
|
|
483 |
|
|
#LC_s<-stack(LC1,LC3,LC4,LC6)
|
484 |
|
|
LC_s<-stack(LC1,LC2,LC3,LC4,LC6)
|
485 |
|
|
layerNames(LC_s)<-c("LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban")
|
486 |
|
|
LC_s<-mask(LC_s,mask_ELEV_SRTM)
|
487 |
|
|
LC_rec_s<-reclass(LC_s,rclmat)
|
488 |
|
|
|
489 |
|
|
#plot average difference per class of forest and LC2
|
490 |
|
|
rast_stack_zones<-LC_rec_s
|
491 |
|
|
|
492 |
|
|
avg_LC1_rec<-zonal(rast_diff,zones=LC1_rec,stat="mean",na.rm=TRUE)
|
493 |
|
|
avg_LC2_rec<-zonal(rast_diff,zones=LC2_rec,stat="mean",na.rm=TRUE)
|
494 |
|
|
avg_LC3_rec<-zonal(rast_diff,zones=LC3_rec,stat="mean",na.rm=TRUE)
|
495 |
|
|
avg_LC4_rec<-zonal(rast_diff,zones=LC4_rec,stat="mean",na.rm=TRUE)
|
496 |
|
|
avg_LC6_rec<-zonal(rast_diff,zones=LC6_rec,stat="mean",na.rm=TRUE)
|
497 |
|
|
|
498 |
|
|
std_LC1_rec<-zonal(rast_diff,zones=LC1_rec,stat="sd",na.rm=TRUE)
|
499 |
|
|
std_LC2_rec<-zonal(rast_diff,zones=LC2_rec,stat="sd",na.rm=TRUE)
|
500 |
|
|
std_LC3_rec<-zonal(rast_diff,zones=LC3_rec,stat="sd",na.rm=TRUE)
|
501 |
|
|
std_LC4_rec<-zonal(rast_diff,zones=LC4_rec,stat="sd",na.rm=TRUE)
|
502 |
|
|
std_LC6_rec<-zonal(rast_diff,zones=LC6_rec,stat="sd",na.rm=TRUE)
|
503 |
|
|
|
504 |
|
|
avg_LC_rec<-zonal(rast_diff,zones=LC_rec_s,stat="mean",na.rm=TRUE)
|
505 |
|
|
std_LC_rec<-zonal(rast_diff,zones=LC_rec_s,stat="sd",na.rm=TRUE)
|
506 |
|
|
|
507 |
|
|
zones_stat_std<-as.data.frame(cbind(std_LC1_rec,std_LC2_rec[,2],std_LC3_rec[,2],std_LC4_rec[,2],std_LC6_rec[,2]))
|
508 |
|
|
zones_stat<-as.data.frame(cbind(avg_LC1_rec,avg_LC2_rec[,2],avg_LC3_rec[,2],avg_LC4_rec[,2],avg_LC6_rec[,2]))
|
509 |
|
|
names(zones_stat)<-c("zones","LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban")
|
510 |
|
|
names(zones_stat_std)<-c("zones","LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban")
|
511 |
|
|
|
512 |
|
|
X11()
|
513 |
|
|
plot(zones_stat$zones,zones_stat$LC1_forest,type="b",ylim=c(-4.5,4.5),
|
514 |
|
|
ylab="difference between CAI and fusion",xlab="land cover percent class/10")
|
515 |
|
|
lines(zones_stat$zones,zones_stat[,3],col="red",type="b")
|
516 |
|
|
lines(zones_stat$zones,zones_stat[,4],col="blue",type="b")
|
517 |
|
|
lines(zones_stat$zones,zones_stat[,5],col="darkgreen",type="b")
|
518 |
|
|
lines(zones_stat$zones,zones_stat[,6],col="purple",type="b")
|
519 |
|
|
legend("topleft",legend=c("LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban"),
|
520 |
|
|
cex=1.2, col=c("black","red","blue","darkgreen","purple"),
|
521 |
|
|
lty=1)
|
522 |
|
|
title(paste("Prediction tmax difference and land cover ",sep=""))
|
523 |
|
|
|
524 |
|
|
savePlot(paste("fig6_diff_prediction_tmax_difference_land cover",date_selected,out_prefix,".png", sep=""), type="png")
|
525 |
|
|
dev.off()
|
526 |
|
|
|
527 |
|
|
avl<-c(0,500,1,500,1000,2,1000,1500,3,1500,2000,4,2000,2500,5,2500,4000,6)
|
528 |
|
|
rclmat<-matrix(avl,ncol=3,byrow=TRUE)
|
529 |
|
|
elev_rec<-reclass(ELEV_SRTM,rclmat) #Loss of layer names when using reclass
|
530 |
|
|
|
531 |
|
|
elev_rec_forest<-elev_rec*LC1_100_m
|
532 |
|
|
avg_elev_rec<-zonal(rast_diff,zones=elev_rec,stat="mean",na.rm=TRUE)
|
533 |
|
|
std_elev_rec<-zonal(rast_diff,zones=elev_rec,stat="sd",na.rm=TRUE)
|
534 |
|
|
avg_elev_rec_forest<-zonal(rast_diff,zones=elev_rec_forest,stat="mean",na.rm=TRUE)
|
535 |
|
|
std_elev_rec_forest<-zonal(rast_diff,zones=elev_rec_forest,stat="sd",na.rm=TRUE)
|
536 |
|
|
|
537 |
|
|
## CREATE plots
|
538 |
|
|
X11()
|
539 |
|
|
plot(avg_elev_rec[,1],avg_elev_rec[,2],type="b",ylim=c(-10,1),
|
540 |
|
|
ylab="difference between CAI and fusion",xlab="elevation classes")
|
541 |
|
|
lines(avg_elev_rec_forest[,1],avg_elev_rec_forest[,2],col="green",type="b") #Elevation and 100% forest...
|
542 |
|
|
legend("topright",legend=c("Elevation", "elev_forest"),
|
543 |
|
|
cex=1.2, col=c("black","darkgreen"),
|
544 |
|
|
lty=1)
|
545 |
|
|
title(paste("Prediction tmax difference and elevation ",sep=""))
|
546 |
|
|
savePlot(paste("fig7_diff_prediction_tmax_difference_elevation",date_selected,out_prefix,".png", sep=""), type="png")
|
547 |
|
|
dev.off()
|
548 |
|
|
#Add plots with std as CI
|
549 |
|
|
|
550 |
|
|
}
|
551 |
|
|
|
552 |
|
|
###################################################################
|
553 |
e9998467
|
Benoit Parmentier
|
################ SPATIAL TRANSECT THROUGH THE IMAGE: ####################
|
554 |
30f84063
|
Benoit Parmentier
|
|
555 |
|
|
#select date
|
556 |
|
|
dates<-c("20100103","20100901")
|
557 |
e9998467
|
Benoit Parmentier
|
#j=2
|
558 |
30f84063
|
Benoit Parmentier
|
|
559 |
|
|
for (j in 1:length(dates)){
|
560 |
|
|
|
561 |
|
|
#Read predicted tmax raster surface and modeling information
|
562 |
|
|
date_selected<-dates[j]
|
563 |
e9998467
|
Benoit Parmentier
|
|
564 |
30f84063
|
Benoit Parmentier
|
oldpath<-getwd()
|
565 |
|
|
setwd(path_data_cai)
|
566 |
|
|
file_pat<-glob2rx(paste("*tmax_predicted*",date_selected,"*_365d_GAM_CAI2_const_all_10312012.rst",sep="")) #Search for files in relation to fusion
|
567 |
|
|
lf_cai2c<-list.files(pattern=file_pat) #Search for files in relation to fusion
|
568 |
|
|
rast_cai2c<-stack(lf_cai2c) #lf_cai2c CAI results with constant sampling over 365 dates
|
569 |
|
|
rast_cai2c<-mask(rast_cai2c,mask_ELEV_SRTM)
|
570 |
|
|
|
571 |
|
|
oldpath<-getwd()
|
572 |
|
|
setwd(path_data_fus)
|
573 |
|
|
file_pat<-glob2rx(paste("*tmax_predicted*",date_selected,"*_365d_GAM_fusion_const_all_lstd_11022012.rst",sep="")) #Search for files in relation to fusion
|
574 |
|
|
lf_fus1c<-list.files(pattern=file_pat) #Search for files in relation to fusion
|
575 |
|
|
rast_fus1c<-stack(lf_fus1c)
|
576 |
|
|
rast_fus1c<-mask(rast_fus1c,mask_ELEV_SRTM)
|
577 |
|
|
|
578 |
|
|
setwd(path)
|
579 |
|
|
rast_fus_pred<-raster(rast_fus1c,1)
|
580 |
|
|
rast_cai_pred<-raster(rast_cai2c,1)
|
581 |
|
|
rast_diff_fc<-rast_fus_pred-rast_cai_pred
|
582 |
e9998467
|
Benoit Parmentier
|
|
583 |
30f84063
|
Benoit Parmentier
|
#Read in data_s and data_v
|
584 |
e9998467
|
Benoit Parmentier
|
|
585 |
|
|
k<-match(date_selected,sampling_date_list)
|
586 |
|
|
names(gam_fus$gam_fus_mod[[k]]) #Show the name structure of the object/list
|
587 |
30f84063
|
Benoit Parmentier
|
|
588 |
e9998467
|
Benoit Parmentier
|
#Extract the training and testing information for the given date...
|
589 |
|
|
data_sf<-gam_fus$gam_fus_mod[[k]]$data_s #object for the first date...20100103 #Make this a function??
|
590 |
|
|
data_vf<-gam_fus$gam_fus_mod[[k]]$data_v #object for the first date...20100103
|
591 |
|
|
data_sc<-gam_cai$gam_CAI_mod[[k]]$data_s #object for the first date...20100103
|
592 |
|
|
data_vc<-gam_cai$gam_CAI_mod[[k]]$data_v #object for the first date...20100103
|
593 |
30f84063
|
Benoit Parmentier
|
|
594 |
|
|
### CREATE A NEW TRANSECT BASED ON LOCATION OF SPECIFIED STATIONS
|
595 |
|
|
|
596 |
|
|
selected_stations<-c("USW00024284","USC00354126","USC00358536","USC00354835",
|
597 |
|
|
"USC00356252","USC00359316","USC00358246","USC00350694",
|
598 |
|
|
"USC00350699","USW00024230","USC00353542")
|
599 |
|
|
#add which one were training and testing
|
600 |
|
|
data_vf$training<-rep(0,nrow(data_vf))
|
601 |
|
|
data_sf$training<-rep(1,nrow(data_sf))
|
602 |
|
|
|
603 |
e9998467
|
Benoit Parmentier
|
data_stat<-rbind(data_vf[,c("id","training")],data_sf[,c("id","training")]) #bringing together data_v and data_s
|
604 |
30f84063
|
Benoit Parmentier
|
m<-match(selected_stations,data_stat$id)
|
605 |
e9998467
|
Benoit Parmentier
|
m<-as.integer(na.omit(m))
|
606 |
30f84063
|
Benoit Parmentier
|
trans4_stations<-transect_from_spdf(data_stat,m)
|
607 |
e9998467
|
Benoit Parmentier
|
point4_stations<-data_stat[m,]
|
608 |
30f84063
|
Benoit Parmentier
|
#tmp<-as.data.frame(data_stat[1,])
|
609 |
|
|
#row.names(tmp)<-rep("X",1)
|
610 |
|
|
#test<-SpatialLinesDataFrame(trans4_stations,data=tmp)
|
611 |
|
|
writeOGR(obj=trans4_stations,layer="t4_line",dsn="t4_line.shp",driver="ESRI Shapefile", overwrite=T)
|
612 |
|
|
## Create list of transect
|
613 |
|
|
|
614 |
|
|
list_transect<-vector("list",nb_transect)
|
615 |
|
|
list_transect[[1]]<-c("t1_line.shp",paste("figure_9_tmax_transect1_OR",date_selected,out_prefix,sep="_"))
|
616 |
|
|
list_transect[[2]]<-c("t2_line.shp",paste("figure_10_tmax_transect2_OR",date_selected,out_prefix,sep="_"))
|
617 |
|
|
list_transect[[3]]<-c("t3_line.shp",paste("figure_11_tmax_transect3_OR",date_selected,out_prefix,sep="_"))
|
618 |
|
|
list_transect[[4]]<-c("t4_line.shp",paste("figure_12_tmax_transect4_OR",date_selected,out_prefix,sep="_"))
|
619 |
|
|
|
620 |
|
|
names(list_transect)<-c("transect_OR1","transect_OR2","transect_OR3","transect_OR4")
|
621 |
|
|
|
622 |
|
|
#now add a transect for elevation
|
623 |
|
|
list_transect2<-vector("list",nb_transect)
|
624 |
|
|
list_transect2[[1]]<-c("t1_line.shp",paste("figure_13_tmax_elevation_transect1_OR",date_selected,out_prefix,sep="_"))
|
625 |
|
|
list_transect2[[2]]<-c("t2_line.shp",paste("figure_14_tmax_elevation_transect2_OR",date_selected,out_prefix,sep="_"))
|
626 |
|
|
list_transect2[[3]]<-c("t3_line.shp",paste("figure_15_tmax_elevation_transect3_OR",date_selected,out_prefix,sep="_"))
|
627 |
e9998467
|
Benoit Parmentier
|
list_transect2[[4]]<-c("t4_line.shp",paste("figure_16_tmax_elevation_transect4_OR",date_selected,out_prefix,sep="_"))
|
628 |
30f84063
|
Benoit Parmentier
|
|
629 |
|
|
names(list_transect2)<-c("transect_OR1","transect_OR2","transect_OR3","transect_OR4")
|
630 |
|
|
|
631 |
|
|
rast_pred<-stack(rast_fus_pred,rast_cai_pred)
|
632 |
|
|
rast_pred2<-stack(rast_fus_pred,rast_cai_pred,ELEV_SRTM)
|
633 |
|
|
layerNames(rast_pred)<-c("fus","CAI")
|
634 |
|
|
layerNames(rast_pred2)<-c("fus","CAI","elev")
|
635 |
|
|
title_plot<-paste(names(list_transect),date_selected)
|
636 |
|
|
title_plot2<-paste(names(list_transect2),date_selected)
|
637 |
|
|
#r_stack<-rast_pred
|
638 |
|
|
|
639 |
|
|
X11(width=9,height=9)
|
640 |
|
|
nb_transect<-length(list_transect)
|
641 |
|
|
s_range<-c(minValue(rast_diff_fc),maxValue(rast_diff_fc)) #stack min and max
|
642 |
|
|
col_breaks <- pretty(s_range, n=50)
|
643 |
|
|
lab_breaks <- pretty(s_range, n=7)
|
644 |
|
|
temp_colors <- colorRampPalette(c('blue', 'white', 'red'))
|
645 |
|
|
plot(rast_diff_fc, breaks=col_breaks, col=temp_colors(length(col_breaks)-1),
|
646 |
|
|
axis=list(at=lab_breaks, labels=lab_breaks))
|
647 |
|
|
for (k in 1:nb_transect){
|
648 |
|
|
trans_file<-list_transect[[k]]
|
649 |
|
|
filename<-sub(".shp","",trans_file) #Removing the extension from file.
|
650 |
|
|
transect<-readOGR(".", filename) #reading shapefile
|
651 |
|
|
plot(transect,add=TRUE)
|
652 |
|
|
}
|
653 |
|
|
|
654 |
|
|
savePlot(paste("fig8_diff_transect_path_tmax_diff_CAI_fusion_",date_selected,out_prefix,".png", sep=""), type="png")
|
655 |
|
|
dev.off()
|
656 |
|
|
|
657 |
|
|
X11(width=18,height=9)
|
658 |
|
|
m_layers_sc<-c(3)
|
659 |
|
|
trans_data<-plot_transect(list_transect,rast_pred,title_plot,disp=TRUE)
|
660 |
|
|
|
661 |
|
|
trans_data2<-plot_transect_m(list_transect2,rast_pred2,title_plot2,disp=TRUE,m_layers_sc)
|
662 |
|
|
dev.off()
|
663 |
|
|
|
664 |
e9998467
|
Benoit Parmentier
|
### PLOT LOCATIONS OF STATION ON FIGURES
|
665 |
30f84063
|
Benoit Parmentier
|
|
666 |
e9998467
|
Benoit Parmentier
|
data_stat<-rbind(data_vf[,c("id","training")],data_sf[,c("id","training")]) #bringing together data_v and data_s
|
667 |
|
|
m<-match(selected_stations,data_stat$id)
|
668 |
|
|
m<-as.integer(na.omit(m))
|
669 |
|
|
trans4_stations<-transect_from_spdf(data_stat,m)
|
670 |
|
|
point4_stations<-data_stat[m,]
|
671 |
|
|
|
672 |
|
|
pos<-match(c("x_OR83M","y_OR83M"),layerNames(s_raster)) #Find column with name "value"
|
673 |
|
|
xy_stack<-subset(s_raster,pos) #Select multiple layers from the stack
|
674 |
|
|
r_stack<-stack(xy_stack, rast_pred2)
|
675 |
|
|
trans4_data<-extract(r_stack,trans4_stations,cellnumbers=TRUE) #This extracts a list
|
676 |
|
|
trans4_data<-as.data.frame(trans4_data[[1]])
|
677 |
|
|
point4_cellID<-cellFromXY(r_stack,coordinates(point4_stations)) #This contains the cell ID the points
|
678 |
|
|
pos<-match(point4_cellID,trans4_data$cell)
|
679 |
|
|
|
680 |
|
|
#Plots lines where there are stations...
|
681 |
|
|
X11(width=18,height=9)
|
682 |
|
|
y<-trans4_data$fus
|
683 |
|
|
x <- 1:length(y)
|
684 |
|
|
plot(x,y,type="l",col="red", #plotting fusion profile
|
685 |
|
|
,xlab="",ylab="tmax (in degree C)")
|
686 |
|
|
y<-trans4_data$CAI
|
687 |
|
|
lines(x,y,col="green")
|
688 |
|
|
abline(v=pos)#addlines whtere the stations area...
|
689 |
|
|
#plot(elev)
|
690 |
|
|
#title(title_plot[i]))
|
691 |
|
|
legend("topleft",legend=c("fus","CAI"),
|
692 |
|
|
cex=1.2, col=c("red","green"),
|
693 |
|
|
lty=1)
|
694 |
|
|
savePlot(paste("fig17_transect_path_tmax_diff_CAI_fusion_",date_selected,out_prefix,".png", sep=""), type="png")
|
695 |
|
|
|
696 |
|
|
|
697 |
|
|
y<-trans4_data$fus[1:150]
|
698 |
|
|
x <- 1:150
|
699 |
|
|
plot(x,y,type="l",col="red", #plotting fusion profile
|
700 |
|
|
,xlab="",ylab="tmax (in degree C)")
|
701 |
|
|
y<-trans4_data$CAI[1:150]
|
702 |
|
|
lines(x,y,col="green")
|
703 |
|
|
abline(v=pos)#addlines whtere the stations area...
|
704 |
|
|
#plot(elev)
|
705 |
|
|
#title(title_plot[i]))
|
706 |
|
|
legend("topleft",legend=c("fus","CAI"),
|
707 |
|
|
cex=1.2, col=c("red","green"),
|
708 |
30f84063
|
Benoit Parmentier
|
lty=1)
|
709 |
e9998467
|
Benoit Parmentier
|
savePlot(paste("fig18a_transect_path_tmax_diff_CAI_fusion_",date_selected,out_prefix,".png", sep=""), type="png")
|
710 |
|
|
|
711 |
|
|
|
712 |
|
|
y<-trans4_data$fus[151:300]
|
713 |
|
|
x <- 151:300
|
714 |
|
|
plot(x,y,type="l",col="red", #plotting fusion profile
|
715 |
|
|
,xlab="",ylab="tmax (in degree C)")
|
716 |
|
|
y<-trans4_data$CAI[151:300]
|
717 |
|
|
lines(x,y,col="green")
|
718 |
|
|
abline(v=pos)#addlines whtere the stations area...
|
719 |
|
|
#plot(elev)
|
720 |
|
|
#title(title_plot[i]))
|
721 |
|
|
legend("topleft",legend=c("fus","CAI"),
|
722 |
|
|
cex=1.2, col=c("red","green"),
|
723 |
|
|
lty=1)
|
724 |
|
|
savePlot(paste("fig18b_transect_path_tmax_diff_CAI_fusion_",date_selected,out_prefix,".png", sep=""), type="png")
|
725 |
|
|
|
726 |
|
|
dev.off()
|
727 |
|
|
|
728 |
|
|
#cor(fus_y,elev_y)
|
729 |
|
|
#cor(cai_y,elev_y)
|
730 |
|
|
#cor(fus_y,cai_y)
|
731 |
30f84063
|
Benoit Parmentier
|
|
732 |
|
|
}
|