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

Question

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)
    block.data = read.csv(file = "path/block.data.csv")
    
    #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,
                  adapt = abw, 
                  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, 
                               auxRatioAdj = 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.

Answer 1

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.