I have a specific question that leads to a more general one.

First, the specific: I have an analysis involving census blocks and point level exposures (see here), where I'm doing work with ~200,000 polygons. The specific question I'm facing involves looking for block-level neighbors with certain qualities, then I'd add a variable to associate that with the block for the purposes of expanding a regression we're working on. In short: I'm finding neighbors, then examining the neighbor set for qualities (for instance, let's say % Black/AA above, say, 40%), then associating those values with the blocks. In this test case, I'm looking for neighbors of "big" polygons.

First let's look at doing that with NC counties.

download.file("ftp://ftp2.census.gov/geo/pvs/tiger2010st/37_North_Carolina/37/tl_2010_37_county10.zip", destfile="NC.counties.zip")
#Read spatial poly data frame. sidenote: odd I *have* to specify dsn here.    
NCpolys.spdf = readOGR(dsn=".", "tl_2010_37_county10", stringsAsFactors = F) 

NCpolys.spdf$landarea.sqmi = NCpolys.spdf$ALAND10*3.8610216E-7 
#ALAND10 is Land area in square meters
NCpolys.spdf$big = as.integer(NCpolys.spdf$landarea.sqmi >  quantile(NCpolys.spdf$landarea.sqmi, probs = c(.75)))
#NCpolys.spdf@data = within(NCpolys.spdf@data, big = landarea.sqmi > quantile(landarea.sqmi, probs = c(.75))) \
#Sidenote I'd like to do above, but within doesn't seem to like spatial objects, and I don't trust @data to reorder it before assigning....

#touch.m = gTouches(NCpolys.spdf, byid = T, returnDense = T) #Creates nxn touch matrix.  Let's not do that, since that might not scale...
touch.l = gTouches(NCpolys.spdf, byid = T, returnDense = F) #<1m for 100 polygons

for (i in 1:length(touch.l)){
  NCpolys.spdf$bigneighbor[i] = max(NCpolys.spdf@data$big[touch.l[[i]]])
#List referencing with spatial and *apply confuses the heck out of me.  

Can anyone help with a sapply/equivalent here?

See below, I'm stumped. #NCpolys.spdf$bigneighbor = sapply(touch.l, function(x) max(NCpolys.spdf@data$big[[[x]]])) #Doesn't work on first post.

#And now, the results.
plot(NCpolys.spdf[NCpolys.spdf$bigneighbor==1,], co="light green", add=T)
invisible(text(coordinates(NCpolys.spdf[NCpolys.spdf$bigneighbor==1,]), labels="n", cex=.6))
plot(NCpolys.spdf[NCpolys.spdf$big==1,], co="light blue", add=T)
invisible(text(coordinates(NCpolys.spdf[NCpolys.spdf$big==1,]), labels="X", cex=.6))

#Now, the real question.  What if this were X00,000 polygons?
download.file("ftp://ftp2.census.gov/geo/pvs/tiger2010st/37_North_Carolina/37/tl_2010_37_tabblock10.zip", destfile="NC.blocks.zip") #202 megs
NCpolys.spdf = readOGR(dsn=".", "tl_2010_37_tabblock10", stringsAsFactors = F) #read blocks
#Try the above again.  gTouches, eek.

enter image description here

But I'm not working with 100 NC counties - I'm working with 2-300,000 census block polygons. In this specific example, even with "returnDense=F" (which seems backwards to me, as a sidenote) so I get a list form instead of a 200,000x200,000 true neighbor matrix, I'm still hitting memory problems. gTouches with ...

  • 10,000 blocks: <1m.
  • 50,000 blocks: 10m.
  • 100,000 blocks: >1h.

Seems like it's O(N^2), instead of a bounding box based method that I could imagine being closer to O(N).

I'm wondering if I could rewrite this as a for-loop and go block by block, avoiding the massive matrix, and if R could keep things a little cleaner and smaller that way... but anytime I find myself breaking down functions and avoiding vectorization (see my sapply avoidance in the code above).

Any specific improvements to the above? And any general guidance?

For instance, I've heard rumors about more and less optimized packages, enabling certain kinds of spatial indexing, efficiently calculating bounding boxes first, etc. I'm pretty sure I could more effectively slice this question "by hand", but it'd be nice to be able to just attack it with a simple, even if expensive "brute force" gTouches (or similar spatial over, gRelate questions). I've heard rumor if I'm hitting these issues I should consider (1) a computer buff, since I'm at a measley 6 gigs; and (2) buying processing power elsewhere. Best spots for a starting place for that?

I'm new to R.

  • did you ever find out how to do this?
    – cschwab98
    Jul 19, 2021 at 20:17
  • @cschwab98 In short, yes. I wrote this a LONG time ago as I was just learning R, now I instruct it. I've done similar work a few times now. Using the sf package and tidyverse, I'd first use a small buffer to subset my comparison group if I needed to. And using purrr - which makes this much more manageable - you could all use future or furrr packages to split the work more efficiently across multiple processors. If you have a specific, new question you're welcome to create it and tag. Jul 22, 2021 at 21:02
  • Thank you for the response! I ended up using poly2nb() and map() to retrieve the population counts and sum them in order to calculate the composite population of census tracts.
    – cschwab98
    Jul 25, 2021 at 1:09
  • Ah! If I'm following, you're adding up populations into tracts from... block groups? I'm not sure if you've checked this out, but sf is compatible with tidyverse, so if you have block groups you can group_by() %>% summarize() and you'll get the polygon area unioned as well. And if you're working with populations, highly suggest checking out the tidycensus package. Best wishes. Jul 26, 2021 at 2:35
  • So I'm working on implementing Oka & Wong (2014)'s measures of residential segregation, which requires calculating the composite population of a census tract - this is defined as the population of the focal tract and the population of its neighbors, so unfortunately aggregation wouldn't be sufficient
    – cschwab98
    Jul 26, 2021 at 21:53


Browse other questions tagged or ask your own question.