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 ###
    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
    #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 ###
    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 point-scale cross-variogram(s) Here are the downscaling raster files for AtAK, AtAK_combine and AtACoK. AtAK_prediction AtAK_combine_prediction AtACoK_prediction

  • 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.
    – Dolby
    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.
    – Dolby
    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.
    – Dolby
    Feb 22, 2022 at 12:32

1 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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.