Revision 9c848d7a
Added by Benoit Parmentier over 12 years ago
| climate/research/oregon/interpolation/GAM_LST_reg.R | ||
|---|---|---|
| 22 | 22 |
path<-"H:/Data/IPLANT_project/data_Oregon_stations" |
| 23 | 23 |
setwd(path) |
| 24 | 24 |
infile2<-"dates_interpolation_03052012.txt" #List of 10 dates for the regression |
| 25 |
infile2<-"list_365_dates_04212012.txt" |
|
| 26 | 25 |
prop<-0.3 #Proportion of testing retained for validation |
| 27 |
out_prefix<-"_04212012_LST"
|
|
| 26 |
out_prefix<-"_04252012_LST_residuals"
|
|
| 28 | 27 |
infile3<-"LST_dates_var_names.txt" |
| 29 | 28 |
infile4<-"models_interpolation_04032012b.txt" |
| 30 | 29 |
|
| ... | ... | |
| 33 | 32 |
###Reading the station data and setting up for models' comparison |
| 34 | 33 |
filename<-sub(".shp","",infile1) #Removing the extension from file.
|
| 35 | 34 |
ghcn<-readOGR(".", filename) #reading shapefile
|
| 36 |
|
|
| 35 |
|
|
| 36 |
CRS<-proj4string(ghcn) |
|
| 37 |
|
|
| 38 |
|
|
| 37 | 39 |
ghcn = transform(ghcn,Northness = cos(ASPECT)) #Adding a variable to the dataframe |
| 38 | 40 |
ghcn = transform(ghcn,Eastness = sin(ASPECT)) #adding variable to the dataframe. |
| 39 | 41 |
ghcn = transform(ghcn,Northness_w = sin(slope)*cos(ASPECT)) #Adding a variable to the dataframe |
| ... | ... | |
| 50 | 52 |
results_GCV<- matrix(1,length(dates),length(models)+2) |
| 51 | 53 |
results_DEV<- matrix(1,length(dates),length(models)+2) |
| 52 | 54 |
results_RMSE<- matrix(1,length(dates),length(models)+2) |
| 53 |
cor_LST_LC1<-matrix(1,length(dates),1) #correlation LST-LC1
|
|
| 54 |
cor_LST_LC3<-matrix(1,length(dates),1) #correlation LST-LC3
|
|
| 55 |
cor_LST_tmax<-matrix(1,length(dates),1) #correlation LST-tmax
|
|
| 55 |
cor_LST_LC1<-matrix(1,10,1) #correlation LST-LC1
|
|
| 56 |
cor_LST_LC3<-matrix(1,10,1) #correlation LST-LC3
|
|
| 57 |
cor_LST_tmax<-matrix(1,10,1) #correlation LST-tmax
|
|
| 56 | 58 |
#Screening for bad values |
| 57 | 59 |
|
| 58 | 60 |
ghcn_all<-ghcn |
| 59 | 61 |
ghcn_test<-subset(ghcn,ghcn$tmax>-150 & ghcn$tmax<400) |
| 60 | 62 |
ghcn_test2<-subset(ghcn_test,ghcn_test$ELEV_SRTM>0) |
| 61 | 63 |
ghcn<-ghcn_test2 |
| 64 |
#coords<- ghcn[,c('x_OR83M','y_OR83M')]
|
|
| 62 | 65 |
|
| 63 | 66 |
month_var<-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")
|
| 64 | 67 |
ghcn.subsets <-lapply(dates, function(d) subset(ghcn, date==d)) #this creates a list of 10 subsets data |
| ... | ... | |
| 91 | 94 |
####Regression part 2: GAM models |
| 92 | 95 |
|
| 93 | 96 |
mod1<- gam(tmax~ s(lat) + s (lon) + s (ELEV_SRTM), data=data_s) |
| 94 |
#mod2<- gam(tmax~ s(lat,lon,ELEV_SRTM), data=data_s) |
|
| 95 |
mod2<- gam(tmax~ s(lat,lon) + s(ELEV_SRTM), data=data_s) |
|
| 97 |
mod2<- gam(tmax~ s(lat,lon,ELEV_SRTM), data=data_s) |
|
| 96 | 98 |
mod3<- gam(tmax~ s(lat) + s (lon) + s (ELEV_SRTM) + s (Northness)+ s (Eastness) + s(DISTOC), data=data_s) |
| 97 | 99 |
mod4<- gam(tmax~ s(lat) + s (lon) + s(ELEV_SRTM) + s(Northness) + s (Eastness) + s(DISTOC) + s(LST), data=data_s) |
| 98 | 100 |
mod5<- gam(tmax~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST), data=data_s) |
| ... | ... | |
| 165 | 167 |
results_RMSE[i,7]<- RMSE_mod5 |
| 166 | 168 |
results_RMSE[i,8]<- RMSE_mod6 |
| 167 | 169 |
results_RMSE[i,9]<- RMSE_mod7 |
| 168 |
#Saving dataset in dataframes |
|
| 170 |
|
|
| 171 |
#Saving dataset in dataframes: residuals from RMSE |
|
| 172 |
|
|
| 173 |
# data_v$mod1<-y_mod1$fit |
|
| 174 |
# data_v$mod2<-y_mod2$fit |
|
| 175 |
# data_v$mod3<-y_mod3$fit |
|
| 176 |
# data_v$mod4<-y_mod4$fit |
|
| 177 |
# data_v$mod5<-y_mod5$fit |
|
| 178 |
# data_v$mod6<-y_mod6$fit |
|
| 179 |
# data_v$mod7<-y_mod7$fit |
|
| 180 |
# |
|
| 181 |
# data_s$mod1<-mod1$fit |
|
| 182 |
# data_s$mod2<-mod2$fit |
|
| 183 |
# data_s$mod3<-mod3$fit |
|
| 184 |
# data_s$mod4<-mod4$fit |
|
| 185 |
# data_s$mod5<-mod5$fit |
|
| 186 |
# data_s$mod6<-mod6$fit |
|
| 187 |
# data_s$mod7<-mod7$fit |
|
| 188 |
# |
|
| 189 |
data_v$res_mod1<-as.numeric(res_mod1) |
|
| 190 |
data_v$res_mod2<-as.numeric(res_mod2) |
|
| 191 |
data_v$res_mod3<-as.numeric(res_mod3) |
|
| 192 |
data_v$res_mod4<-as.numeric(res_mod4) |
|
| 193 |
data_v$res_mod5<-as.numeric(res_mod5) |
|
| 194 |
data_v$res_mod6<-as.numeric(res_mod6) |
|
| 195 |
data_v$res_mod7<-as.numeric(res_mod7) |
|
| 196 |
|
|
| 197 |
data_s$res_mod1<-as.numeric(mod1$residuals) |
|
| 198 |
data_s$res_mod2<-as.numeric(mod2$residuals) |
|
| 199 |
data_s$res_mod3<-as.numeric(mod3$residuals) |
|
| 200 |
data_s$res_mod4<-as.numeric(mod4$residuals) |
|
| 201 |
data_s$res_mod5<-as.numeric(mod5$residuals) |
|
| 202 |
data_s$res_mod6<-as.numeric(mod6$residuals) |
|
| 203 |
data_s$res_mod7<-as.numeric(mod7$residuals) |
|
| 204 |
|
|
| 169 | 205 |
data_name<-paste("ghcn_v_",dates[[i]],sep="")
|
| 170 | 206 |
assign(data_name,data_v) |
| 207 |
write.table(data_v, file= paste(path,"/",data_name,".txt",sep=""), sep=" ") |
|
| 208 |
#write out a new shapefile (including .prj component) |
|
| 209 |
coords<- data_v[,c('x_OR83M','y_OR83M')]
|
|
| 210 |
coordinates(data_v)<-coords |
|
| 211 |
proj4string(data_v)<-CRS #Need to assign coordinates... |
|
| 212 |
outfile<-sub(".shp","",data_name) #Removing extension if it is present
|
|
| 213 |
writeOGR(data_v,".", outfile, driver ="ESRI Shapefile") |
|
| 214 |
|
|
| 171 | 215 |
data_name<-paste("ghcn_s_",dates[[i]],sep="")
|
| 172 | 216 |
assign(data_name,data_s) |
| 173 |
#ghcn_v<-ls(pattern="ghcn_v_") |
|
| 217 |
write.table(data_s, file= paste(path,"/",data_name,".txt",sep=""), sep=" ") |
|
| 218 |
#write out a new shapefile (including .prj component) |
|
| 219 |
coords<- data_s[,c('x_OR83M','y_OR83M')]
|
|
| 220 |
coordinates(data_s)<-coords |
|
| 221 |
proj4string(data_s)<-CRS #Need to assign coordinates.. |
|
| 222 |
outfile<-sub(".shp","",data_name) #Removing extension if it is present
|
|
| 223 |
writeOGR(data_s,".", outfile, driver ="ESRI Shapefile") |
|
| 174 | 224 |
|
| 175 | 225 |
# end of the for loop #1 |
| 176 | 226 |
} |
| ... | ... | |
| 200 | 250 |
f<-file(paste(path,"/","results_GAM_Assessment",out_prefix,"2.txt",sep=""),"w") |
| 201 | 251 |
write(results_table_RMSE, file=f) |
| 202 | 252 |
close(f) |
| 253 |
|
|
| 203 | 254 |
# End of script########## |
| 204 | 255 |
|
| 205 | 256 |
#Selecting dates and files based on names |
Also available in: Unified diff
GAM LST, added section to save training and test residuals in shapefiles, task #361