Revision 2ba35352
Added by Benoit Parmentier about 11 years ago
| climate/research/oregon/interpolation/results_covariates_database_stations_output_analyses.R | ||
|---|---|---|
| 4 | 4 |
#Part 1: Script produces summary information about stations used in the interpolation |
| 5 | 5 |
#Part 2: Script produces plots of input covariates for study region |
| 6 | 6 |
#AUTHOR: Benoit Parmentier |
| 7 |
#DATE: 03/27/2013
|
|
| 7 |
#DATE: 04/01/2013
|
|
| 8 | 8 |
|
| 9 | 9 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#???-- |
| 10 | 10 |
|
| ... | ... | |
| 20 | 20 |
#5) monthly_query_ghcn_data: ghcn daily data from monthly query before application of quality flag |
| 21 | 21 |
#6) monthly_covar_ghcn_data: ghcn monthly averaged data with covariates for the year range of interpolation (locally projected) |
| 22 | 22 |
#7) var |
| 23 |
#8) range_years |
|
| 24 |
#9) range_years_clim |
|
| 23 | 25 |
#8) infile_covariate: s_raster brick file |
| 24 | 26 |
#9) covar_names: variable mames |
| 25 | 27 |
#10) raster_prediction |
| ... | ... | |
| 27 | 29 |
#12) region_outline |
| 28 | 30 |
#13) out_prefix |
| 29 | 31 |
|
| 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 | 32 |
## Functions used in the script |
| 37 | 33 |
|
| 38 | 34 |
load_obj <- function(f) |
| ... | ... | |
| 43 | 39 |
} |
| 44 | 40 |
|
| 45 | 41 |
extract_number_obs<-function(list_param){
|
| 42 |
#Function to extract the number of observations used in modeling |
|
| 46 | 43 |
|
| 47 | 44 |
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)) |
|
| 45 |
|
|
| 46 |
#Change to results_mod_obj[[i]]$data_s to make it less specific!!!! |
|
| 47 |
|
|
| 52 | 48 |
#number of observations |
| 49 |
list_nrow_data_s<-lapply(1:length(method_obj),function(k) nrow(method_mod_obj[[k]]$data_s)) |
|
| 50 |
list_nrow_data_v<-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 |
list_nrow_data_month<-lapply(1:length(gamclim_fus_mod),function(k) nrow(gamclim_fus_mod[[k]]$data_month)) |
|
| 53 |
|
|
| 54 |
#number of valid observations |
|
| 55 |
list_nrow_valid_data_s<-lapply(1:length(method_obj),function(k) length(method_mod_obj[[k]]$data_s$[[y_var_name]]) |
|
| 56 |
list_nrow_data_v<-lapply(1:length(method_obj),function(k) length(method_mod_obj[[k]]$data_v$[[y_var_name]]) |
|
| 57 |
list_nrow_valid_data_month<-lapply(1:length(gamclim_fus_mod),function(k) length(gamclim_fus_mod[[k]]$data_month$y_var)) |
|
| 58 |
|
|
| 59 |
c1<-do.call(rbind,list_nrow_data_s) |
|
| 60 |
c2<-do.call(rbind,list_nrow_data_v) |
|
| 61 |
c3<-do.call(rbind,list_nrow_data_month) |
|
| 53 | 62 |
|
| 54 |
return() |
|
| 63 |
c1v<-do.call(rbind,list_valid_nrow_data_s) |
|
| 64 |
c2v<-do.call(rbind,list_valid_nrow_data_v) |
|
| 65 |
c3v<-do.call(rbind,list_valid_nrow_data_month) |
|
| 66 |
|
|
| 67 |
n_data_s<-cbdind(c1,c1v) |
|
| 68 |
n_data_v<-cbdind(c2,c2v) |
|
| 69 |
n_data_month<-cbdind(c3,c3vv) |
|
| 70 |
list_observations<-list(n_data_s,n_data_v,n_data_month) |
|
| 71 |
return(list_observations) |
|
| 55 | 72 |
} |
| 56 | 73 |
|
| 74 |
########################################################################### |
|
| 75 |
########################## BEGIN SCRIPT/FUNCTION ########################## |
|
| 76 |
|
|
| 77 |
|
|
| 78 |
######################################## |
|
| 57 | 79 |
#### STEP 1: read in data |
| 58 | 80 |
|
| 81 |
#Stations in the processing region/study area |
|
| 59 | 82 |
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)))
|
|
| 83 |
#Stations available before screening the data query: ghcn daily data for the year range of interpolation |
|
| 84 |
data_d <- readOGR(dsn=dirname(list_outfiles$loc_stations_ghcn),sub(".shp","",basename(list_outfiles$loc_stations_ghcn)))
|
|
| 85 |
#Stations available after screening the data query: ghcn daily data for the year range of interpolation |
|
| 86 |
data_reg <- readOGR(dsn=dirname(list_outfiles$daily_query_ghcn_data),sub(".shp","",basename(list_outfiles$daily_query_ghcn_data)))
|
|
| 87 |
#Covariates data available after screening:ghcn daily data with covariates for the year range of interpolation |
|
| 61 | 88 |
data_RST_SDF <-readOGR(dsn=dirname(list_outfiles$daily_covar_ghcn_data),sub(".shp","",basename(list_outfiles$daily_covar_ghcn_data)))
|
| 89 |
#Stations before screening monthly_query_ghcn_data: ghcn daily data from monthly query before application of quality flag |
|
| 62 | 90 |
data_m <- readOGR(dsn=dirname(list_outfiles$monthly_query_ghcn_data),sub(".shp","",basename(list_outfiles$monthly_query_ghcn_data)))
|
| 91 |
#Stations after screening monthly_query_ghcn_data, extraction of covariates and monthly averages |
|
| 63 | 92 |
dst<-readOGR(dsn=dirname(list_outfiles$monthly_covar_ghcn_data),sub(".shp","",basename(list_outfiles$monthly_covar_ghcn_data)))
|
| 64 | 93 |
|
| 94 |
### Load data used in fitting the model at monthly scale... |
|
| 95 |
data_month<-gamclim_fus_mod[[1]]$data_month |
|
| 96 |
list_data_month<-lapply(1:length(gamclim_fus_mod),function(k) gamclim_fus_mod[[k]]$data_month) |
|
| 97 |
|
|
| 98 |
#Loading covariates raster images |
|
| 65 | 99 |
s_raster<-brick(infile_covariates) |
| 66 | 100 |
names(s_raster)<-covar_names |
| 67 | 101 |
rast_ref<-subset(s_raster,"mm_01") # mean month for January |
| 68 | 102 |
|
| 69 |
######## PART I |
|
| 103 |
names(raster_prediction_obj) |
|
| 104 |
var<-list_param$var |
|
| 105 |
|
|
| 106 |
raster_prediction_obj<-load_obj(list_param$raster_prediction_obj) |
|
| 107 |
#method_mod_obj<-raster_prediction_obj$method_mod_obj |
|
| 108 |
method_mod_obj<-raster_prediction_obj$gam_fus_mod #change later for any model type |
|
| 109 |
#validation_obj<-raster_prediction_obj$validation_obj |
|
| 110 |
validation_obj<-raster_prediction_obj$gam_fus_validation_mod #change later for any model type |
|
| 111 |
#clim_obj<-raster_prediction_obj$clim_obj |
|
| 112 |
clim_obj<-raster_prediction_obj$gamclim_fus_mod #change later for any model type |
|
| 113 |
|
|
| 114 |
names(raster_prediction_obj$method_obj[[1]]) |
|
| 115 |
data_s<-validation_obj[[1]]$data_s |
|
| 116 |
summary_data_v<-validation_obj[[1]]$summary_data_v |
|
| 117 |
names(validation_obj[[1]]) |
|
| 118 |
################################################################### |
|
| 119 |
######## PART I: Script produces summary information about stations used in the interpolation ######## |
|
| 120 |
|
|
| 70 | 121 |
### Figue 1: stations in the study area/processing tile |
| 71 | 122 |
png(paste("Total_number_of_stations_in_study_area_",out_prefix,".png", sep=""))
|
| 72 | 123 |
plot(rast_ref) |
| 73 | 124 |
plot(stat_reg,add=TRUE) |
| 74 | 125 |
nb_point1<-paste("n_station=",length(stat_reg$STAT_ID))
|
| 75 | 126 |
#Add the number of data points on the plot |
| 127 |
title("Stations located in the study area")
|
|
| 76 | 128 |
legend("topleft",legend=c(nb_point1),bty="n",cex=1.2)
|
| 77 | 129 |
dev.off() |
| 78 | 130 |
|
| 79 |
#add number of stations+ name of region |
|
| 80 |
#title() |
|
| 81 |
|
|
| 82 | 131 |
### Figue 2: stations in the study area/processing tile: from query without flag screening |
| 83 |
png(paste("Studay_area_",out_prefix,".png", sep=""))
|
|
| 132 |
png(paste("Stations_for_range_",range_years[1],"_",range_years[2],"_no_screening",out_prefix,".png", sep=""))
|
|
| 84 | 133 |
plot(rast_ref) |
| 85 | 134 |
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="")
|
|
| 135 |
nb_point<-paste("nrow=",nrow(data_d),sep="")
|
|
| 136 |
nb_point2<-paste("ns_stations=",length(unique(data_d$station)),sep="")
|
|
| 89 | 137 |
#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)
|
|
| 138 |
legend("topleft",legend=c(nb_point,nb_point2),bty="n",cex=1)
|
|
| 139 |
title(paste("Stations available for year ",range_years[1],sep=""))
|
|
| 91 | 140 |
dev.off() |
| 92 | 141 |
|
| 93 |
### Figue 3: stations in the study area/processing tile: after screening |
|
| 142 |
### Figue 3: stations in the study area/processing tile: after screening |
|
| 143 |
png(paste("Stations_for_range_",range_years[1],"_",range_years[2],"_after_screening",out_prefix,".png", sep=""))
|
|
| 94 | 144 |
plot(rast_ref) |
| 95 | 145 |
plot(data_reg,add=TRUE) |
| 146 |
nb_point<-paste("nrow=",nrow(data_reg),sep="")
|
|
| 147 |
nb_point2<-paste("ns_stations=",length(unique(data_reg$station)),sep="")
|
|
| 148 |
#Add the number of data points on the plot |
|
| 149 |
legend("topleft",legend=c(nb_point,nb_point2),bty="n",cex=1)
|
|
| 150 |
title(paste("Stations available for year ",range_years[1]," after screening",sep=""))
|
|
| 151 |
#Add the number of data points on the plot |
|
| 152 |
dev.off() |
|
| 96 | 153 |
|
| 97 |
### Figue 4: stations in the study area/processing tile: from monthly climatology query without flag screening |
|
| 154 |
### Figue 4: stations in the study area/processing tile: after screening and covar extraction |
|
| 155 |
png(paste("Stations_for_range_",range_years[1],"_",range_years[2],"_after_screening",out_prefix,".png", sep=""))
|
|
| 98 | 156 |
plot(rast_ref) |
| 99 | 157 |
plot(data_RST_SDF,add=TRUE) |
| 100 |
|
|
| 101 |
### Figue 5: stations in the study area/processing tile: after screening |
|
| 158 |
nb_point<-paste("nrow=",nrow(data_RST_SDF),sep="")
|
|
| 159 |
nb_point2<-paste("ns_stations=",length(unique(data_RST_SDF$station)),sep="")
|
|
| 160 |
#Add the number of data points on the plot |
|
| 161 |
legend("topleft",legend=c(nb_point,nb_point2),bty="n",cex=1)
|
|
| 162 |
title(paste("Stations year ",range_years[1]," after screening and covar extraction",sep=""))
|
|
| 163 |
#Add the number of data points on the plot |
|
| 164 |
dev.off() |
|
| 165 |
|
|
| 166 |
### Figue 5: stations in the study area/processing tile: monthly query for specified range of years before screening |
|
| 167 |
png(paste("Stations_monthly_for_range_",range_years[1],"_",range_years[2],"_before_screening",out_prefix,".png", sep=""))
|
|
| 102 | 168 |
plot(rast_ref) |
| 103 | 169 |
plot(data_m,add=TRUE) |
| 170 |
nb_point<-paste("nrow=",nrow(data_m),sep="")
|
|
| 171 |
nb_point2<-paste("ns_stations=",length(unique(data_m$station)),sep="")
|
|
| 172 |
#Add the number of data points on the plot |
|
| 173 |
legend("topleft",legend=c(nb_point,nb_point2),bty="n",cex=1)
|
|
| 174 |
title(paste("Stations ",range_years[1],"-",range_years[2]," after screening and covar extraction",sep=""))
|
|
| 175 |
#Add the number of data points on the plot |
|
| 176 |
dev.off() |
|
| 104 | 177 |
|
| 105 | 178 |
### Figue 6: histogram |
| 106 |
#plot(dst) |
|
| 107 | 179 |
|
| 108 | 180 |
#dst$nobs_station |
| 181 |
png(paste("Stations_data_month_modeled_for_range_",range_years_clim[1],"_",range_years_clim[2],out_prefix,".png", sep=""))
|
|
| 109 | 182 |
dst$nobs_station<-dst$nbs_stt |
| 110 | 183 |
hist(dst$nobs_station) |
| 184 |
dev.off() |
|
| 111 | 185 |
|
| 112 |
### Figue 7: LST and TMax |
|
| 186 |
### Figure 7: data month |
|
| 187 |
|
|
| 188 |
png(paste("Stations_data_month_modeled_for_range_",range_years_clim[1],"_",range_years_clim[2],out_prefix,".png", sep=""))
|
|
| 189 |
par(mfrow=c(3,4)) |
|
| 190 |
for (j in 1:12){
|
|
| 191 |
|
|
| 192 |
data_month<-list_data_month[[j]] |
|
| 193 |
plot(rast_ref) |
|
| 194 |
plot(data_month,add=TRUE) |
|
| 195 |
nb_point<-paste("nrow=",nrow(data_month),sep="")
|
|
| 196 |
nb_point2<-paste("ns_stations=",length(unique(data_month$station)),sep="")
|
|
| 197 |
nb_point3<-paste("ns_non_na_stations=",length(unique(data_month$y_var)),sep="")
|
|
| 198 |
#Add the number of data points on the plot |
|
| 199 |
legend("topleft",legend=c(nb_point,nb_point2,nb_point3),bty="n",cex=0.9)
|
|
| 200 |
title(paste("Stations ",range_years_clim[1],"-",range_years_clim[2]," used for modeling on ",sep=""))
|
|
| 201 |
#data_NA<-subset(data=data_month,is.na(data_month$y_var)) |
|
| 202 |
data_NA<-data_month[is.na(data_month$y_var),] |
|
| 203 |
plot(data_NA,add=TRUE,pch=2,col=c("red"))
|
|
| 204 |
legend("topright",legend=c("NA value"),pch=2,col=c("red"),bty="n",cex=0.9)
|
|
| 205 |
} |
|
| 206 |
dev.off() |
|
| 207 |
|
|
| 208 |
### Figue 8: LST and TMax |
|
| 113 | 209 |
|
| 114 | 210 |
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 | 211 |
LST_s<-subset(s_raster,names_tmp) |
| ... | ... | |
| 117 | 213 |
names_tmp<-c("nobs_01","nobs_02","nobs_03","nobs_04","nobs_05","nobs_06","nobs_07","nobs_08",
|
| 118 | 214 |
"nobs_09","nobs_10","nobs_11","nobs_12") |
| 119 | 215 |
LST_nobs<-subset(s_raster,names_tmp) |
| 120 |
|
|
| 121 |
## Function...note differnces in patternin agricultural areas |
|
| 122 |
|
|
| 216 |
|
|
| 123 | 217 |
list_statistic_values_LST_s<-mclapply(c("min","max","mean","sd"),FUN=cellStats,x=LST_s,mc.preschedule=FALSE,mc.cores = 4)
|
| 124 | 218 |
list_statistic_values_LST_nobs<-mclapply(c("min","max","mean","sd"),FUN=cellStats,x=LST_nobs,mc.preschedule=FALSE,mc.cores = 4)
|
| 125 | 219 |
|
| ... | ... | |
| 131 | 225 |
LSTnobs_stat_data<-as.data.frame(statistics_LSTnobs_s) |
| 132 | 226 |
names(LSTnobs_stat_data)<-c("min","max","mean","sd")
|
| 133 | 227 |
|
| 134 |
# X11(width=12,height=12) |
|
| 135 |
# #Plot statiscs (mean,min,max) for monthly LST images |
|
| 228 |
png(paste("Stations_data_month_modeled_for_range_",range_years_clim[1],"_",range_years_clim[2],out_prefix,".png", sep=""))
|
|
| 136 | 229 |
plot(1:12,LST_stat_data$mean,type="b",ylim=c(-15,70),col="black",xlab="month",ylab="tmax (degree C)") |
| 137 | 230 |
lines(1:12,LST_stat_data$min,type="b",col="blue") |
| 138 | 231 |
lines(1:12,LST_stat_data$max,type="b",col="red") |
| ... | ... | |
| 141 | 234 |
lty=1,lwd=1.4,bty="n") |
| 142 | 235 |
title(paste("LST statistics for Study area",sep=" "))
|
| 143 | 236 |
# savePlot("lst_statistics_OR.png",type="png")
|
| 144 |
|
|
| 145 |
#### Fig8...##### |
|
| 237 |
dev.off() |
|
| 238 |
|
|
| 239 |
#### Fig9...##### |
|
| 146 | 240 |
|
| 147 | 241 |
# #Plot number of valid observations for LST |
| 242 |
png(paste("Stations_data_month_modeled_for_range_",range_years_clim[1],"_",range_years_clim[2],out_prefix,".png", sep=""))
|
|
| 148 | 243 |
plot(1:12,LSTnobs_stat_data$mean,type="b",ylim=c(0,280),col="black",xlab="month",ylab="tmax (degree C)") |
| 149 | 244 |
lines(1:12,LSTnobs_stat_data$min,type="b",col="blue") |
| 150 | 245 |
lines(1:12,LSTnobs_stat_data$max,type="b",col="red") |
| 151 | 246 |
text(1:12,LSTnobs_stat_data$mean,rownames(LSTnobs_stat_data),cex=1,pos=2) |
| 152 |
# |
|
| 153 | 247 |
legend("topleft",legend=c("min","mean","max"), cex=1.5, col=c("blue","black","red"),lty=1)
|
| 154 | 248 |
title(paste("LST number of valid observations for Oregon", "2010",sep=" "))
|
| 155 | 249 |
# savePlot("lst_nobs_OR.png",type="png")
|
| 156 |
# |
|
| 157 |
#### Fig 9...##### |
|
| 250 |
dev.off() |
|
| 251 |
|
|
| 252 |
#### Fig 10...##### |
|
| 158 | 253 |
|
| 159 | 254 |
#Use data_month!!! |
| 255 |
|
|
| 160 | 256 |
d_month<-aggregate(TMax~month, data=dst, mean) #Calculate monthly mean for every station in OR |
| 161 | 257 |
d_month<-aggregate(TMax~month, data=dst, length) #Calculate monthly mean for every station in OR |
| 162 | 258 |
|
| ... | ... | |
| 229 | 325 |
# |
| 230 | 326 |
### MAP3: Majority land cover for every pixels in the study region |
| 231 | 327 |
|
| 232 |
|
|
| 328 |
#Add barplot of majority map... |
|
| 233 | 329 |
|
| 234 | 330 |
###Map 4: Elevation and LST in January |
| 235 | 331 |
# tmp_s<-stack(LST,elev_1) |
Also available in: Unified diff
results stations and covar outputs modifications