Project

General

Profile

« Previous | Next » 

Revision defe53e5

Added by Benoit Parmentier over 12 years ago

GWR assessment by visualizing residuals using kriging

View differences:

climate/research/oregon/interpolation/gwr_reg.R
8 8 
library(spdep)
9 9 
library(rgdal)
10 10 
library(spgwr)
11 

  		




  
  11
  
    
library(gpclib)


  		




  
  12
  
    
library(PBSmapping)


  		




  
  13
  
    
library(maptools)


  		




  
  14
  
    
library(gstat)


  		




  12
  15
  
    
###Parameters and arguments


  		




  13
  16
  
    

  		




  14
  17
  
    
infile1<-"ghcn_or_tmax_b_03032012_OR83M.shp" 


  		




  ......




  17
  20
  
    
#infile2<-"dates_interpolation_03012012.txt"  # list of 10 dates for the regression


  		




  18
  21
  
    
infile2<-"dates_interpolation_03052012.txt"


  		




  19
  22
  
    
prop<-0.3


  		




  20
  
  
    
out_prefix<-"_03132012_0"


  		




  
  23
  
    
out_prefix<-"_03272012_Res_fit"


  		




  
  24
  
    

  		




  
  25
  
    
###Reading the shapefile and raster image from the local directory


  		




  21
  26
  
    

  		




  22
  
  
    
###Reading the shapefile from the local directory


  		




  
  27
  
    
mean_LST<- readGDAL("mean_day244_rescaled.rst")  #This reads the whole raster in memory and provide a grid for kriging


  		




  23
  28
  
    
ghcn<-readOGR(".", "ghcn_or_tmax_b_03032012_OR83M") 


  		




  
  29
  
    
proj4string(ghcn) #This retrieves the coordinate system for the SDF


  		




  
  30
  
    
CRS_ghcn<-proj4string(ghcn) #this can be assigned to mean_LST!!!


  		




  
  31
  
    
proj4string(mean_LST)<-CRS_ghcn #Assigning coordinates information


  		




  
  32
  
    

  		




  
  33
  
    
# Creating state outline from county


  		




  
  34
  
    

  		




  
  35
  
    
orcnty<-readOGR(".", "orcnty24_OR83M")


  		




  
  36
  
    
proj4string(orcnty) #This retrieves the coordinate system for the SDF


  		




  
  37
  
    
lps <-getSpPPolygonsLabptSlots(orcnty)  #Getting centroids county labels


  		




  
  38
  
    
IDOneBin <- cut(lps[,1], range(lps[,1]), include.lowest=TRUE)  #Creating one bin var


  		




  
  39
  
    
gpclibPermit() #Set the gpclib to True to allow union


  		




  
  40
  
    
OR_state <- unionSpatialPolygons(orcnty ,IDOneBin) #Dissolve based on bin var


  		




  
  41
  
    

  		




  
  42
  
    
# Adding variables for the regression


  		




  24
  43
  
    

  		




  25
  44
  
    
ghcn$Northness<- cos(ghcn$ASPECT) #Adding a variable to the dataframe


  		




  26
  45
  
    
ghcn$Eastness <- sin(ghcn$ASPECT)  #adding variable to the dataframe.


  		




  ......




  35
  54
  
    
results <- matrix(1,length(dates),3)            #This is a matrix containing the diagnostic measures from the GAM models.


  		




  36
  55
  
    

  		




  37
  56
  
    
#Screening for bad values


  		




  38
  
  
    
#tmax~ lon + lat + ELEV_SRTM + Eastness + Northness + DISTOC


  		




  39
  57
  
    
#tmax range: min max)


  		




  40
  58
  
    
ghcn_test<-subset(ghcn,ghcn$tmax>-150 & ghcn$tmax<400)


  		




  41
  59
  
    
ghcn_test2<-subset(ghcn_test,ghcn_test$ELEV_SRTM>0)


  		




  42
  60
  
    
ghcn<-ghcn_test2


  		




  43
  
  
    
#lon range


  		




  44
  
  
    
#lat range


  		




  45
  
  
    
#ELEV_SRTM


  		




  46
  
  
    
#Eastness


  		




  47
  
  
    
#Northness


  		




  48
  61
  
    

  		




  49
  62
  
    
#ghcn.subsets <-lapply(dates, function(d) subset(ghcn, date==as.numeric(d)))#this creates a list of 10 subsets data


  		




  50
  63
  
    
ghcn.subsets <-lapply(dates, function(d) subset(ghcn, ghcn$date==as.numeric(d)))


  		




  51
  64
  
    

  		




  52
  
  
    
#ghcn.subsets <-subset(ghcn,ghcn$date=="20100101")


  		




  53
  
  
    
#summary(lm(y~x,data=df,weights=(df$wght1)^(3/4))


  		




  54
  
  
    
#ggwr.sel: find the bandwith from the data


  		




  55
  
  
    

  		




  56
  
  
    

  		




  57
  
  
    

  		




  58
  65
  
    
###Regression part 1: Creating a validation dataset by creating training and testing datasets


  		




  59
  66
  
    
for(i in 1:length(dates)){            # start of the for loop #1


  		




  60
  67
  
    
  


  		




  ......




  86
  93
  
    
  results[i,1]<- dates[i]  #storing the interpolation dates in the first column


  		




  87
  94
  
    
  results[i,2]<- ns        #number of stations used in the training stage


  		




  88
  95
  
    
  results[i,3]<- RMSE_f


  		




  
  96
  
    
  


  		




  
  97
  
    
  #Kriging residuals!!


  		




  
  98
  
    
  X11()


  		




  
  99
  
    
  hscat(residuals~1,data_s,(0:9)*20000) # 9 lag classes with 20,000m width


  		




  
  100
  
    
  v<-variogram(residuals~1, data_s)


  		




  
  101
  
    
  plot(v)


  		




  
  102
  
    
  tryCatch(v.fit<-fit.variogram(v,vgm(1,"Sph", 150000,1)),error=function()next)


  		




  
  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 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))


  		




  
  110
  
    

  		




  
  111
  
    
  plot(OR_state, axes = TRUE, add=TRUE)


  		




  
  112
  
    
  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")


  		




  
  116
  
    

  		




  
  117
  
    
  savePlot(paste(data_name,out_prefix,".png", sep=""), type="png")


  		




  
  118
  
    
  dev.off()


  		




  89
  119
  
    
  }


  		




  90
  120
  
    
  


  		




  91
  121
  
    
## Plotting and saving diagnostic measures


  		




  ......




  99
  129
  
    

  		




  100
  130
  
    
write.csv(results_table, file= paste(path,"/","results_GWR_Assessment",out_prefix,".txt",sep=""))


  		




  101
  131
  
    

  		




  
  132
  
    

  		




  102
  133
  
    
# End of script##########


  		




  103
  134
  
    

  		




  104
  135
  
    
# ###############################


  		












Also available in:  Unified diff