1 |
2a0267ad
|
Benoit Parmentier
|
##################################### METHODS COMPARISON part 3 ##########################################
|
2 |
|
|
#################################### Spatial Analysis ############################################
|
3 |
|
|
#This script utilizes the R ojbects created during the interpolation phase. #
|
4 |
|
|
#At this stage the script produces figures of various accuracy metrics and compare methods: #
|
5 |
|
|
#- rank station in terms of residuals and metrics (MAE and RMSE) #
|
6 |
|
|
#- calculate accuary metrics for climtology predictions... #
|
7 |
|
|
#- spatial density of station network and accuracy metric #
|
8 |
|
|
#- visualization of maps of prediction and difference for comparison #
|
9 |
|
|
#AUTHOR: Benoit Parmentier #
|
10 |
|
|
#DATE: 10/23/2012 #
|
11 |
|
|
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#491 -- #
|
12 |
|
|
###################################################################################################
|
13 |
|
|
|
14 |
|
|
###Loading R library and packages
|
15 |
|
|
library(gtools) # loading some useful tools such as mixedsort
|
16 |
|
|
library(mgcv) # GAM package by Wood 2006 (version 2012)
|
17 |
|
|
library(sp) # Spatial pacakge with class definition by Bivand et al. 2008
|
18 |
|
|
library(spdep) # Spatial package with methods and spatial stat. by Bivand et al. 2012
|
19 |
|
|
library(rgdal) # GDAL wrapper for R, spatial utilities (Keitt et al. 2012)
|
20 |
|
|
library(gstat) # Kriging and co-kriging by Pebesma et al. 2004
|
21 |
|
|
library(automap) # Automated Kriging based on gstat module by Hiemstra et al. 2008
|
22 |
|
|
library(spgwr)
|
23 |
|
|
library(gpclib)
|
24 |
|
|
library(maptools)
|
25 |
|
|
library(graphics)
|
26 |
|
|
library(parallel) # Urbanek S. and Ripley B., package for multi cores & parralel processing
|
27 |
|
|
library(raster)
|
28 |
|
|
library(rasterVis)
|
29 |
|
|
library(plotrix) #Draw circle on graph
|
30 |
|
|
library(reshape)
|
31 |
|
|
## Functions
|
32 |
|
|
#loading R objects that might have similar names
|
33 |
|
|
load_obj <- function(f)
|
34 |
|
|
{
|
35 |
|
|
env <- new.env()
|
36 |
|
|
nm <- load(f, env)[1]
|
37 |
|
|
env[[nm]]
|
38 |
|
|
}
|
39 |
|
|
|
40 |
|
|
###Parameters and arguments
|
41 |
|
|
|
42 |
|
|
infile1<- "ghcn_or_tmax_covariates_06262012_OR83M.shp" #GHCN shapefile containing variables for modeling 2010
|
43 |
|
|
#infile2<-"list_10_dates_04212012.txt" #List of 10 dates for the regression
|
44 |
|
|
infile2<-"list_365_dates_04212012.txt" #list of dates
|
45 |
|
|
infile3<-"LST_dates_var_names.txt" #LST dates name
|
46 |
|
|
infile4<-"models_interpolation_05142012.txt" #Interpolation model names
|
47 |
|
|
infile5<-"mean_day244_rescaled.rst" #mean LST for day 244
|
48 |
|
|
inlistf<-"list_files_05032012.txt" #list of raster images containing the Covariates
|
49 |
|
|
infile6<-"OR83M_state_outline.shp"
|
50 |
|
|
#stat_loc<-read.table(paste(path,"/","location_study_area_OR_0602012.txt",sep=""),sep=",", header=TRUE)
|
51 |
|
|
|
52 |
c352e9e1
|
Benoit Parmentier
|
out_prefix<-"methods_11072012_"
|
53 |
2a0267ad
|
Benoit Parmentier
|
#output path
|
54 |
|
|
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_07192012_GAM"
|
55 |
c352e9e1
|
Benoit Parmentier
|
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_07192012_GAM"
|
56 |
|
|
|
57 |
2a0267ad
|
Benoit Parmentier
|
setwd(path)
|
58 |
|
|
|
59 |
|
|
##### LOAD USEFUL DATA
|
60 |
|
|
|
61 |
|
|
obj_list<-"list_obj_08262012.txt" #Results of fusion from the run on ATLAS
|
62 |
c352e9e1
|
Benoit Parmentier
|
path<-"/home/parmentier/Data/IPLANT_project/methods_interpolation_comparison_10242012" #Jupiter LOCATION on Atlas for kriging #Jupiter Location on XANDERS
|
63 |
2a0267ad
|
Benoit Parmentier
|
#path<-"/Users/benoitparmentier/Dropbox/Data/NCEAS/Oregon_covariates/" #Local dropbox folder on Benoit's laptop
|
64 |
|
|
setwd(path)
|
65 |
|
|
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";
|
66 |
|
|
#User defined output prefix
|
67 |
|
|
|
68 |
|
|
sampling_CAI<-load_obj("results2_CAI_sampling_obj_09132012_365d_GAM_CAI2_multisampling2.RData")
|
69 |
|
|
sampling_fus<-load_obj("results2_fusion_sampling_obj_10d_GAM_fusion_multisamp4_09192012.RData")
|
70 |
|
|
fus_CAI_mod<-load_obj("results2_CAI_Assessment_measure_all_09132012_365d_GAM_CAI2_multisampling2.RData")
|
71 |
|
|
gam_fus_mod1<-load_obj("results2_fusion_Assessment_measure_all_10d_GAM_fusion_multisamp4_09192012.RData")
|
72 |
|
|
|
73 |
|
|
filename<-sub(".shp","",infile1) #Removing the extension from file.
|
74 |
|
|
ghcn<-readOGR(".", filename) #reading shapefile
|
75 |
c352e9e1
|
Benoit Parmentier
|
filename<-sub(".shp","",infile6) #Removing the extension from file.
|
76 |
|
|
#reg_outline<-readOGR(".", filename) #reading shapefile
|
77 |
|
|
reg_outline<-readOGR(dsn=".", layer=filename) #reading shapefile
|
78 |
2a0267ad
|
Benoit Parmentier
|
|
79 |
|
|
#CRS<-proj4string(ghcn) #Storing projection information (ellipsoid, datum,etc.)
|
80 |
|
|
lines<-read.table(paste(path,"/",inlistf,sep=""), sep="") #Column 1 contains the names of raster files
|
81 |
|
|
inlistvar<-lines[,1]
|
82 |
|
|
inlistvar<-paste(path,"/",as.character(inlistvar),sep="")
|
83 |
|
|
covar_names<-as.character(lines[,2]) #Column two contains short names for covaraites
|
84 |
|
|
|
85 |
|
|
s_raster<- stack(inlistvar) #Creating a stack of raster images from the list of variables.
|
86 |
|
|
layerNames(s_raster)<-covar_names #Assigning names to the raster layers
|
87 |
|
|
projection(s_raster)<-proj_str
|
88 |
|
|
|
89 |
|
|
#Create mask using land cover data
|
90 |
|
|
pos<-match("LC10",layerNames(s_raster)) #Find the layer which contains water bodies
|
91 |
|
|
LC10<-subset(s_raster,pos)
|
92 |
|
|
LC10[is.na(LC10)]<-0 #Since NA values are 0, we assign all zero to NA
|
93 |
|
|
#LC10_mask[LC10_mask<100]<-1
|
94 |
|
|
mask_land<-LC10==100 #All values below 100% water are assigned the value 1, value 0 is "water"
|
95 |
|
|
mask_land_NA<-mask_land
|
96 |
|
|
mask_land_NA[mask_land_NA==0]<-NA #Water bodies are assigned value 1
|
97 |
|
|
|
98 |
|
|
data_name<-"mask_land_OR"
|
99 |
|
|
raster_name<-paste(data_name,".rst", sep="")
|
100 |
|
|
writeRaster(mask_land, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI)
|
101 |
|
|
#writeRaster(r2, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI)
|
102 |
|
|
|
103 |
|
|
pos<-match("mm_01",layerNames(s_raster))
|
104 |
|
|
mm_01<-subset(s_raster,pos)
|
105 |
|
|
mm_01<-mm_01-273.15
|
106 |
|
|
mm_01<-mask(mm_01,mask_land_NA) #Raster image used as backround
|
107 |
|
|
#mention this is the last... files
|
108 |
|
|
|
109 |
c352e9e1
|
Benoit Parmentier
|
################# PART 1: COMPARING RESULTS USING ALL MONTHNLY DATA AND SAMPLED FROM DAILY ###################
|
110 |
|
|
######## Visualize data: USING ALL MONTHLY DATA...
|
111 |
2a0267ad
|
Benoit Parmentier
|
|
112 |
|
|
file_pat<-glob2rx("GAMCAI_tmax_pred_predicted_mod*_365d_GAM_CAI2_all_lstd_10262012.rst")
|
113 |
|
|
lf_cai_all<-list.files(pattern=file_pat)
|
114 |
|
|
image_file<-grep(file_pat,lf_ck,value=TRUE) # using grep with "value" extracts the matching names
|
115 |
|
|
raster_pred_k<-raster(image_file)
|
116 |
|
|
|
117 |
|
|
#### CHECKING RESULTS ON 10/28 BY VISUALIZING PREDICTIONS... DO THAT FOR ALL MODELS AND TWO DATES...
|
118 |
|
|
|
119 |
|
|
date_selected<-"20100103"
|
120 |
|
|
lf_fus<-list.files(pattern=paste("GAM_bias_tmax_predicted_mod8_",date_selected,"_30_1_365d_GAM_fusion_all_lstd_10242012.rst$",sep="")) #Search for files in relation to fusion
|
121 |
|
|
lf_cai<-list.files(pattern=paste("*CAI_tmax_pred.*",date_selected,"*.08072012_365d_GAM_CAI2.rst$",sep="")) #Search for files in relation to fusion
|
122 |
|
|
|
123 |
|
|
lf_cai<-list.files(pattern=paste("GAMCAI_tmax_pred_predicted_mod8_20100103_30_1_365d_GAM_CAI2_all_lstd_10262012.rst",sep=""))
|
124 |
|
|
|
125 |
|
|
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_GAM"
|
126 |
|
|
setwd(path)
|
127 |
|
|
tmax_predicted_mod8<-raster("GAM_bias_tmax_predicted_mod8_20100103_30_1_365d_GAM_fusion_all_lstd_10242012.rst")
|
128 |
|
|
tmax_bias_mod8<-raster("GAM_bias_predicted_mod8_20100103_30_1_365d_GAM_fusion_all_lstd_10242012.rst")
|
129 |
|
|
daily_delta_tmax<-raster("fusion_daily_delta_20100103_30_1_365d_GAM_fusion_all_lstd_10242012.rst")
|
130 |
|
|
fusion<-raster("fusion_tmax_predicted_20100103_30_1_365d_GAM_fusion_all_lstd_10242012.rst")
|
131 |
|
|
LST<-(raster("mean_month1_rescaled.rst"))-273.16
|
132 |
|
|
|
133 |
|
|
tmp_rast<-LST+daily_delta_tmax-tmax_bias_mod8
|
134 |
|
|
dif_rast<-tmp_rast-tmax_predicted_mod8
|
135 |
|
|
eq<-tmp_rast==tmax_predicted_mod8
|
136 |
|
|
plot(tmax_predicted_mod8)
|
137 |
|
|
|
138 |
|
|
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_CAI"
|
139 |
|
|
setwd(path)
|
140 |
|
|
cai_tmax_mod8_rast<-raster("GAMCAI_tmax_pred_predicted_mod8_20100103_30_1_365d_GAM_CAI2_all_lstd_10262012.rst")
|
141 |
|
|
cai_delta_rast<-raster("CAI_deltaclim_20100103_30_1_365d_GAM_CAI2_all_lstd_10262012.rst")
|
142 |
|
|
cai_clim_rast<-raster("CAI_clim_20100103_30_1_365d_GAM_CAI2_all_lstd_10262012.rst")
|
143 |
|
|
cai_climod8_rast<-raster("GAMCAI_predicted_mod8_20100103_30_1_365d_GAM_CAI2_all_lstd_10262012.rst")
|
144 |
|
|
|
145 |
|
|
# s.range <- c(min(minValue(diff)), max(maxValue(diff)))
|
146 |
|
|
# col.breaks <- pretty(s.range, n=50)
|
147 |
|
|
# lab.breaks <- pretty(s.range, n=5)
|
148 |
|
|
temp.colors <- colorRampPalette(c('blue', 'white', 'red'))
|
149 |
|
|
# plot(diff, breaks=col.breaks, col=temp.colors(length(col.breaks)-1)
|
150 |
|
|
#
|
151 |
|
|
pred_cai_mod8<-stack(cai_climod8_rast,cai_delta_rast,cai_tmax_mod8_rast)
|
152 |
|
|
layerNames(pred_cai_mod8)<-paste(c("cai_climod8_rast","cai_delta_rast","cai_tmax_mod8_rast"),"20100103",sep=" ")
|
153 |
|
|
pred_fus_mod8<-stack(tmax_bias_mod8,daily_delta_tmax,LST,tmax_predicted_mod8)
|
154 |
|
|
layerNames(pred_fus_mod8)<-paste(c("tmax_bias_mod8","daily_delta_tmax","LST","tmax_predicted_mod8"),"20100103",sep=" ")
|
155 |
|
|
pred_fus_cai_comp<-stack(pred_fus_mod8,pred_cai_mod8)
|
156 |
|
|
|
157 |
|
|
X11(12,12)
|
158 |
|
|
#for (i in 1:)
|
159 |
|
|
#plot(pred_cai_mod8,col=temp.colors)
|
160 |
|
|
plot(pred_cai_mod8,col=temp.colors(25))
|
161 |
|
|
#title("prediction for 20100103")
|
162 |
|
|
plot(pred_fus_mod8,col=temp.colors(25))
|
163 |
|
|
plot(pred_fus_cai_comp,col=temp.colors(25))
|
164 |
|
|
dev.off
|
165 |
|
|
tmp_rast2<-cai_delta_rast + cai_climod8_rast
|
166 |
|
|
dif2_rast<-tmp_rast2-cai_tmax_mod8_rast
|
167 |
|
|
|
168 |
|
|
dif2<-cai_rast-tmax_predicted_mod8 #This seems to be the same!!!
|
169 |
|
|
plot(cai_climod8_rast) #THe bias surface is different when compared to the climatology based on mod8
|
170 |
|
|
plot(tmax_bias_mod8)
|
171 |
|
|
plot(tmax_predicted_mod8)
|
172 |
|
|
plot(cai_tmax_rast)
|
173 |
|
|
|
174 |
|
|
cellStats(pred_fus_cai_comp)
|
175 |
|
|
|
176 |
|
|
#t<-"results2_CAI_Assessment_measure_all_08072012_365d_GAM_CAI2.RData"
|
177 |
|
|
|
178 |
|
|
## Examining data with simplified CAI3 models in comparison to other: 11/03/2012
|
179 |
|
|
|
180 |
|
|
oldpath<-getwd()
|
181 |
|
|
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_CAI"
|
182 |
|
|
setwd(path)
|
183 |
|
|
|
184 |
|
|
date_selected<-"20100103"
|
185 |
|
|
#lf_fus<-list.files(pattern=paste("GAM_bias_tmax_predicted_mod8_",date_selected,"_30_1_365d_GAM_fusion_all_lstd_10242012.rst$",sep="")) #Search for files in relation to fusion
|
186 |
|
|
file_pat<-glob2rx(paste("*tmax_pred*",date_selected,"*_365d_GAM_CAI3_all_10302012.rst",sep="")) #Search for files in relation to fusion
|
187 |
|
|
lf_cai3<-list.files(pattern=file_pat) #Search for files in relation to fusion
|
188 |
|
|
|
189 |
|
|
cai3_mod<-load_obj("results_mod_obj__365d_GAM_CAI3_all_10302012.RData")
|
190 |
|
|
names(cai3_mod$gam_CAI_mod[[1]]) #check names within object
|
191 |
|
|
formulas<-as.character(cai3_mod$gam_CAI_mod[[1]]$formulas)
|
192 |
|
|
rmse<-(cai3_mod$gam_CAI_mod[[1]]$tb_metrics1)
|
193 |
|
|
nmodels<-6
|
194 |
|
|
for(i in 4:(nmodels+3)){ # start of the for loop #1
|
195 |
|
|
rmse[,i]<-as.numeric(as.character(rmse[,i]))
|
196 |
|
|
}
|
197 |
|
|
rmse_val<-as.numeric(c(rmse[1,9],rmse[1,4:8]))
|
198 |
|
|
rmse_val<-format(rmse_val,digits=3)
|
199 |
|
|
rmse_val<-paste("(",rmse_val,")",sep="")
|
200 |
|
|
rname<-c("CAI_Kr",formulas)
|
201 |
|
|
rname<-paste(rname,rmse_val,sep=" ")
|
202 |
|
|
|
203 |
|
|
rast_cai3<-stack(lf_cai3)
|
204 |
|
|
rast_cai3<-mask(rast_cai3,mask_ELEV_SRTM)
|
205 |
|
|
#s_range <- c(min(minValue(diff)), max(maxValue(diff)))
|
206 |
|
|
s_range<-c(minValue(rast_cai3),maxValue(rast_cai3)) #stack min and max
|
207 |
|
|
s_range<-c(min(s_range),max(s_range))
|
208 |
|
|
col_breaks <- pretty(s_range, n=50)
|
209 |
|
|
lab_breaks <- pretty(s_range, n=5)
|
210 |
|
|
temp_colors <- colorRampPalette(c('blue', 'white', 'red'))
|
211 |
|
|
par(mfrow=c(2,3))
|
212 |
|
|
for (k in 1:length(lf_cai3)){
|
213 |
|
|
cai3_r<-raster(rast_cai3,k)
|
214 |
|
|
plot(cai3_r, breaks=col_breaks, col=temp_colors(length(col_breaks)-1),
|
215 |
|
|
axis=list(at=lab_breaks, labels=lab_breaks),main=rname[k])
|
216 |
|
|
}
|
217 |
|
|
|
218 |
|
|
#Wihtout setting range
|
219 |
|
|
s_range<-c(-12,18)
|
220 |
|
|
#col_breaks <- pretty(s_range, n=50)
|
221 |
|
|
#lab_breaks <- pretty(s_range, n=5)
|
222 |
|
|
col_breaks<-49
|
223 |
|
|
temp_colors <- colorRampPalette(c('blue', 'white', 'red'))
|
224 |
|
|
par(mfrow=c(2,3))
|
225 |
|
|
for (k in 1:length(lf_cai3)){
|
226 |
|
|
cai3_r<-raster(rast_cai3,k)
|
227 |
|
|
plot(cai3_r, col=temp_colors(col_breaks),
|
228 |
|
|
,main=rname[k])
|
229 |
|
|
}
|
230 |
|
|
|
231 |
|
|
###### PART 2: ################
|
232 |
|
|
#### EXAMINING LST AND ELEVATION DATASETS TO CHECK FOR EXTREMES AND SELECT SCREENING VALUES....
|
233 |
|
|
|
234 |
|
|
# Examining ELEV_SRTM and screening out...
|
235 |
|
|
|
236 |
|
|
pos<-match("ELEV_SRTM",layerNames(s_raster))
|
237 |
|
|
ELEV_SRTM<-raster(s_raster,pos)
|
238 |
|
|
head(freq(ELEV_SRTM)) #Show how many pixels are below 0
|
239 |
|
|
|
240 |
|
|
# mask out values below zero
|
241 |
|
|
elev<-ELEV_SRTM
|
242 |
|
|
elev[0<elev]<-NA #Remove all negative elevation
|
243 |
|
|
#Show how many were NA...
|
244 |
|
|
|
245 |
|
|
mask_ELEV_SRTM<-elev >0
|
246 |
|
|
quartz()
|
247 |
|
|
plot(elev, main="Elevation in Oregon")
|
248 |
|
|
savePlot(paste("elevation_Oregon",out_prefix,".png", sep=""), type="png")
|
249 |
|
|
|
250 |
|
|
# Examining LST data: extremes values and number of valid observation
|
251 |
|
|
|
252 |
|
|
l<-list.files(pattern="mean_month.*rescaled.rst")
|
253 |
|
|
ln<-list.files(pattern="mean_month_valid_obs_.*_Sum.rst") #number of observation missing
|
254 |
|
|
l<-mixedsort(l)
|
255 |
|
|
ln<-mixedsort(ln)
|
256 |
|
|
|
257 |
|
|
#l <-readLines(paste(path,"/",infile6, sep=""))
|
258 |
|
|
molst<-stack(l) #Creating a raster stack...
|
259 |
|
|
nobslst<-stack(ln) #Creating a raster stack...
|
260 |
|
|
|
261 |
|
|
#setwd(old)
|
262 |
|
|
molst<-molst-273.16 #K->C #LST stack of monthly average...
|
263 |
|
|
idx <- seq(as.Date('2010-01-15'), as.Date('2010-12-15'), 'month')
|
264 |
|
|
molst <- setZ(molst, idx)
|
265 |
|
|
layerNames(molst) <- month.abb
|
266 |
|
|
layerNames(nobslst)<-month.abb
|
267 |
|
|
|
268 |
|
|
X11(width=18,height=12)
|
269 |
|
|
#PLOT THE STACK AND SAVE IMAGE
|
270 |
|
|
plot(molst,col=temp.colors(25))
|
271 |
|
|
savePlot("lst_climatology_OR.png",type="png")
|
272 |
|
|
levelplot(molst,col.regions=temp.colors(25))
|
273 |
|
|
plot(nobslst)
|
274 |
|
|
savePlot("lst_climatology_OR_nobs.png",type="png")
|
275 |
|
|
dev.off)
|
276 |
|
|
|
277 |
c352e9e1
|
Benoit Parmentier
|
#note differences in patternin agricultural areas and
|
278 |
2a0267ad
|
Benoit Parmentier
|
extremeValues(molst) #not yet implemented...
|
279 |
|
|
min_values<-cellStats(molst,"min")
|
280 |
|
|
max_values<-cellStats(molst,"max")
|
281 |
|
|
mean_values<-cellStats(molst,"mean")
|
282 |
|
|
sd_values<-cellStats(molst,"sd")
|
283 |
|
|
#median_values<-cellStats(molst,"median") Does not extist
|
284 |
|
|
statistics_LST_s<-cbind(min_values,max_values,mean_values,sd_values) #This shows that some values are extremes...especially in October
|
285 |
|
|
LST_stat_data<-as.data.frame(statistics_LST_s)
|
286 |
|
|
names(LST_stat_data)<-c("min","max","mean","sd")
|
287 |
|
|
# Statistics for number of valid observation stack
|
288 |
|
|
min_values<-cellStats(nobslst,"min")
|
289 |
|
|
max_values<-cellStats(nobslst,"max")
|
290 |
|
|
mean_values<-cellStats(nobslst,"mean")
|
291 |
|
|
sd_values<-cellStats(nobslst,"sd")
|
292 |
|
|
statistics_LSTnobs_s<-cbind(min_values,max_values,mean_values,sd_values) #This shows that some values are extremes...especially in October
|
293 |
|
|
LSTnobs_stat_data<-as.data.frame(statistics_LSTnobs_s)
|
294 |
|
|
|
295 |
|
|
X11(width=12,height=12)
|
296 |
|
|
#Plot statiscs (mean,min,max) for monthly LST images
|
297 |
|
|
plot(1:12,LST_stat_data$mean,type="b",ylim=c(-15,60),col="black",xlab="month",ylab="tmax (degree C)")
|
298 |
|
|
lines(1:12,LST_stat_data$min,type="b",col="blue")
|
299 |
|
|
lines(1:12,LST_stat_data$max,type="b",col="red")
|
300 |
|
|
text(1:12,LST_stat_data$mean,rownames(LST_stat_data),cex=1,pos=2)
|
301 |
|
|
|
302 |
|
|
legend("topleft",legend=c("min","mean","max"), cex=1.5, col=c("blue","black","red"),
|
303 |
|
|
lty=1)
|
304 |
|
|
title(paste("LST statistics for Oregon", "2010",sep=" "))
|
305 |
|
|
savePlot("lst_statistics_OR.png",type="png")
|
306 |
|
|
|
307 |
|
|
#Plot number of valid observations for LST
|
308 |
|
|
plot(1:12,LSTnobs_stat_data$mean,type="b",ylim=c(0,280),col="black",xlab="month",ylab="tmax (degree C)")
|
309 |
|
|
lines(1:12,LSTnobs_stat_data$min,type="b",col="blue")
|
310 |
|
|
lines(1:12,LSTnobs_stat_data$max,type="b",col="red")
|
311 |
|
|
text(1:12,LSTnobs_stat_data$mean,rownames(LSTnobs_stat_data),cex=1,pos=2)
|
312 |
|
|
|
313 |
|
|
legend("topleft",legend=c("min","mean","max"), cex=1.5, col=c("blue","black","red"),
|
314 |
|
|
lty=1)
|
315 |
|
|
title(paste("LST number of valid observations for Oregon", "2010",sep=" "))
|
316 |
|
|
savePlot("lst_nobs_OR.png",type="png")
|
317 |
|
|
|
318 |
|
|
mm_10<-raster(molst,10)
|
319 |
|
|
tab_mm_10<-freq(mm_10) #This shows that there are only 6 pixels with values between -86 and -89.. and 8 pixels with values between 38 and 31
|
320 |
|
|
head(tab_mm_10)
|
321 |
|
|
#screening maybe good...?
|
322 |
|
|
|
323 |
|
|
mm_01<-raster(molst,1)
|
324 |
|
|
tab_mm_10<-freq(mm_01) #This shows that there are only 6 pixels with values between -86 and -89.. and 8 pixels with values between 38 and 31
|
325 |
|
|
#head(tab_mm_10)
|
326 |
|
|
|
327 |
|
|
mm_07<-raster(molst,7)
|
328 |
|
|
tab_mm_10<-freq(mm_07) #This shows that there are only 6 pixels with values between -86 and -89.. and 8 pixels with values between 38 and 31
|
329 |
|
|
#head(tab_mm_10)
|
330 |
|
|
|
331 |
|
|
LST_comp<-stack(mm_01,mm_07)
|
332 |
|
|
|
333 |
|
|
X11(width=18,height=10)
|
334 |
|
|
#PLOT THE STACK AND SAVE IMAGE
|
335 |
|
|
plot(LST_comp,col=temp.colors(25))
|
336 |
|
|
#title("LST spatial pattern in January and August")
|
337 |
|
|
savePlot("lst_climatology_OR.png",type="png")
|
338 |
|
|
|
339 |
|
|
|
340 |
|
|
###### PART 3: EXAMANING MODEL OUTPUTS USING MOD OBJECTS TO CHECK AND COMPARE RESULTS ############
|
341 |
|
|
############ LOOKING AT SPECIFIC STATIONS...
|
342 |
|
|
path_data_cai<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_CAI"
|
343 |
|
|
path_data_fus<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_GAM"
|
344 |
|
|
setwd(path_data) #output data
|
345 |
|
|
gam_fus<-load_obj(file.path(path_data_fus,
|
346 |
|
|
"results2_fusion_Assessment_measure_all_365d_GAM_fusion_all_lstd_10242012.RData"))
|
347 |
|
|
gam_cai<-load_obj(file.path(path_dat_cai,
|
348 |
|
|
"results2_CAI_Assessment_measure_all_365d_GAM_CAI2_all_lstd_10262012.RData") #This contains all the info
|
349 |
|
|
sampling_obj_cai<-load_obj(file.path(path_data_cai,
|
350 |
|
|
"results2_CAI_sampling_obj_365d_GAM_CAI2_all_lstd_10262012.RData"))
|
351 |
|
|
sampling_date_list<-sampling_obj_cai$sampling_dat$date
|
352 |
|
|
date_selected="20100103"
|
353 |
|
|
|
354 |
|
|
k<-match(date_selected,sampling_date_list)
|
355 |
|
|
|
356 |
|
|
fus_d<-gam_fus[[k]] #object for the first date...20100103
|
357 |
|
|
fus_d$tb_metrics1 #This is a summary for validation metrics...(RMSE and MAE)
|
358 |
|
|
|
359 |
|
|
cai_d<-gam_cai[[k]] #object for the first date...20100103
|
360 |
|
|
cai_d$tb_metrics1 #This is a summary for validation metrics...(RMSE and MAE)
|
361 |
|
|
|
362 |
|
|
names(fus_d$mod_obj) #list models, nine models: mod1, mod2,...,mod9a, mod9b
|
363 |
|
|
names(cai_d$mod_obj)
|
364 |
|
|
|
365 |
|
|
mod1_f<-fus_d$mod_obj$mod1 # extract model 1 for fusion on
|
366 |
|
|
mod4_f<-fus_d$mod_obj$mod4
|
367 |
|
|
mod7_f<-fus_d$mod_obj$mod7
|
368 |
|
|
data_sf<-fus_d$data_s
|
369 |
|
|
data_vf<-fus_d$data_v
|
370 |
|
|
|
371 |
|
|
plot(mod1_f) #Examining the fit
|
372 |
|
|
plot(mod1) #Examing the fit
|
373 |
|
|
#data_monthf<-fus_d1$data_month
|
374 |
|
|
|
375 |
|
|
mod1<-cai_d1$mod_obj$mod1
|
376 |
|
|
mod4<-cai_d1$mod_obj$mod4
|
377 |
|
|
mod7<-cai_d1$mod_obj$mod7
|
378 |
|
|
data_s<-cai_d1$data_s
|
379 |
|
|
data_v<-cai_d1$data_v
|
380 |
|
|
data_month<-cai_d1$data_month
|
381 |
|
|
|
382 |
|
|
#### COMPARE LST AND TMAX averages over the year... ###############
|
383 |
|
|
#### Screening necessary?
|
384 |
|
|
|
385 |
|
|
LSTD_bias_month<-vector("list",length(gam_cai)) # list to store bias information
|
386 |
|
|
data_month_list<-vector("list",length(gam_cai)) # list to store data_month used for training
|
387 |
|
|
|
388 |
|
|
for (i in 1:length(gam_cai)){
|
389 |
|
|
data_month_list[[i]]<-gam_cai[[i]]$data_month
|
390 |
|
|
#LSTD_bias_month[[i]]<-aggregate(LSTD_bias~month,data=data_month_list[[i]],mean)
|
391 |
|
|
}
|
392 |
|
|
|
393 |
|
|
#tmp<-do.call(rbind,LSTD_bias_month)
|
394 |
|
|
data_month_all<-do.call(rbind,data_month_list)
|
395 |
|
|
|
396 |
|
|
data_month_all$LST<-data_month_all$LST-273.16
|
397 |
|
|
LSTD_bias_avgm<-aggregate(LSTD_bias~month,data=data_month_all,mean)
|
398 |
|
|
LSTD_bias_sdm<-aggregate(LSTD_bias~month,data=data_month_all,sd)
|
399 |
|
|
|
400 |
|
|
LST_avgm<-aggregate(LST~month,data=data_month_all,mean)
|
401 |
|
|
TMax_avgm<-aggregate(TMax~month,data=data_month_all,mean)
|
402 |
|
|
|
403 |
|
|
#plot(LST_avgm,type="b")
|
404 |
|
|
#lines(TMax_avgm,col="blue",add=T)
|
405 |
|
|
plot(data_month_all$month,data_month_all$LST, xlab="month",ylab="LST at station")
|
406 |
|
|
title(paste("Monthly LST at stations in Oregon", "2010",sep=" "))
|
407 |
|
|
savePlot(paste("LST_at_stations_per_month_",out_prefix,".png", sep=""), type="png")
|
408 |
|
|
|
409 |
|
|
tab_freq<-cbind(hist_LST$breaks,hist_LST$counts)
|
410 |
|
|
median(data_month_all$LST,na.rm=T)
|
411 |
|
|
median(data_month_all$TMax,na.rm=T)
|
412 |
|
|
|
413 |
|
|
statistics_LST_s<-as.data.frame(statistics_LST_s)
|
414 |
|
|
plot(TMax_avgm,type="b",ylim=c(0,35),col="red",xlab="month",ylab="tmax (degree C)")
|
415 |
|
|
lines(1:12,statistics_LST_s$mean_values,type="b",col="blue")
|
416 |
|
|
text(TMax_avgm[,1],TMax_avgm[,2],rownames(statistics_LST_s),cex=0.8,pos=2)
|
417 |
|
|
|
418 |
|
|
legend("topleft",legend=c("TMax","LST"), cex=1, col=c("red","blue"),
|
419 |
|
|
lty=1)
|
420 |
|
|
title(paste("Monthly mean tmax and LST at stations in Oregon", "2010",sep=" "))
|
421 |
|
|
savePlot(paste("Temporal_profile_res",id[i],out_prefix,".png", sep=""), type="png")
|
422 |
|
|
|
423 |
|
|
#Account for forest
|
424 |
|
|
data_month_all_forest<-data_month_all[data_month_all$LC1>=50,]
|
425 |
|
|
data_month_all_grass<-data_month_all[data_month_all$LC3>=10,]
|
426 |
|
|
LST_avgm_forest<-aggregate(LST~month,data=data_month_all_forest,mean)
|
427 |
|
|
LST_avgm_grass<-aggregate(LST~month,data=data_month_all_grass,mean)
|
428 |
|
|
|
429 |
|
|
plot(TMax_avgm,type="b",ylim=c(-7,42))
|
430 |
|
|
lines(LST_avgm,col="blue",add=T)
|
431 |
|
|
lines(LST_avgm_forest,col="green",add=T)
|
432 |
|
|
lines(LST_avgm_grass,col="red",add=T)
|
433 |
|
|
legend("topleft",legend=c("TMax","LST","LST_forest","LST_grass"), cex=0.8, col=c("black","blue","green","red"),
|
434 |
|
|
lty=1)
|
435 |
|
|
title(paste("Monthly average tmax for stations in Oregon ", "2010",sep=""))
|
436 |
|
|
|
437 |
|
|
savePlot(paste("Temporal_profile_res",id[i],out_prefix,".png", sep=""), type="png")
|
438 |
|
|
|
439 |
|
|
#Check where forest >50 and grass>10
|
440 |
|
|
#Create mask using land cover data
|
441 |
|
|
pos<-match("LC1",layerNames(s_raster)) #Find the layer which contains water bodies
|
442 |
|
|
LC1<-subset(s_raster,pos)
|
443 |
|
|
LC1[is.na(LC1)]<-0 #Since NA values are 0, we assign all zero to NA
|
444 |
|
|
LC1_50<-LC1>= 50
|
445 |
|
|
plot(LC1_50)
|
446 |
|
|
pos<-match("LC3",layerNames(s_raster)) #Find the layer which contains water bodies
|
447 |
|
|
LC3<-subset(s_raster,pos)
|
448 |
|
|
LC3[is.na(LC3)]<-0 #Since NA values are 0, we assign all zero to NA
|
449 |
|
|
LC3_10<-LC3>= 10
|
450 |
|
|
|
451 |
|
|
LC1_mean <- cellStats(LC3, mean)
|
452 |
|
|
LC1_mean <- cellStats(LC3, mean)
|
453 |
|
|
#
|
454 |
|
|
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)
|
455 |
|
|
rclmat<-matrix(avl,ncol=3,byrow=TRUE)
|
456 |
|
|
LC3_rec<-reclass(LC3,rclmat) #Loss of layer names when using reclass
|
457 |
|
|
tab_freq<-freq(LC3_rec)
|
458 |
|
|
|
459 |
|
|
crosstab(LC1_50,LC3_10) #Very little overalp between classes...
|
460 |
|
|
X11(12,12)
|
461 |
|
|
LC_rec<-stack(LC1_50,LC3_10)
|
462 |
|
|
layerNames(LC_rec)<-c("forest: LC1>50%","grass: LC3>10%")
|
463 |
|
|
plot(LC_rec,col=c("black","red"))
|
464 |
|
|
dev.off()
|
465 |
|
|
#legend("topleft",legend=c("fo","LST","LST_forest","LST_grass"), cex=0.8, col=c("black","blue","green","red"),
|
466 |
|
|
#lty=1)
|
467 |
|
|
|
468 |
c352e9e1
|
Benoit Parmentier
|
###################### Visualization of results ######################
|
469 |
2a0267ad
|
Benoit Parmentier
|
|
470 |
|
|
"tmax_predicted*20100103_30_1_365d_GAM_fusion_all_lstd_10242012.rst"
|
471 |
|
|
oldpath<-getwd()
|
472 |
|
|
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_10242012_GAM"
|
473 |
|
|
setwd(path)
|
474 |
|
|
date_selected<-"20100103"
|
475 |
|
|
#lf_fus<-list.files(pattern=paste("GAM_bias_tmax_predicted_mod8_",date_selected,"_30_1_365d_GAM_fusion_all_lstd_10242012.rst$",sep="")) #Search for files in relation to fusion
|
476 |
|
|
file_pat<-glob2rx(paste("*tmax_pred*",date_selected,"*_365d_GAM_fusion_all_lstd_10242012.rst",sep="")) #Search for files in relation to fusion
|
477 |
|
|
lf_fus1a<-list.files(pattern=file_pat) #Search for files in relation to fusion
|
478 |
|
|
|
479 |
|
|
rast_fus1a<-stack(lf_fus1a)
|
480 |
|
|
rast_fus1a<-mask(rast_fus1a,mask_ELEV_SRTM)
|
481 |
|
|
|
482 |
|
|
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations_07192012_GAM"
|
483 |
|
|
setwd(path)
|
484 |
|
|
date_selected<-"20100103"
|
485 |
|
|
#lf_fus<-list.files(pattern=paste("GAM_bias_tmax_predicted_mod8_",date_selected,"_30_1_365d_GAM_fusion_all_lstd_10242012.rst$",sep="")) #Search for files in relation to fusion
|
486 |
|
|
#lf_fus1<-list.files(pattern=paste("*.tmax_predicted.*.",date_selected,".*._365d_GAM_fusion_lstd_10062012.rst$",sep=""))
|
487 |
|
|
file_pat<-glob2rx(paste("*tmax_predicted*",date_selected,"*_365d_GAM_fusion_lstd_10062012.rst",sep="")) #Search for files in relation to fusion
|
488 |
|
|
lf_fus1s<-list.files(pattern=file_pat) #Search for files in relation to fusion
|
489 |
|
|
|
490 |
|
|
rast_fus1s<-stack(lf_fus1s)
|
491 |
|
|
rast_fus1s<-mask(rast_fus1s,mask_ELEV_SRTM)
|
492 |
|
|
|
493 |
|
|
s_range<-c(minValue(rast_fusa1),maxValue(rast_fusa1)) #stack min and max
|
494 |
|
|
s_range<-c(min(s_range),max(s_range))
|
495 |
|
|
col_breaks <- pretty(s_range, n=50)
|
496 |
|
|
lab_breaks <- pretty(s_range, n=5)
|
497 |
|
|
temp_colors <- colorRampPalette(c('blue', 'white', 'red'))
|
498 |
|
|
X11(width=18,height=12)
|
499 |
|
|
par(mfrow=c(3,3))
|
500 |
|
|
for (k in 1:length(lf_fus1a)){
|
501 |
|
|
fus1a_r<-raster(rast_fus1a,k)
|
502 |
|
|
plot(fus1a_r, breaks=col_breaks, col=temp_colors(length(col_breaks)-1),
|
503 |
|
|
axis=list(at=lab_breaks, labels=lab_breaks))
|
504 |
|
|
}
|
505 |
|
|
|
506 |
|
|
#Wihtout setting range
|
507 |
|
|
s_range<-c(-12,18)
|
508 |
|
|
#col_breaks <- pretty(s_range, n=50)
|
509 |
|
|
#lab_breaks <- pretty(s_range, n=5)
|
510 |
|
|
col_breaks<-49
|
511 |
|
|
temp_colors <- colorRampPalette(c('blue', 'white', 'red'))
|
512 |
|
|
lab_breaks <- pretty(s_range, n=5)
|
513 |
|
|
X11(width=18,height=12)
|
514 |
|
|
par(mfrow=c(3,3))
|
515 |
|
|
mod_name<-paste("mod",1:8, sep="")
|
516 |
|
|
rname<-c("FUS_kr",mod_name)
|
517 |
|
|
for (k in 1:length(lf_fus1a)){
|
518 |
|
|
fus1a_r<-raster(rast_fus1a,k)
|
519 |
|
|
plot(fus1a_r, breaks=col_breaks, col=temp_colors(col_breaks-1),
|
520 |
|
|
axis=list(at=lab_breaks, labels=lab_breaks),main=rname[k])
|
521 |
|
|
}
|
522 |
|
|
|
523 |
|
|
#Wihtout setting range
|
524 |
|
|
s_range<-c(-12,23)
|
525 |
|
|
#col_breaks <- pretty(s_range, n=50)
|
526 |
|
|
#lab_breaks <- pretty(s_range, n=5)
|
527 |
|
|
col_breaks<-49
|
528 |
|
|
temp_colors <- colorRampPalette(c('blue', 'white', 'red'))
|
529 |
|
|
lab_breaks <- pretty(s_range, n=5)
|
530 |
|
|
X11(width=18,height=12)
|
531 |
|
|
par(mfrow=c(3,3))
|
532 |
|
|
mod_name<-paste("mod",1:8, sep="")
|
533 |
|
|
rname<-c("FUS_kr",mod_name)
|
534 |
|
|
for (k in 1:length(lf_fus1s)){
|
535 |
|
|
fus1s_r<-raster(rast_fus1s,k)
|
536 |
|
|
plot(fus1s_r, breaks=col_breaks, col=temp_colors(col_breaks-1),
|
537 |
|
|
axis=list(at=lab_breaks, labels=lab_breaks),main=rname[k])
|
538 |
|
|
}
|
539 |
|
|
|
540 |
|
|
|
541 |
|
|
|
542 |
|
|
|