Project

General

Profile

« Previous | Next » 

Revision 0847d47c

Added by Benoit Parmentier over 11 years ago

Initial commit summary outputs for covariates and stations used in modeling

View differences:

climate/research/oregon/interpolation/results_covariates_database_stations_output_analyses.R
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 ####

Also available in: Unified diff