1
|
################## Covariates and stations analyses #######################################
|
2
|
############################ Covariate production for a given tile/region ##########################################
|
3
|
#This script examines inputs and outputs from the interpolation step.
|
4
|
#Part 1: Script produces summary information about stations used in the interpolation
|
5
|
#Part 2: Script produces plots of input covariates for study region
|
6
|
#AUTHOR: Benoit Parmentier
|
7
|
#DATE: 03/27/2013
|
8
|
|
9
|
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#???--
|
10
|
|
11
|
##Comments and TODO:
|
12
|
#Add interpolation analyses of number of stations used for every day of the year...
|
13
|
|
14
|
## Input arguments
|
15
|
#The output is a list of four shapefile names produced by the function:
|
16
|
#1) loc_stations: locations of stations as shapefile in EPSG 4326
|
17
|
#2) loc_stations_ghcn: ghcn daily data for the year range of interpolation (locally projected)
|
18
|
#3) daily_query_ghcn_data: ghcn daily data from daily query before application of quality flag
|
19
|
#4) daily_covar_ghcn_data: ghcn daily data with covariates for the year range of interpolation (locally projected)
|
20
|
#5) monthly_query_ghcn_data: ghcn daily data from monthly query before application of quality flag
|
21
|
#6) monthly_covar_ghcn_data: ghcn monthly averaged data with covariates for the year range of interpolation (locally projected)
|
22
|
#7) var
|
23
|
#8) infile_covariate: s_raster brick file
|
24
|
#9) covar_names: variable mames
|
25
|
#10) raster_prediction
|
26
|
#11) world_countries
|
27
|
#12) region_outline
|
28
|
#13) out_prefix
|
29
|
|
30
|
#list_param <-list_outfiles$loc_stations
|
31
|
#list_param <-list_outfiles$loc_stations_ghcn
|
32
|
#list_param <-list_outfiles$loc_stations_ghcn_data
|
33
|
#list_param <-list_outfiles$loc_stations
|
34
|
#list_param <-list_outfiles$loc_stations
|
35
|
|
36
|
## Functions used in the script
|
37
|
|
38
|
load_obj <- function(f)
|
39
|
{
|
40
|
env <- new.env()
|
41
|
nm <- load(f, env)[1]
|
42
|
env[[nm]]
|
43
|
}
|
44
|
|
45
|
extract_number_obs<-function(list_param){
|
46
|
|
47
|
method_mod_obj<-list_param$method_mod_obj
|
48
|
#Change to results_mod_obj[[i]]$data_s to make it less specific
|
49
|
lapply(1:length(method_obj),function(k) nrow(method_mod_obj[[k]]$data_s))
|
50
|
lapply(1:length(method_obj),function(k) nrow(method_mod_obj[[k]]$data_v))
|
51
|
lapply(1:length(clim_obj),function(k) nrow(method_mod_obj[[k]]$data_month))
|
52
|
#number of observations
|
53
|
|
54
|
return()
|
55
|
}
|
56
|
|
57
|
#### STEP 1: read in data
|
58
|
|
59
|
stat_reg <- readOGR(dsn=dirname(list_outfiles$loc_stations),sub(".shp","",basename(list_outfiles$loc_stations)))
|
60
|
data_reg <- readOGR(dsn=dirname(list_outfiles$loc_stations_ghcn),sub(".shp","",basename(list_outfiles$loc_stations_ghcn)))
|
61
|
data_RST_SDF <-readOGR(dsn=dirname(list_outfiles$daily_covar_ghcn_data),sub(".shp","",basename(list_outfiles$daily_covar_ghcn_data)))
|
62
|
data_m <- readOGR(dsn=dirname(list_outfiles$monthly_query_ghcn_data),sub(".shp","",basename(list_outfiles$monthly_query_ghcn_data)))
|
63
|
dst<-readOGR(dsn=dirname(list_outfiles$monthly_covar_ghcn_data),sub(".shp","",basename(list_outfiles$monthly_covar_ghcn_data)))
|
64
|
|
65
|
s_raster<-brick(infile_covariates)
|
66
|
names(s_raster)<-covar_names
|
67
|
rast_ref<-subset(s_raster,"mm_01") # mean month for January
|
68
|
|
69
|
######## PART I
|
70
|
### Figue 1: stations in the study area/processing tile
|
71
|
png(paste("Total_number_of_stations_in_study_area_",out_prefix,".png", sep=""))
|
72
|
plot(rast_ref)
|
73
|
plot(stat_reg,add=TRUE)
|
74
|
nb_point1<-paste("n_station=",length(stat_reg$STAT_ID))
|
75
|
#Add the number of data points on the plot
|
76
|
legend("topleft",legend=c(nb_point1),bty="n",cex=1.2)
|
77
|
dev.off()
|
78
|
|
79
|
#add number of stations+ name of region
|
80
|
#title()
|
81
|
|
82
|
### Figue 2: stations in the study area/processing tile: from query without flag screening
|
83
|
png(paste("Studay_area_",out_prefix,".png", sep=""))
|
84
|
plot(rast_ref)
|
85
|
plot(data_d, add=TRUE)
|
86
|
nb_point<-paste("ns=",length(data_d$TMax),sep="")
|
87
|
nb_point2<-paste("ns_obs=",length(data_s$TMax)-sum(is.na(data_s[[y_var_name]])),sep="")
|
88
|
nb_point3<-paste("n_month=",length(data_month$TMax),sep="")
|
89
|
#Add the number of data points on the plot
|
90
|
legend("topleft",legend=c(nb_point,nb_point2,nb_point3),bty="n",cex=0.8)
|
91
|
dev.off()
|
92
|
|
93
|
### Figue 3: stations in the study area/processing tile: after screening
|
94
|
plot(rast_ref)
|
95
|
plot(data_reg,add=TRUE)
|
96
|
|
97
|
### Figue 4: stations in the study area/processing tile: from monthly climatology query without flag screening
|
98
|
plot(rast_ref)
|
99
|
plot(data_RST_SDF,add=TRUE)
|
100
|
|
101
|
### Figue 5: stations in the study area/processing tile: after screening
|
102
|
plot(rast_ref)
|
103
|
plot(data_m,add=TRUE)
|
104
|
|
105
|
### Figue 6: histogram
|
106
|
#plot(dst)
|
107
|
|
108
|
#dst$nobs_station
|
109
|
dst$nobs_station<-dst$nbs_stt
|
110
|
hist(dst$nobs_station)
|
111
|
|
112
|
### Figue 7: LST and TMax
|
113
|
|
114
|
names_tmp<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12")
|
115
|
LST_s<-subset(s_raster,names_tmp)
|
116
|
names(LST_s)<-names_tmp
|
117
|
names_tmp<-c("nobs_01","nobs_02","nobs_03","nobs_04","nobs_05","nobs_06","nobs_07","nobs_08",
|
118
|
"nobs_09","nobs_10","nobs_11","nobs_12")
|
119
|
LST_nobs<-subset(s_raster,names_tmp)
|
120
|
|
121
|
## Function...note differnces in patternin agricultural areas
|
122
|
|
123
|
list_statistic_values_LST_s<-mclapply(c("min","max","mean","sd"),FUN=cellStats,x=LST_s,mc.preschedule=FALSE,mc.cores = 4)
|
124
|
list_statistic_values_LST_nobs<-mclapply(c("min","max","mean","sd"),FUN=cellStats,x=LST_nobs,mc.preschedule=FALSE,mc.cores = 4)
|
125
|
|
126
|
statistics_LST_s<-do.call(cbind,list_statistic_values_LST_s)
|
127
|
statistics_LST_nobs<-do.call(cbind,list_statistic_values_LST_nobs)
|
128
|
LST_stat_data<-as.data.frame(statistics_LST_s)
|
129
|
names(LST_stat_data)<-c("min","max","mean","sd")
|
130
|
statistics_LSTnobs_s<-cbind(min_values,max_values,mean_values,sd_values) #This shows that some values are extremes...especially in October
|
131
|
LSTnobs_stat_data<-as.data.frame(statistics_LSTnobs_s)
|
132
|
names(LSTnobs_stat_data)<-c("min","max","mean","sd")
|
133
|
|
134
|
# X11(width=12,height=12)
|
135
|
# #Plot statiscs (mean,min,max) for monthly LST images
|
136
|
plot(1:12,LST_stat_data$mean,type="b",ylim=c(-15,70),col="black",xlab="month",ylab="tmax (degree C)")
|
137
|
lines(1:12,LST_stat_data$min,type="b",col="blue")
|
138
|
lines(1:12,LST_stat_data$max,type="b",col="red")
|
139
|
text(1:12,LST_stat_data$mean,rownames(LST_stat_data),cex=1,pos=2)
|
140
|
legend("topleft",legend=c("min","mean","max"), cex=0.8, col=c("blue","black","red"),
|
141
|
lty=1,lwd=1.4,bty="n")
|
142
|
title(paste("LST statistics for Study area",sep=" "))
|
143
|
# savePlot("lst_statistics_OR.png",type="png")
|
144
|
|
145
|
#### Fig8...#####
|
146
|
|
147
|
# #Plot number of valid observations for LST
|
148
|
plot(1:12,LSTnobs_stat_data$mean,type="b",ylim=c(0,280),col="black",xlab="month",ylab="tmax (degree C)")
|
149
|
lines(1:12,LSTnobs_stat_data$min,type="b",col="blue")
|
150
|
lines(1:12,LSTnobs_stat_data$max,type="b",col="red")
|
151
|
text(1:12,LSTnobs_stat_data$mean,rownames(LSTnobs_stat_data),cex=1,pos=2)
|
152
|
#
|
153
|
legend("topleft",legend=c("min","mean","max"), cex=1.5, col=c("blue","black","red"),lty=1)
|
154
|
title(paste("LST number of valid observations for Oregon", "2010",sep=" "))
|
155
|
# savePlot("lst_nobs_OR.png",type="png")
|
156
|
#
|
157
|
#### Fig 9...#####
|
158
|
|
159
|
#Use data_month!!!
|
160
|
d_month<-aggregate(TMax~month, data=dst, mean) #Calculate monthly mean for every station in OR
|
161
|
d_month<-aggregate(TMax~month, data=dst, length) #Calculate monthly mean for every station in OR
|
162
|
|
163
|
plot(d_month,type="l")
|
164
|
# plot(data_month$TMax,add=TRUE)
|
165
|
|
166
|
#data_month #-> plot for everymonth
|
167
|
|
168
|
#number of stations
|
169
|
## Summarize information for the day: write out textfile...
|
170
|
|
171
|
#Number of station per month
|
172
|
#Number of station per day (training, testing,NA)
|
173
|
#metrics_v,metrics_s
|
174
|
#
|
175
|
|
176
|
# ################
|
177
|
# #PART 2: Region Covariate analyses ###
|
178
|
# ################
|
179
|
#
|
180
|
# # This should be in a separate script to analyze covariates from region.
|
181
|
#
|
182
|
# #MAP1:Study area with LC mask and tiles/polygon outline
|
183
|
|
184
|
#LC_mask<-subset(s_raster,"LC12")
|
185
|
#LC_mask[LC_mask==100]<-NA
|
186
|
#LC_mask <- LC_mask < 100
|
187
|
#LC_mask_rec<-LC_mask
|
188
|
#LC_mask_rec[is.na(LC_mask_rec)]<-0
|
189
|
|
190
|
#Add proportion covered by study area+ total of image pixels
|
191
|
#tmp_tb<-freq(LC_mask_rec)
|
192
|
#tmp_tb[2,2]/sum(tmp_tb[,2])*100
|
193
|
#png(paste("Study_area_",
|
194
|
# out_prefix,".png", sep=""))
|
195
|
#plot(LC_mask_rec,legend=FALSE,col=c("black","red"))
|
196
|
#legend("topright",legend=c("Outside","Inside"),title="Study area",
|
197
|
# pt.cex=0.9,fill=c("black","red"),bty="n")
|
198
|
# title("Study area")
|
199
|
# dev.off()
|
200
|
|
201
|
#
|
202
|
# #MAP 2: plotting land cover in the study region:
|
203
|
#
|
204
|
# l1<-"LC1,Evergreen/deciduous needleleaf trees"
|
205
|
# l2<-"LC2,Evergreen broadleaf trees"
|
206
|
# l3<-"LC3,Deciduous broadleaf trees"
|
207
|
# l4<-"LC4,Mixed/other trees"
|
208
|
# l5<-"LC5,Shrubs"
|
209
|
# l6<-"LC6,Herbaceous vegetation"
|
210
|
# l7<-"LC7,Cultivated and managed vegetation"
|
211
|
# l8<-"LC8,Regularly flooded shrub/herbaceous vegetation"
|
212
|
# l9<-"LC9,Urban/built-up"
|
213
|
# l10<-"LC10,Snow/ice"
|
214
|
# l11<-"LC11,Barren lands/sparse vegetation"
|
215
|
# l12<-"LC12,Open water"
|
216
|
# lc_names_str<-c(l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12)
|
217
|
#
|
218
|
# names(lc_reg_s)<-lc_names_str
|
219
|
#
|
220
|
# png(paste("LST_TMax_scatterplot_",sampling_dat$date[i],"_",sampling_dat$prop[i],"_",sampling_dat$run_samp[i], out_prefix,".png", sep=""))
|
221
|
# plot(modst$TMax,sta_tmax_from_lst,xlab="Station mo Tmax",ylab="LST mo Tmax",main=paste("LST vs TMax for",datelabel,sep=" "))
|
222
|
# abline(0,1)
|
223
|
# nb_point<-paste("n=",length(modst$TMax),sep="")
|
224
|
# mean_bias<-paste("LST bigrasas= ",format(mean(modst$LSTD_bias,na.rm=TRUE),digits=3),sep="")
|
225
|
# #Add the number of data points on the plot
|
226
|
# legend("topleft",legend=c(mean_bias,nb_point),bty="n")
|
227
|
# dev.off()
|
228
|
#
|
229
|
#
|
230
|
### MAP3: Majority land cover for every pixels in the study region
|
231
|
|
232
|
|
233
|
|
234
|
###Map 4: Elevation and LST in January
|
235
|
# tmp_s<-stack(LST,elev_1)
|
236
|
# png(paste("LST_elev_",sampling_dat$date[i],"_",sampling_dat$prop[i],"_",sampling_dat$run_samp[i], out_prefix,".png", sep=""))
|
237
|
# plot(tmp_s)
|
238
|
#
|
239
|
###Histogram elevation?
|
240
|
#
|
241
|
# #Map 5: mean TMIN/TMAX: LST climatology per month
|
242
|
#
|
243
|
# names_tmp<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12")
|
244
|
# LST_s<-subset(s_raster,names_tmp)
|
245
|
# names_tmp<-c("nobs_01","nobs_02","nobs_03","nobs_04","nobs_05","nobs_06","nobs_07","nobs_08",
|
246
|
# "nobs_09","nobs_10","nobs_11","nobs_12")
|
247
|
# LST_nobs<-subset(s_raster,names_tmp)
|
248
|
#
|
249
|
# #Map 6: number of obs- LST climatology per month
|
250
|
#
|
251
|
# names_tmp<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12")
|
252
|
# LST_s<-subset(s_raster,names_tmp)
|
253
|
# names_tmp<-c("nobs_01","nobs_02","nobs_03","nobs_04","nobs_05","nobs_06","nobs_07","nobs_08",
|
254
|
# "nobs_09","nobs_10","nobs_11","nobs_12")
|
255
|
# LST_nobs<-subset(s_raster,names_tmp)
|
256
|
#
|
257
|
|
258
|
# LST_nobs<-mask(LST_nobs,LC_mask,filename="test2.tif")
|
259
|
# LST_s<-mask(LST_s,LC_mask,filename="test3.tif")
|
260
|
# c("Jan","Feb")
|
261
|
# plot(LST_s)
|
262
|
# plot(LST_nobs)
|
263
|
#
|
264
|
# #Map 7: LST-TMAX/TMIN fit for all 12 months...
|
265
|
#
|
266
|
#
|
267
|
#
|
268
|
#### End of script ####
|