# Geographically weighted area-to-area Regression CoKriging for raster downscaling in R?

this will be probable a series of questions, so in this one I will try to set the big picture. I want to perform `geographically weighted area-to-area regression Cokriging` (GWATARCoK) to downscale 1 satellite image using 2 covariates. The general workflow is:

1. upscale the covariates to match the resolution of the dependent variable
2. perform GWR and extract the residuals of the GWR
3. downscale the residuals using ATACoK
4. Add the downscaled residuals (from step 3) back to the GWR (not sure for this part yet). The code I have used so far (skipping the 1st bullet):
``````    ### part 2 ###
library(spgwr)
library(sf)

#create separate df for the x & y coords
x = as.data.frame(block.data\$x)
y = as.data.frame(block.data\$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 ~ ndvi + ndbi

# find optimal ADAPTIVE kernel bandwidth using cross validation
abw <- gwr.sel(eq1, data = block.data, adapt = T, gweight = gwr.Gauss);

# fit a gwr based on adaptive bandwidth
ab_gwr <- gwr(eq1,
data = block.data,
gweight = gwr.Gauss,
hatmatrix = T,
se.fit = T);
#print the results of the model
ab_gwr

#attach the coefficients to a dataframe
sp <- ab_gwr\$SDF
sf <- st_as_sf(sp)

#convert the residuals of the GWR to a raster file and export it
map.resids <- as.data.frame(sf\$gwr.e)
map.resids <- SpatialPointsDataFrame(data=map.resids, coords=cbind(x,y))
gridded(map.resids) <- TRUE
r <- raster(map.resids)
writeRaster(r, filename = "path/gwr_resids.tif", format = "GTiff")
### end of part 2 ###

### part 3 ###
library(atakrig)
library(raster)
library(beepr)

gwr_resids = raster("path/gwr_resids.tif") #dependent var
ndbi = raster("path/ndbi1.tif") # covariate_1
ndvi = raster("path/ndvi1.tif") #covariate_2

#discretization of raster
gwr_resids.d <- discretizeRaster(gwr_resids, 100);
ndbi.d <- discretizeRaster(ndbi, 100);
ndvi.d = discretizeRaster(ndvi, 100);
grid.pred <- discretizeRaster(ndvi, 100, type = "value"); #discretized grid to be predicted

# point-scale cross-variogram
aod.list <- list(ndvi = ndvi.d, ndbi = ndbi.d, gwr_resids = gwr_resids.d)
sv.ck <- deconvPointVgmForCoKriging(aod.list,
model = "Sph",
rd = 0.8); beep(7)

### area-to-area CoKriging
pred.atack <- ataCoKriging(aod.list,
unknownVarId = "gwr_resids",
unknown = grid.pred,
ptVgms = sv.ck,
oneCondition = T,
showProgress = T); beep(7)

#convert the result to raster format
pred.atack.r <- rasterFromXYZ(pred.atack[,-1])
writeRaster(pred.atack.r, 'path/atack.tif')
### end of part 3 ###
``````

The value range of the resulting raster is completely wrong (-2 millions to 9 millions). At this stage I'd like to get recommendations about the procedure I am following and if someone has any experience in this topic. In case someone wants the data, here. I cropped the study area to an extent that even `deconvPointVgmForCoKriging` function should not take more than 5 mins to execute (I am using a 10 y.o. laptop).

In case it helps, this is how the `point-scale cross-variogram(s)` before the `prediction` looks like Here are the downscaling raster files for `AtAK`, `AtAK_combine` and `AtACoK`.

• When performing `area-to-area kriging` (using 1 explanatory variable), the resulting downscaling raster is fine. In `AtACoK` I tried to log-transform the data (like the authors of the package `atakrig` did), I do not have extreme negative values now (actually the min value is 0), but I do get +infinity. Feb 17, 2022 at 17:21
• Also, I changed the NaN values (-3.4e+38) to 0 in case the algorithm was "accidentally" reading them as actual numbers, but no luck. Feb 21, 2022 at 19:32
• Moreover, I changed the data type of the data (16bit, Float 32 etc) but the results was the same, astronomical values for the CoKriging. Feb 22, 2022 at 12:32

Check the correlation of covariates with your target variable. According to Rossiter, 2018, (Technical Note: Co-kriging with the gstat package of the R environment for statistical computing): candidates with co-variable (target variable) must have:

1. a feature-space correlation with the target variable (fancy way to say pearson correlation and scatter-plot)
2. a spatial co-variance with the target variable. From your (cross)variogram one can understand the covariance between your target and co-variables is very small or even zero at small distances.

Also, check the histogram of the downscaled raster (atack.tif). Sometimes 2 pixels with infinite values (-Inf:+Inf) can stretch it and cause this issue.