-
Notifications
You must be signed in to change notification settings - Fork 0
/
GWR_2019_fok.R
55 lines (46 loc) · 1.36 KB
/
GWR_2019_fok.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
#### Geographically Weighted Regression (GWR)
### Author: Ying-Jung Deweese
### Date:2019-03-24
## Load library
library(spgwr)
library(sf)
library(sp)
library(spData)
library(Rcpp)
library(spdep)
library(ggplot2)
library(rgdal)
library(lattice)
library(maptools)
library(raster)
library(robustbase)
library(GWmodel)
library(RColorBrewer)
### Read the rainfall information
P=read.csv("allok_WY2005_fok.csv",header=T)
### Read Digital Elevation Model(DEM)
DEM <- raster("all_ok.tif")
##check the file dimension
S=dim(P)[2]-7
year=2005
## generate the output matrix
out=matrix(NA,S,1,byrow=T)
## generate spatialpoint dataframe
locations=cbind(P$lon,P$lat)
P.spdf=SpatialPointsDataFrame(P[,3:4],P)
##Compute basic GWR model and export the output
for (i in 1:length(S)){
S_index = i + 7
best.bw =gwr.sel(P[,S_index]~elv,data=P,coords=locations)
gwr.res=gwr.basic(P[,S_index]~elv,data=P.spdf,bw=best.bw,kernel='gaussian')
### estimated coeffient
elv.coefs <- gwr.res$SDF$elv
int.coefs=gwr.res$SDF$Intercept
ME=median(elv.coefs,)
MEIN=median(int.coefs,)
### Apply 10m DEM to the regression model results
rain_day=DEM*ME+MEIN
ras2 <- as(rain_day, "SpatialGridDataFrame")
ras2_name = sprintf("All_%d_S%d_GWR.tiff", year, i)
writeGDAL(ras2,ras2_name, drivername="GTiff", options=NULL)
}