If it's not too much, I have a specific question that leads to a more general one: I'd love thoughts on both if it's possible. 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 metersNCpolys.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) #
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 metersNCpolys.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) #