Apply geographically weighted regression's model parameters to a finer spatial scale

Question

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.

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

Edit

Because the problem has been solved the csv is available upon request

You can download the csv from here. The raster I want to apply GWR's coefficients with:

tirs = raster(ncols=407, nrows=342, xmn=509600, xmx=550300, ymn=161800, ymx=196000, crs='+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs')

please provide a *minimal* self-contained example. It should be possible to do set up your question with much less, code and without relying on external files. And do not include screenshots, copy the console output instead (and show what generated that). — Robert Hijmans, Dec 04 '22 at 06:32
I removed the screenshot and parts of the code for linear regression as well as I posted my raster data. Unfortunately, I couldn't replicate your example, using my data, for linear regression because your data have the same spatial resolution before the regression while my data have different pixel size. But I believe the principle is the same when I am using my data. — Nikos, Dec 04 '22 at 12:25

score 0 · Answer 1 · answered Dec 04 '22 at 06:52

This is how you can do global regression and predict to a higher resolution (downscale)

library(terra)
r <- rast(system.file("ex/logo.tif", package="terra"))
a <- aggregate(r, 10, mean)

model <- lm(formula = red ~ green, data=a)
p <- predict(r, model)

And

d <- as.data.frame(a[[1:2]], xy=TRUE)

Perhaps this helps to write a better example in your question.

score 0 · Accepted Answer · answered Dec 05 '22 at 18:13

The solution was to use the regression.point argument in the gwr.basic function. The steps were:

load the high resolution raster and convert it so SpatilPointsDataFrame (SPDF)
run GWR and apply the model to the SPDF

The code:

library(GWmodel)
library(sp)

tirs = raster("path/tirs.tif") # high resolution raster
regpoints <- as(tirs, "SpatialPoints")

block.data = read.csv(file = "path/block.data.psf.csv")

coordinates(block.data) <- c("x", "y")
proj4string(block.data) <- "EPSG:27700"

eq1 <- ntl ~ tirs000 # tirs000 is the coarse version of the high res raster
abw = bw.gwr(eq1, 
             data = block.data, 
             approach = "AIC", 
             kernel = "gaussian",
             adaptive = TRUE, 
             p = 2,
             parallel.method = "omp",
             parallel.arg = "omp")

ab_gwr = gwr.basic(eq1, 
                   data = block.data, 
                   regression.points = regpoints,
                   bw = abw, 
                   kernel = "gaussian", 
                   adaptive = TRUE, 
                   p = 2, 
                   F123.test = FALSE,
                   cv = FALSE,
                   parallel.method = "omp",
                   parallel.arg = "omp")

ab_gwr

sp <- ab_gwr$SDF
sf <- st_as_sf(sp)

# intercept
intercept = as.data.frame(sf$Intercept)
intercept = SpatialPointsDataFrame(data = intercept, coords = regpoints)
gridded(intercept) <- TRUE
intercept <- raster(intercept)
raster::crs(intercept) <- "EPSG:27700"

# slope
slope = as.data.frame(sf$tirs000)
slope = SpatialPointsDataFrame(data = slope, coords = regpoints)
gridded(slope) <- TRUE
slope <- raster(slope)
raster::crs(slope) <- "EPSG:27700"

gwr_pred000 = intercept + slope * tirs

writeRaster(gwr_pred000, 
            "path/gwr_pred000.tif", 
            overwrite = TRUE)

Apply geographically weighted regression's model parameters to a finer spatial scale

2 Answers2