3

I'm creeping into building surfaces - useful for a number of public health projects - and can't seem to make sense of what's happening with the resolution and raster building process to get what I want. I've got a hopefully trivial and reproducible example below to build from. Of course, for projection reasons, this might not be sensible at this scale (what I'm really doing is more at the block/block group/tract level), but I think this makes for a digestible talking point.

I'd like to get the population within a big circle in the US. To that end, I've got a shapefile, some population data, and a ~1000 km circle smack dab in the heartland.

How can I figure the population in that circle?

Picture and broken code below.

Thoughts?

Map of shapefile over US raster with circle in the middle

  #################################################################
  # "I don't get rasters"
  # an example by Mike Dolan Fliss. NC Epidemiology.
  #################################################################

  #NOTE: You'll need to install.packages() these if you don't use them.
  library(rgdal)
  library(sp)
  library(dplyr)
  library(ggplot2) # Could get better plots with ggplot, but for speed we'll use base/sp
  library(rgeos)
  library(raster)
  library(xlsx) #For the description table, an excel file

  # Download, read, and touch up the geographic data ##############
  #################################################################
  download.file("http://www2.census.gov/geo/tiger/TIGER2010DP1/State_2010Census_DP1.zip", "states.zip") #~9meg
  unzip("states.zip")
  states.spdf = readOGR(dsn = ".", layer = "State_2010Census_DP1", stringsAsFactors = F)
  states.spdf = states.spdf[!(states.spdf$NAME10 %in% c("Alaska", "Hawaii", "Puerto Rico")), ] #get continental

  descriptions.df = read.xlsx("DP_TableDescriptions.xls", 1)
  head(descriptions.df) # Total pop is in : DP0010001
  states.spdf$total.pop = states.spdf$DP0010001 #For readability

  plot(states.spdf) #Plot in base
  spplot(states.spdf, "total.pop") #Plot in sp, gorgeous Lisa Frank colors
  spplot(states.spdf, "total.pop", col.regions=rev(heat.colors(20))) #heat

  # Project & build some buffers ##################################
  #################################################################
  epsg_codes = make_EPSG()# Let's go Albers equal area. I'm not great with projections.  
  us.equalarea.prj = epsg_codes$prj4[epsg_codes$note=="# US National Atlas Equal Area"]
  states.spdf = spTransform(states.spdf, CRS(us.equalarea.prj)) #unit is now m.

  center = gCentroid(states.spdf)
  little.buffer.spdf = gBuffer(center, width=1000*1000) #1000km radius "circle"
  big.buffer.spdf = gBuffer(center, width=2700*1000) #2700 radius - all of US
  plot(states.spdf);plot(big.buffer.spdf, add=T); plot(little.buffer.spdf, add=T)

  # Now, I fail to make the right raster... :) ####################
  #################################################################
  pop.raster = raster(extent(states.spdf))
  e = extent(states.spdf); (e@xmax-e@xmin)/1000 #US extent is 4560km across. ok...
  projection(pop.raster) = proj4string(states.spdf)
  raster.res = 1000*100 #makes 100km x 100km raster squares (?)
  res(pop.raster) = c(raster.res, raster.res) 
  states.spdf$frac.pop = states.spdf$total.pop / (states.spdf$ALAND10+states.spdf$AWATER10)*raster.res
  #head(bgs.spdf$frac.pop)
  pop.raster = rasterize(states.spdf, pop.raster, "frac.pop")

  plot(pop.raster) #Plot - looks good
  plot(states.spdf, border="black", add=T)

  cellStats(pop.raster, sum) #3093? Clearly I'm confused
  sum(states.spdf$total.pop) #US pop = 307 million

  # ... and subsequently fail to get the right area ###############
  #################################################################
  plot(rp)
  plot(states.spdf, border="black", add=T)
  plot(little.buffer.spdf, border="blue", add=T)
  # What's the area in this blue circle?

  e = extract(pop.raster, big.buffer.spdf)
  sum(unlist(e), na.rm=T) #3093... somethings
  #^ Here we are, back to the same as cellStats sum, above.

  e = extract(pop.raster, little.buffer.spdf)
  sum(unlist(e)) #663... 
  #^ I'd like this to represent the population contained in the circle
2

There are two problems with your code.

The first and most crucial is shown below:

states.spdf$frac.pop = states.spdf$total.pop / (states.spdf$ALAND10+states.spdf$AWATER10)*raster.res

You are assuming that population is equally distributed in space, however you treat density wrong. Instead of using raster.res to distribute population raster.res^2 as you are dealing with area. Correct code line is given next:

states.spdf$frac.pop = states.spdf$total.pop / (states.spdf$ALAND10+states.spdf$AWATER10)*raster.res^2

Implementing this line would get you closer to the real number:

sum(states.spdf$total.pop) #US pop = 307 million

[1] 306675006

cellStats(pop.raster, sum)

[1] 309351269

The second problem relates to the resolution in its effect on boundaries:

A 100 KM cell will probably overlay two states. Reduced cell size also reduces the occurrence of such overlays. Thus it improves population allocation to cells along boundaries. See results for 10 KM below:

cellStats(pop.raster, sum) # 10 KM cell size

[1] 305963050

sum(states.spdf$total.pop) #US pop = 307 million

[1] 306675006

In addition it has a significant effect on visual results. 10 KM Density map

Finally, note that there is a trade-off between resolution and run-time, so resolution should be chosen carefully.

  • This is a great start, thank you for taking a shot at this. Could you explain a little more about the res function? I wondered about whether I needed a square term, but isn't the res function the area (given I passed it an x and y length for the rectangle - in this case a square)? Second question: I appreciate that more resolution provides better... resolution! :) To be clear: the lack of accuracy to the total pop is due to the boundaries between regions (where populations are mixed within the raster unit). Higher resolution provides benefits there (and only there, really), yes? – Mike Dolan Fliss Jul 10 '16 at 17:46
  • res function is used to get / set the x and / or y resolution of a raster object. You can get more information in the raster package documentation. You used it well. Your object raster.res stand for either x or y length, thus when multiplying to distribute population across cells you need to use its square. Regarding your second question; as far as I can see it - yes. The lack of accuracy is due to boundaries and higher resolution will help you only there. – dof1985 Jul 10 '16 at 19:03
  • 1
    Hmm. Thanks for this. I did (of course!) read the package documentation on raster before posting, but I'm still confused on why I'm squaring when the ... oh shoot, nevermind. I wouldn't be squaring if I genuinely used the true cell area, res(pop.raster), but I used its (in this case, square) components. Which is why, duh, I'm squaring. More robust would be to use re(pop.raster) given I might use units that aren't square in the future. I think that's the click I needed... thanks a ton. Follow-up coming... – Mike Dolan Fliss Jul 10 '16 at 20:07
  • Can I ask a small follow-up? How does raster decide what the resolution will be (x and y) when I just set the unit area with res(r) directly? Last I did that the x and y were different. I couldn't find that logic in the package helpfiles, but must have missed it. That is, if I set res(pop.raster) = x directly, sometimes it's square, sometimes its not. – Mike Dolan Fliss Jul 10 '16 at 20:10
  • @MikeDolanFliss, I'm not sure about this question. It is better however if you can create another reproducible example and post it as another one; Additionally I'm not quite sure that this forum fits best such question, maybe stackoverflow. – dof1985 Jul 11 '16 at 6:27

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.