Search code examples
rrasterterra

disaggregate raster value to finer cells equally

I have a SpatRaster in the following dimension

  layer1

  class       : SpatRaster
  dimensions  : 3600, 3600, 1  (nrow, ncol, nlyr)
  resolution  : 0.0002777778, 0.0002777778  (x, y)
  extent      : -81.00014, -80.00014, 24.99986, 25.99986  (xmin, xmax, ymin, ymax)
  coord. ref. : lon/lat WGS 84 (EPSG:4326)
  source      : nn.tif
  name        : nn_1

I have another SpatRaster which has data for revenue (in dollar) at a different resolution than layer1

  layer2
  class       : SpatRaster
  dimensions  : 21600, 43200, 1  (nrow, ncol, nlyr)
  resolution  : 0.008333333, 0.008333333  (x, y)
  extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
  coord. ref. : lon/lat WGS 84
  source      : XX.nc
  varname     : XX
  name        : XX
  unit        : US dollar

What I want to do is to take the average of layer1 using layer2 as weights. This is my approach:

  1. Give them the same projections.

    layer2 <- project(layer2, crs(layer1))

  2. Clip layer2 to match the extent of layer1

     r <- terra::crop(layer2, ext(layer1))
 r <- terra::mask(r, layer1) # I get error in this step
 [mask] number of rows and/or columns do not match

  3. Suppose step 2 worked. Now I want to disaggregate clipped layer2 to match the resolution of layer1 i.e. 30 times finer. This means revenue of each disaggregated cell should be 1/30 of the revenue of the parent cell.
    Is the below the right way to do this:

     r_disagg <- disaggregate(r, fact = 30)/30

  4. Divide each cell value in r_disagg with total sum of r_disagg to get a weight raster r_disagg_wt

  5. Multiply r_disagg_wt with layer1 and add the resulting raster to get the weighted average of layer1

I am stuck in step 2 and 3

Solution

  • When asking an R question, please include example data like this

    Example data

    library(terra)
r1 <- rast(nrow=360, ncol=360, ext=c(-81.00014, -80.00014, 24.99986, 25.99986))
r2 <- rast(nrow=2160, ncol=4320)
values(r1) <- 1:ncell(r1)
set.seed(0)
values(r2) <- runif(ncell(r2))

    Solution

    re2 <- resample(r2, r1)
global(r1, fun="mean", weights=re2)
#      weighted_mean
#lyr.1      66190.15

    In cases where the two input rasters align, it could be preferable to use disagg instead of resample (in this case, the extent of r1 is a bit odd, it seems to be shifted with half a cell).

    There is no need to adjust the dollar values if they are only used as weights. In cases that the sum of the cell values should remain constant, when using resample or project, you could compute the density ($/km2) before resampling by dividing the values with the area of each cell (see cellSize) and multiply the values again with the area of the new cells after resampling.