Revision 15a39f25
Added by Benoit Parmentier over 12 years ago
| climate/research/oregon/interpolation/gwr_reg.R | ||
|---|---|---|
| 1 | 1 |
####################GWR of Tmax for 10 dates.##################### |
| 2 |
#This script generate station values for the Oregon case study. This program loads the station data from a shp file |
|
| 3 |
#and performs a GWR regression.
|
|
| 4 |
#Script created by Benoit Parmentier on March 13, 2012.
|
|
| 2 |
#This script generates predicted values from station values for the Oregon case study. This program loads the station data from a shp file
|
|
| 3 |
#and performs a GWR regression and a kriged surface of residuals.
|
|
| 4 |
#Script created by Benoit Parmentier on April 3, 2012.
|
|
| 5 | 5 |
|
| 6 | 6 |
###Loading r library and packages |
| 7 | 7 |
library(sp) |
| ... | ... | |
| 9 | 9 |
library(rgdal) |
| 10 | 10 |
library(spgwr) |
| 11 | 11 |
library(gpclib) |
| 12 |
library(PBSmapping)
|
|
| 12 |
data
|
|
| 13 | 13 |
library(maptools) |
| 14 | 14 |
library(gstat) |
| 15 | 15 |
###Parameters and arguments |
| ... | ... | |
| 101 | 101 |
plot(v) |
| 102 | 102 |
tryCatch(v.fit<-fit.variogram(v,vgm(1,"Sph", 150000,1)),error=function()next) |
| 103 | 103 |
gwr_res_krige<-krige(residuals~1, data_s,mean_LST, v.fit)#mean_LST provides the data grid/raster image for the kriging locations. |
| 104 |
|
|
| 105 |
# GWR visualization of Residuals using histograms and over space |
|
| 106 |
X11() |
|
| 107 |
title=paste("Histogram of residuals of ",data_name, sep="")
|
|
| 108 |
hist(data_s$residuals,main=title) |
|
| 109 |
savePlot(paste("Histogram_",data_name,out_prefix,".png", sep=""), type="png")
|
|
| 110 |
dev.off() |
|
| 104 | 111 |
|
| 105 |
# GWR visualization of Residuals fit over space |
|
| 106 |
grays = gray.colors(5,0.45, 0.95) |
|
| 107 |
image(gwr_res_krige,col=grays) #needs to change to have a bipolar palette !!! |
|
| 108 |
|
|
| 109 |
#image(mean_LST, col=grays,breaks = c(185,245,255,275,315,325)) |
|
| 112 |
X11(width=20,height=20) |
|
| 113 |
topo = cm.colors(9) |
|
| 114 |
image(gwr_res_krige,col=topo) #needs to change to have a bipolar palette !!! |
|
| 110 | 115 |
|
| 111 | 116 |
plot(OR_state, axes = TRUE, add=TRUE) |
| 112 | 117 |
plot(data_s, pch=1, col="red", cex= abs(data_s$residuals)/10, add=TRUE) #Taking the absolute values because residuals are |
| 113 |
LegVals<- c(0,20,40,80,110) |
|
| 114 |
legend("topleft", legend=LegVals,pch=1,col="red",pt.cex=LegVals/10,bty="n",title= "residuals")
|
|
| 115 |
legend("left", legend=c("275-285","285-295","295-305", "305-315","315-325"),fill=grays, bty="n", title= "LST mean DOY=244")
|
|
| 118 |
LegVals<- c(0,10,20,30,40,50,110) |
|
| 119 |
legend(-98000,510000, legend=LegVals,pch=1,col="red",pt.cex=LegVals/10,bty="n",title= "residuals",cex=1.6) |
|
| 120 |
#legend("left", legend=c("275-285","285-295","295-305", "305-315","315-325"),fill=grays, bty="n", title= "LST mean DOY=244")
|
|
| 121 |
legend(-98000,290000, legend=c("-60 -30","-30 -20","-20 -10", "-10 0"," 0 10"," 10 20", " 20 30"," 30 60"),fill=topo, bty="n", title= "Kriged RMSE",cex=1.6)
|
|
| 122 |
title(paste("Kriging of residuals of ",data_name, sep=""),cex=2)
|
|
| 123 |
|
|
| 124 |
krig_raster_name<-paste("Kriged_res_",data_name,out_prefix,".tif", sep="")
|
|
| 125 |
#writeGDAL(gwr_res_krige,fname=krig_raster_name, driver="GTiff", type="Float32",options ="INTERLEAVE=PIXEL") |
|
| 126 |
|
|
| 127 |
savePlot(paste("Kriged_res_",data_name,out_prefix,".png", sep=""), type="png")
|
|
| 128 |
dev.off() |
|
| 116 | 129 |
|
| 117 |
savePlot(paste(data_name,out_prefix,".png", sep=""), type="png") |
|
| 118 |
dev.off() |
|
| 119 | 130 |
} |
| 120 | 131 |
|
| 121 | 132 |
## Plotting and saving diagnostic measures |
| ... | ... | |
| 134 | 145 |
|
| 135 | 146 |
# ############################### |
| 136 | 147 |
|
| 148 |
# # GWR visualization of Residuals fit over space |
|
| 149 |
# X11() |
|
| 150 |
# title=paste("Histogram of residuals of ",data_name, sep="")
|
|
| 151 |
# hist(data_s$residuals,main=title) |
|
| 152 |
# savePlot(paste("Histogram_",data_name,out_prefix,".png", sep=""), type="png")
|
|
| 153 |
# dev.off() |
|
| 154 |
# |
|
| 155 |
# |
|
| 156 |
data_s<-ghcn_s20100901 |
|
| 157 |
X11(width=20,height=20) |
|
| 158 |
topo = cm.colors(9) |
|
| 159 |
image(gwr_res_krige,col=topo) #needs to change to have a bipolar palette !!! |
|
| 137 | 160 |
|
| 138 |
|
|
| 161 |
#image(mean_LST, col=grays,breaks = c(185,245,255,275,315,325)) |
|
| 162 |
|
|
| 163 |
plot(OR_state, axes = TRUE, add=TRUE) |
|
| 164 |
plot(data_s, pch=1, col="red", cex= abs(data_s$residuals)/10, add=TRUE) #Taking the absolute values because residuals are |
|
| 165 |
LegVals<- c(0,10,20,30,40,50,110) |
|
| 166 |
legend(-98000,510000, legend=LegVals,pch=1,col="red",pt.cex=LegVals/10,bty="n",title= "residuals",cex=1.6) |
|
| 167 |
#legend("left", legend=c("275-285","285-295","295-305", "305-315","315-325"),fill=grays, bty="n", title= "LST mean DOY=244")
|
|
| 168 |
legend(-98000,290000, legend=c("-60 -30","-30 -20","-20 -10", "-10 0"," 0 10"," 10 20", " 20 30"," 30 60"),fill=topo, bty="n", title= "Kriged RMSE",cex=1.6)
|
|
| 169 |
title(paste("Kriging of residuals of ",data_name, sep=""),cex=2)
|
|
| 170 |
|
|
| 171 |
krig_raster_name<-paste("Kriged_res_",data_name,out_prefix,".rst", sep="")
|
|
| 172 |
writeGDAL(gwr_res_krige,fname="test_krige.tif", driver="GTiff", type="Float32",options ="INTERLEAVE=PIXEL") |
|
| 173 |
|
|
| 174 |
savePlot(paste("Kriged_res_",data_name,out_prefix,".png", sep=""), type="png")
|
|
| 175 |
dev.off() |
|
| 176 |
|
|
| 177 |
|
|
| 178 |
|
|
| 179 |
|
|
| 139 | 180 |
#Compare the coefficients and residuals using both 30 and 100% |
| 140 | 181 |
#coefficients are stored in gwrG$SDF$lon |
| 141 | 182 |
#write out a new shapefile (including .prj component) |
Also available in: Unified diff
GWR modified code to save plots and images, task #364