[R-sig-Geo] Apply geographically weighted regression's model parameters to a finer spatial scale

[Roger Bivand]

On Sun, 4 Dec 2022, Nikolaos Tziokas wrote:

> I have two raster layers, one coarse resolution and one fine resolution. My
> goal is to extract GWR's coefficients (intercept and slope) and apply them
> to my fine resolution raster.

You are making life too complicated, just use the regression.points = 
argument to gwr.basic():

library(GWmodel)
library(raster)
tirs0 = terra::rast(ncols=407, nrows=342, nlyrs=1, xmin=509600,
  xmax=550300, ymin=161800, ymax=196000, names=c('tirs0'),
  crs='EPSG:27700')
tirs1 <- raster(tirs0)
regpoints <- as(tirs1, "SpatialPoints")
summary(regpoints)
block.data = read.csv(file = "block.data00.csv")
coordinates(block.data) <- c("x", "y")
proj4string(block.data) <- "EPSG:27700"
summary(block.data)
eq1 <- ntl ~ tirs
abw = bw.gwr(eq1, data = block.data, approach = "AIC", kernel = "tricube",
  adaptive = TRUE, p = 2)
ab_gwr = gwr.basic(eq1, data = block.data, regression.points = regpoints,
  bw = abw, kernel = "tricube", adaptive = TRUE, p = 2, F123.test = FALSE,
  cv = FALSE)
ab_gwr
summary(ab_gwr$SDF)

The final multiplication can be done with your actual tirs raster, which 
was not provided.

With the small number of input points, a geostatistical approach to 
prediction would be preferable to GWR, giving well-founded prediction 
standard errors.

Hope this clarifies,

Roger

>
> I can do this easily when I perform simple linear regression. For example:
>
> library(terra)
> library(sp)
> # focal terra
> tirs = rast("path/tirs.tif") # fine res raster
> ntl = rast("path/ntl.tif") # coarse res raster
>    # fill null values
> tirs = focal(tirs,
>             w = 9,
>             fun = mean,
>             na.policy = "only",
>             na.rm = TRUE)
>
> gf <- focalMat(tirs, 0.10*400, "Gauss", 11)
> r_gf <- focal(tirs, w = gf, na.rm = TRUE)
>
> r_gf = resample(r_gf, ntl, method = "bilinear")
>
> s = c(ntl, r_gf)names(s) = c('ntl', 'r_gf')
>
> model <- lm(formula = ntl ~ tirs, data = s)
> # apply the lm coefficients to the fine res raster
> lm_pred = model$coefficients[1] + model$coefficients[2] * tirs
>
> But when I run GWR, the slope and intercept are not just two numbers (like
> in linear model) but it's a range. For example, below are the results of
> the GWR:
> *Summary of GWR coefficient estimates*:
>
>                Min.     1st Qu.      Median     3rd Qu.     Max.
>
> Intercept -1632.61196   -55.79680   -15.99683    15.01596 1133.299
>
> tirs20      -42.43020     0.43446     1.80026     3.75802   70.987
>
>
> My question is how can extract GWR model parameters (intercept and slope)
> and apply them to my fine resolution raster? In the end I would like to do
> the same thing as I did with the linear model, that is, *GWR_intercept +
> GWR_slope * fine resolution raster*.
>
> Here is the code of GWR:
>
> library(GWmodel)
> library(raster)
>
> block.data = read.csv(file = "path/block.data00.csv")
> #create mararate df for the x & y coords
> x = as.data.frame(block.data$x)
> y = as.data.frame(block.data$y)
> sint = as.matrix(cbind(x, y))
> #convert the data to spatialPointsdf and then to spatialPixelsdf
> coordinates(block.data) = c("x", "y")#gridded(block.data) <- TRUE
> # specify a model equation
> eq1 <- ntl ~ tirs
>
> dist = GWmodel::gw.dist(dp.locat = sint, focus = 0, longlat = FALSE)
>
> abw = bw.gwr(eq1,
>       data = block.data,
>       approach = "AIC",
>       kernel = "tricube",
>       adaptive = TRUE,
>       p = 2,
>       longlat = F,
>       dMat = dist,
>       parallel.method = "omp",
>       parallel.arg = "omp")
>
> ab_gwr = gwr.basic(eq1,
>          data = block.data,
>          bw = abw,
>          kernel = "tricube",
>          adaptive = TRUE,
>          p = 2,
>          longlat = FALSE,
>          dMat = dist,
>          F123.test = FALSE,
>          cv = FALSE,
>          parallel.method = "omp",
>          parallel.arg = "omp")
>
> ab_gwr
>
> You can download the csv from here
> <https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdrive.google.com%2Fdrive%2Ffolders%2F1V115zpdU2-5fXssI6iWv_F6aNu4E5qA7%3Fusp%3Dsharing&data=05%7C01%7CRoger.Bivand%40nhh.no%7C6855c5767a894cdb4edf08dad609d991%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638057635036324587%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=ROLhw75M%2BvWhGOIT7mlyPfNZX2W0Ub4gMHgrTBVsd8A%3D&reserved=0>.
> The fine resolution raster I am using:
>
> tirs = rast(ncols=407, nrows=342, nlyrs=1, xmin=509600, xmax=550300,
> ymin=161800, ymax=196000, names=c('tirs'), crs='EPSG:27700')
>
>

-- 
Roger Bivand
Emeritus Professor
Department of Economics, Norwegian School of Economics,
Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
e-mail: Roger.Bivand using nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en