rrastercoefficientsgwr

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

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

Solution

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

    1. load the high resolution raster and convert it so SpatilPointsDataFrame (SPDF)
    2. 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)