1

I am attempting to (a) combine several non-overlapping pop density raster files, (b) create a velox object from this raster, (c) and then extract the mean population value for each overlapping polygon from a shapefile I provide that covers the big raster.

I'm stuck on (a) with the error (R):

Warning 6: gdalbuildvrt does not support heterogeneous band characteristics.

I don't know enough about raster files to know whether I need to somehow standardize the input rasters or pass an additional parameter to one of these functions.

library(gdalUtils)
library(rgdal)

# rasters
# 42: https://data.humdata.org/dataset/southasia_as42-high-resolution-population-density-maps
# 43: https://data.humdata.org/dataset/southasia_as43-high-resolution-population-density-maps
# 44: https://data.humdata.org/dataset/southasia_as434-high-resolution-population-density-maps
# 47: https://data.humdata.org/dataset/southasia_as47-high-resolution-population-density-maps

# zip with all, ~600MB: https://www.dropbox.com/s/26tc77upgpr51sr/rasterExample.zip?dl=0

# example https://stackoverflow.com/a/39554006/841405
  all_my_rasts <- list.files(pattern="\\.tif$")
  e <- extent(61.87125, 97.73931, 5.91125, 37.27042)
  template <- raster(e)
  projection(template) <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
  writeRaster(template, file="MyBigNastyRasty.tif", 
              format="GTiff",
              overwrite=TRUE)
  mosaic_rasters(gdalfile=all_my_rasts,
                 dst_dataset="MyBigNastyRasty.tif",
                 of="GTiff")
  gdalinfo("MyBigNastyRasty.tif")

Update 1:

I verified that all rasters have type "FLT8S". Removing rest of old update.

Update 2:

Building on @Jeffrey_Evans answer below, here's my first attempt at breaking up the processing. Basically, extract 2 rasters and rbind, then rbind that result to the third raster, and then to the forth, in the process removing things we don't need.

Here's where I stand at the end of step 1.

enter image description here

enter image description here

Is there a better way to approach this chunking?

# step 1
  e1 <- exact_extract(raster(r[1]), p) 
  e2 <- exact_extract(raster(r[2]), p) 
  results12 <- list()
  for(i in 1:nrow(p)) {
    results12[[i]] <- na.omit(rbind(e1[[i]], e2[[i]]))
  }
  remove(e1)
  remove(e2)
  gc()

# step 2
  e3 <- exact_extract(raster(r[3]), p) 
  results123 <- list()
  for(i in 1:nrow(p)) {
    results123[[i]] <- na.omit(rbind(results12[[i]], e3[[i]]))
  }
  remove(e3)
  remove(results12)
  gc()

# step 3
  e4 <- exact_extract(raster(r[4]), p) 
  results1234 <- list()
  for(i in 1:nrow(p)) {
    results1234[[i]] <- na.omit(rbind(results123[[i]], e4[[i]]))
  }
  remove(e4)
  remove(results123)
  gc()

  lapply(results1234, FUN=function(x) mean(x[,1]))
  p$rmean <- unlist(lapply(results1234, FUN=function(x) mean(x[,1])))
7
  • Your example makes no sense, the templet raster has no values associated with it. Can you please clarify, do your rasters share a consistent overlapping extent or do you want to combine non-overlapping rasters into a new single raster covering the cumulative extents? Commented May 13, 2020 at 14:47
  • Maybe I've misunderstood the SO answer I'm trying to follow (link provided). I understood the process to start with creating a blank template that has the full extent of the non-overlapping rasters.
    – Eric Green
    Commented May 13, 2020 at 14:59
  • Have you tried raster::merge or raster::mosaic? Commented May 13, 2020 at 15:22
  • I tried raster::merge but it never finished overnight. I'll look into raster::mosaic. Commenters elsewhere on SO have noted that the solution I'm trying to follow is a lot faster.
    – Eric Green
    Commented May 13, 2020 at 15:24
  • Try assigning values to your template raster before writing it to disk eg., template[] <- rep(1, ncell(template)) will write a raster with a single value. Also, check to see if the input rasters are the same bit type ie., float verses integer. If you create raster objects from each one you can check this or use gdalinfo. If they are different bit types this could account for the error. Commented May 13, 2020 at 15:34

1 Answer 1

1

So, being that there is always more than one way to skin a cat. if you really just need to get polygon statistics across all of the rasters, I can offer a very efficient alternative. Although, there may be other reasons that you need a raster mosaic.

You can extract values from each raster separately and then combine the resulting list objects, dropping NA values for times that there is no overlap in the polygon set. Here is an example on your data. Please note that I am using the exactextractr package for raster extraction, which is just as fast as velox. It provides two distinct advantages over velox. First, it does not require coercing to a velox object, saving overhead. Second, and most important, it is not a depreciated package (no longer on CRAN). This example clocked in at 5.23 seconds.

library(raster)
library(exactextractr)
library(sf)

setwd(...)
r <- list.files(getwd(), pattern="tif$")

Here we create two polygons where polygon 1 crosses the boundary of AS42 and AS43 and polygon 2 crosses the boundary of AS44 and AS47. In this way, for our example, we have multiple polygons that have partial values across multiple rasters. We coerce it to an sf object for the exact_extract function.

p <- rbind(SpatialPolygonsDataFrame(as(extent(72.7976,73.56307,
       29.31558,29.86756), "SpatialPolygons"), data.frame(ID=1)), 
       SpatialPolygonsDataFrame(as(extent(84.31119, 84.44135, 23.83541,
       23.94003 ), "SpatialPolygons"), data.frame(ID=2))) 
  proj4string(p) <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'  
p <- as(p, "sf")

Now, we extract values for each separate raster. This will result in list objects containing data.frames of raster values and cell fractions for each polygon.

e1 <- exact_extract(raster(r[1]), p) 
e2 <- exact_extract(raster(r[2]), p) 
e3 <- exact_extract(raster(r[3]), p) 
e4 <- exact_extract(raster(r[4]), p)

Now comes the part where we combine everything together. First, create and empty list to hold the results. Then we create a for loop that iterates through the polygon raster values and combines the data.frame from each list, dropping NA values. Now we have a list the would be the same as if you extracted values from a single mosaic raster. You can then use lapply to get summary statistics for the polygons. Just remember that the list contains data.frames so, you need to index the value column (see example).

results <- list()
  for(i in 1:nrow(p)) {
    results[[i]] <- na.omit(rbind(e1[[i]], e2[[i]], e3[[i]], e4[[i]]))
  }

lapply(results, FUN=function(x) mean(x[,1]))
  p$rmean <- unlist(lapply(results, FUN=function(x) mean(x[,1])))

Given the size of your data, if you have numerous polygons covering large areas you could run out of RAM in short order. In this case you could structure things in a way that you are iterating through polygons. This will do one polygon at a time but, you could combine this with the above solution to process multiple polygons at a time.

results <- vector()
  for(i in 1:nrow(p)) {
    p.sub <- p[i,]
    d <- na.omit(rbind(exact_extract(raster(r[1]), p.sub)[[1]],
            exact_extract(raster(r[2]), p.sub)[[1]], 
            exact_extract(raster(r[3]), p.sub)[[1]],
            exact_extract(raster(r[4]), p.sub)[[1]]))
    results[i] <- mean(d[,1])
  }
5
  • Incredibly helpful. Yes, this is what I am trying to do. Made it to the na.omit step. R crashed because I ran out of memory, but I can work on solving that.
    – Eric Green
    Commented May 13, 2020 at 19:17
  • 1
    For big polygon data you may have to loop through the polygons themselves. Commented May 13, 2020 at 19:54
  • I edited my post to show how I've tried to adapt your answer. Could you explain more about the looping you have in mind? I thought for(i in 1:nrow(p)) was looping already.
    – Eric Green
    Commented May 13, 2020 at 21:05
  • 1
    I will expand my answer in a bit. Yes there is already a loop but it is operating on values extract from all of the polygons in one cell swoop. For large data it is necessary to break up the problem. This was the next issue you were going to face anyway so, may as well deal with it now. Commented May 13, 2020 at 23:03
  • Thanks, @Jeffrey_Evans! That would be excellent. I really appreciate the help.
    – Eric Green
    Commented May 14, 2020 at 10:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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