Aggregate raster by a non-integer factor with arbitrary function
Asked Answered
D

1

6

I would like to aggregate a population raster by a factor of 1.5, summing the values of the cells.

While aggregate() allows me to sum values when aggregating, its factor parameter accepts only integer values. projectRaster() and resample() allow me to adjust the resolution precisely, but (as far as I know) I am restricted to the prepackaged bilinear-interpolation and nearest-neighbor computation methods.

Is there a way to aggregate a raster by a non-integer factor AND specify the function to use when aggregating?

library(raster)
set.seed(10)

proj <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
r <- raster(resolution = 1, nrow = 100, crs = proj)
r[] <- round(rnorm(ncell(r), 100, 10))

# Doesn't accept non-integer factors
aggregate(r, fact = 1.5, fun = sum)

template <- raster(extent(r), crs = crs(r), resolution = 1.5)

# Correct resolution, but incorrect / impossible values for population
projectRaster(r, to = template, method = "ngb")
projectRaster(r, to = template, method = "bilinear")

Possible workaround

So far, the only method I've been able to come up with is to coerce the template to a SpatialPoints object; extract values from the original, higher-resolution raster; and rasterize() the result:

pts <- as(template, "SpatialPoints")
vals <- extract(r, pts)
pts2 <- SpatialPointsDataFrame(pts, data.frame(vals))

rasterize(pts2, template, field = "vals", fun = sum)

However, if the points are created at the centroids of the raster cells, I'm not sure how they are handled when extracting with a resolution of 1.5x the original raster. My preferred method would be to create a SpatialPolygonsDataFrame and rasterize with fun = mean, but (in my experience) extracting raster values using polygons is very inefficient.

Diversity answered 12/3, 2017 at 14:8 Comment(2)
just for a better understanding. why do you want the sum, while resampling from 1 to 1.5 resolution ?Corps
In 1 to 1.5 for a population raster, maybe an average weighted by area would make more sense? I guess I just meant that, if we were to think of this raster as actually reporting the distribution of people over space, I would like an estimate of the "sum total" population living in the new (1.5x-area) cells.Diversity
C
1

Here a workaround:

#first resample to higher resolution
template    <- raster(extent(r), crs = crs(r), resolution = .5)  
detailedRas <- projectRaster(r, to = template, method = "ngb")

#then use an integer as a factor (in this case 3)
aggRas      <- aggregate(detailedRas, fact=3, fun=sum)

Note however that sum in this case won't return the sum of people who are living in a certain aggregated area.

I.e.: Let's say we have four cells and these values with a resolution of 1 m:

10  15  
12  18  

after resampling to 0.5 using NN:

10  10  15  15  
10  10  15  15  
12  12  18  18  
12  12  18  18  

Then aggregating by sum to 1.5m we get for the first pixel:

10+10+15+10+10+15+12+12+18 = 112

While in fact it should be something like: 10 + 15/2 + 12/2 + 18/4 = 28 (if we assume an equal population distribution over each pixel.)

I would recommend using the focal raster function with a custom / user defined function for summing up population values as you wish.

Or you divide the resampled raster by 4 and then take the sum:

2.5  2.5  3.75  3.75  
2.5  2.5  3.75  3.75   
3    3    4.5   4.5  
3    3    4.5   4.5  

2.5 + 2.5 + 3.75 + 2.5 + 2.5 + 3.75 + 3 + 3 + 4.5 = 28

Corps answered 14/3, 2017 at 19:40 Comment(3)
Thank you for your answer. Are there any circumstances in which bilinear interpolation would be better for creating the detailed raster? (I'd like to write a function for this, and I'm wondering if I should give that function a method parameter to use in the projectRaster() step.)Diversity
I would stick with nn "interpolation". Its designed for discrete variables like population. Bilinear is used for continuous variables and the interpolated values are not as easily reproducible as the nn interpolation.Corps
In case anyone would like to see this in the form of a function: gist.github.com/coletl/e67620b8e059cc85538346a5863bc19c. Contributions and feedback are more than welcome.Diversity

© 2022 - 2024 — McMap. All rights reserved.