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.
require(maptools)
require(sp)
require(rgdal)
require(rgeos)
download.file("ftp://ftp2.census.gov/geo/pvs/tiger2010st/37_North_Carolina/37/tl_2010_37_county10.zip", destfile="NC.counties.zip")
unzip("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(NC.counties.sp)
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
unzip("NC.counties.zip")
NCpolys.spdf = readOGR(dsn=".", "tl_2010_37_tabblock10", stringsAsFactors = F) #read blocks
#Try the above again. gTouches, eek.
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.
