R – Extracting Raster Cell Values from Gridded World Population Dataset

cell-statisticsextractrrastervector

I am using the Gridded World Population dataset which gives population density values for each cell in a raster. I want to find the population of shapefiles by extracting the population density values from this raster. Data loaded here:

library(sf)
library(dplyr)
library(raster)
library(geodata)
tdir=tempdir()
# Gather Lower 48 / DC USA sf shapefile
stateurl = "https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_state_500k.zip"
non_L48 = c("72","78","15","66","69","60","02")
if(file.exists(paste(tdir,"/cb_2018_us_state_500k.shp",sep=""))==F){
  download.file(stateurl, destfile = file.path(tdir, "States.zip"))
  unzip(file.path(tdir,"States.zip"),exdir=tdir)}
USA = read_sf(paste(tdir,"/cb_2018_us_state_500k.shp",sep="")) %>%
  filter(!c(STATEFP%in%non_L48))

# Convert sf back to sp:
USAsp = as(USA,Class="Spatial")
#Download GWP dataset
geodata::population(year=2015,res=.5,path=tdir)
#load via raster:
popdat = raster(x = paste(tdir,"/pop/gpw_v4_population_density_rev11_2015_30s.tif",sep=""))

#Crop to USA borders
USA_transform = st_transform(USA,4326)
cr_USA = crop(popdat,extent(USA_transform))
USA_popdat = mask(cr_USA,USA_transform)

# Rename:
names(USA_popdat)[1]="pop_dens"

I now have a raster dataset cropped to the Lower 48 US States and DC called USA_popdat. As mentioned, I would like to find the total population within each shapefile by extracting via raster::extract(). So far I have this:

extracted = raster::extract(USA_popdat,
                            USA_transform,
                            fun=sum,
                            na.rm=TRUE,
                            p = T)

However, this only returns the sum of population density values within each cell contained by each respective shapefile. The sum of population density using the above code is 490,010,296, which is much larger than the ~350,000,000 in the U.S. population.

Intuitively, to achieve a total population count rather than a sum of population density, we should sum the population density and divide it by the number of cells. I am unsure about how to enter this operation into the fun= argument within raster::extract()

Further, when dividing the sum of population density (490,010,296) by the total number of cells in the U.S. we achieve 23.67:

sum(extracted)/ncells(USA_popdat)

which leads me to believe I am doing something incorrectly with the gridded world population dataset: https://sedac.ciesin.columbia.edu/data/set/gpw-v4-population-density-rev11

Best Answer

A regular grid in lat-long has cells of varying area. So adding a population density while ignoring the cell area will lead to incorrect values of population.

Luckily there's a handy raster::area function that computes the area of each cell in a lat-long raster.

> USA_area = raster::area(USA_popdat)
> plot(USA_area)

This is varying from 0.75ish to 0.55ish (km^2), which is a substantial range and needs to be considered.

If we multiply the area in km^2 by the density in people/km^2 we get a map of population count per cell:

> USA_count = USA_popdat * USA_area
> plot(USA_count)

And if we sum that up we get:

> sum(USA_count[], na.rm=TRUE)
[1] 330292666

...a population of 330 million which seems believable.

You can then do your extraction (I recommend using the exactextractr package which is fast and considers where polygons clip edge cells) on the population count and if you want density then divide the count by the polygon areas.

exactextractr is really fast. I got fed up with waiting for raster::extract to do its stuff.

> extracts = exact_extract(USA_count, USA_transform, fun="sum")
  |======================================================================| 100%
> head(extracts)
[1]  3079368 10562649  3963717  8755188  1934866  4732980
> USA_transform$pop = extracts
> plot(USA_transform[,"pop"])

Related Solutions

[GIS] Extract Raster from Raster using Polygon shapefile in R

The extract function is behaving exactly as it should. You can force the crop function to use the extent of the polygon and then mask the object to return the exact raster representing the polygon area. If you continue to receive the error it means that your data, in fact, does not overlap.

Please keep in mind that R does not perform "on the fly" projection so, check your projections. You can check if your extents overlap using the "extent()" function.

Here is an example of cropping using the polygon extent then masking the resulting raster using the "rasterized" polygon.

# Add required packages
require(raster)
require(rgeos)
require(sp)

# Create some data using meuse 
data(meuse)
  coordinates(meuse) <- ~x+y
    proj4string(meuse) <- CRS("+init=epsg:28992")
data(meuse.grid)
  coordinates(meuse.grid) = ~x+y
    proj4string(meuse.grid) <- CRS("+init=epsg:28992")
      gridded(meuse.grid) = TRUE    
        r <- raster(meuse.grid) 
          r[] <- runif(ncell(r))

# Create a polygon
f <- gBuffer(meuse[10,], byid=FALSE, id=NULL, width=250, 
                         joinStyle="ROUND", quadsegs=10)   

# Plot full raster and polygon                       
plot(r)
  plot(f,add=T)

# Crop using extent, rasterize polygon and finally, create poly-raster
#          **** This is the code that you are after ****  
cr <- crop(r, extent(f), snap="out")                    
  fr <- rasterize(f, cr)   
    lr <- mask(x=cr, mask=fr)

# Plot results
plot(lr)
  plot(f,add=T)

R Raster Processing – Increasing Speed of Crop, Mask, & Extract by Polygons

I have finally gotten around to improving this function. I found that for my purposes, it was fastest to rasterize() the polygon first and use getValues() instead of extract(). The rasterizing isn't much faster than the original code for tabulating raster values in small polygons, but it shines when it came to large polygon areas that had large rasters to be cropped and the values extracted. I also found getValues() was much faster than the extract() function.

I also figured out the multi-core processing using foreach().

I hope this is useful for other people who want an R solution for extracting raster values from a polygon overlay. This is similar to the "Tabulate Intersection" of ArcGIS, which did not work well for me, returning empty outputs after hours of processing, like this user.

#initiate multicore cluster and load packages
library(foreach)
library(doParallel)
library(tcltk)
library(sp)
library(raster)

cores<- 7
cl <- makeCluster(cores, output="") #output should make it spit errors
registerDoParallel(cl)

Here's the function:

multicore.tabulate.intersect<- function(cores, polygonlist, rasterlayer){ 
  foreach(i=1:cores, .packages= c("raster","tcltk","foreach"), .combine = rbind) %dopar% {

    mypb <- tkProgressBar(title = "R progress bar", label = "", min = 0, max = length(polygonlist[[i]]), initial = 0, width = 300) 

    foreach(j = 1:length(polygonlist[[i]]), .combine = rbind) %do% {
      final<-data.frame()
      tryCatch({ #not sure if this is necessary now that I'm using foreach, but it is useful for loops.

        single <- polygonlist[[i]][j,] #pull out individual polygon to be tabulated

        dir.create (file.path("c:/rtemp",i,j,single@data$OWNER), showWarnings = FALSE) #creates unique filepath for temp directory
        rasterOptions(tmpdir=file.path("c:/rtemp",i,j, single@data$OWNER))  #sets temp directory - this is important b/c it can fill up a hard drive if you're doing a lot of polygons

        clip1 <- crop(rasterlayer, extent(single)) #crop to extent of polygon
        clip2 <- rasterize(single, clip1, mask=TRUE) #crops to polygon edge & converts to raster
        ext <- getValues(clip2) #much faster than extract
        tab<-table(ext) #tabulates the values of the raster in the polygon

        mat<- as.data.frame(tab)
        final<-cbind(single@data$OWNER,mat) #combines it with the name of the polygon
        unlink(file.path("c:/rtemp",i,j,single@data$OWNER), recursive = TRUE,force = TRUE) #delete temporary files
        setTkProgressBar(mypb, j, title = "number complete", label = j)

      }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")}) #trycatch error so it doesn't kill the loop

      return(final)
    }  
    #close(mypb) #not sure why but closing the pb while operating causes it to return an empty final dataset... dunno why. 
  }
}

So to use it, adjust the single@data$OWNER to fit with the column name of your identifying polygon (guess that could have been built into the function...) and put in:

myoutput <- multicore.tabulate.intersect(cores, polygonlist, rasterlayer)

Best Answer

Related Solutions

[GIS] Extract Raster from Raster using Polygon shapefile in R

R Raster Processing – Increasing Speed of Crop, Mask, & Extract by Polygons

Related Question