4

I'm trying to extract population density data from a raster by shapefile. I get NaN values, but I think my problem is further upstream at the raster level. I don't do much GIS work, so I'm likely making a rookie mistake here.

The population data source is interesting. It comes from an effort by Facebook to create population density maps from satellite imagery.

What I do below is grab the tif, use the raster package to load the file, and extract and average the density values by constituency (an administrative unit) with a shapefile I load.

The first problem is that the raster plot is blank (but the legend appears).

enter image description here

The second problem probably flows from the first: I am not able to extract any data.

library(raster)
library(velox)
library(sf)

# download the raster data 
# https://data.humdata.org/dataset/highresolutionpopulationdensitymaps-ken
  download.file("https://data.humdata.org/dataset/2964b369-c10c-4b55-94a8-495de3fc9858/resource/1271b53b-30be-4c4b-aa27-1bd104506636/download/population_ken_2018-10-01.zip", 
                destfile = 'popden.zip')

# unzip the file
  unzip(zipfile = "popden.zip", 
        exdir = 'popden')

# load raster
  r <- raster("popden/population_ken_2018-10-01.tif")
  plot(r)

# download the constituency shapefile
# http://www.opendata.go.ke/datasets/constituencies-2009/data
  ken_con_sf <- 
    read_sf("https://raw.githubusercontent.com/mikelmaron/kenya-election-data/master/data/constituencies.geojson") %>%
    st_transform(4326)
  plot(ken_con_sf)

# extract
  r_v <- velox(r)
  ken_con_den_mean <- r_v$extract(ken_con_sf, fun=mean)

@Spacedman correctly pointed out that my initial raster plot is not blank. Here's my new attempt to set breaks and visualize.

# get breaks as quantiles
# very slow so commenting out and giving values
  #rdf <- rasterToPoints(r) # very slow
  #rdf <- data.frame(rdf)
  #colnames(rdf) <- c("X","Y","den")
  #rb <- as.numeric(quantile(rdf$den, prob = seq(0, 1, length = 11)))

  rb <- c(0.002665759, 0.003351455, 0.005186487, 0.017159642, 0.048804275, 2.909804751, 3.859476157, 4.869144348, 5.943431054     7.815839071, 13091.488689846)

# plot
  plot(r, breaks = rb,  col = rainbow(100))

enter image description here

It's a little better, but still seems off to me.

To @Spacedman's second suggestion, I updated r_v$extract to run r_v$extract(ken_con_sf, fun=function(x){mean(x,na.rm=TRUE)}). This appears to have done the trick.

3
  • 2
    The raster is not blank. I can see some gray on your plot. The pixels probably have very skewed values. If most of them are around 0 to 10, and only one pixel is 700, then everything except one pixel will be white or gray, and you won't see that one green pixel because its like 30,000x20,000 and your screen is probably only 3000x2000 max.
    – Spacedman
    Aug 15 '19 at 22:30
  • 1
    And the NaNs are probably due to some missing data pixels over your polygons. Try: r_v$extract(ken_con_sf, fun=function(x){mean(x,na.rm=TRUE)})
    – Spacedman
    Aug 16 '19 at 10:21
  • Whoa! Good eyes @Spacedman. I thought it was a dirty screen.
    – Eric Green
    Aug 16 '19 at 14:33
2

Here I propose a function that works in parallel and speeds up this extraction process. I just tried with some of your polygons but, as you can see in the last part of the code, in my pc it will take more than an hour to extract all the mean values for all the polygons (in a ryzen 7 3700x using all the cores...).

In resume, this is what the function does:

  1. For each polygon --> crop the raster
  2. Extract de values
  3. Mean of the values (getting rid of the NaN)
  4. Create a data.frame with two columns: a selected ID field (specified in the function) and the mean value for the polygon (r_mean).
  5. Join the created data.frame with the original polygons by the ID.

Here is the code (yours and the function):

# rm(list = ls()) #uncomment to clean all environment
library(raster)
library(velox)
library(sf)
library(mapview)
library(doParallel)
library(dplyr)

# download the raster data 
# https://data.humdata.org/dataset/highresolutionpopulationdensitymaps-ken
download.file("https://data.humdata.org/dataset/2964b369-c10c-4b55-94a8-495de3fc9858/resource/1271b53b-30be-4c4b-aa27-1bd104506636/download/population_ken_2018-10-01.zip", 
              destfile = 'popden.zip')

# unzip the file
unzip(zipfile = "popden.zip", exdir = 'popden')

# load raster
r <- raster("popden/population_ken_2018-10-01.tif")

# download the constituency shapefile
# http://www.opendata.go.ke/datasets/constituencies-2009/data
pols <- 
    read_sf("https://raw.githubusercontent.com/mikelmaron/kenya-election-data/master/data/constituencies.geojson") %>%
    st_transform(4326)

# view layers
# mapview(r)+pols

# view pols data to get the ID
pols

#-------------------------------------------------------------------------------
# FUNCTION
#-------------------------------------------------------------------------------
# Function to extract mean raster values in polygons (using parallel processing)
MeanRasterInPolygons <- function(r, p, id = 'OBJECTID', cores=16){
    # register parallel
    cl <- parallel::makeCluster(cores, type="FORK")
    doParallel::registerDoParallel(cl)
    
    # extract mean raster value within each polygon
    plist <- foreach(x=seq(1, nrow(p))) %dopar% {
        # precrop raster to reduce raster size
        rcrop <- crop(r, p[x,])
        ex <- raster::extract(rcrop, p[x,], method = "simple")
        # calculate mean with numeric values (not na)
        exmean <- mean(ex[[1]][ex[[1]]!='NaN'])
        # compound df with pol id and mean
        return(Map(data.frame, ID=p[x,][[id]], r_mean=exmean)[[1]])
    }
    
    # stop parallel
    stopCluster(cl)
    
    # merge results after parallel process
    pi <- do.call(rbind, plist)
    
    # sort by id
    pi <- pi %>% arrange(ID)
    
    # rename df to match ID later on
    names(pi)=c(id, "r_mean")
    
    return(pi)
}
#-------------------------------------------------------------------------------
# PROCESS
#-------------------------------------------------------------------------------
# optional --> subset randomly your pols to check if this works and estimate how it will take for the whole process
pols.sel <- pols[sample(nrow(pols),20),]
# extract values and print processing time...
system.time({df <- MeanRasterInPolygons(r, pols.sel, id ='OBJECTID', cores = 16)})
# Stimated time in AMD ryzen 7 3700x and 64 RAM
# time for 20 polygons: 330 s --> +- 5 min
# stimated time for 295 features in pols --> +- 74 min

# join mean raster value to polygons and show in map
pols.join <- left_join(pols,df,by='OBJECTID')

enter image description here

When you are ready to process it all, change pols.sel by pols and remove that subseting at the end.

Here is what it looks (just showing a random sample of results for 20 polygons):

enter image description here

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.