The first round reveals block groups that are:
However, this omits adjacent blocks that could be spatially merged into one or more organizational units. To identify these other areas, we isolate clusters of majority black blocks with fewer than 50 BHH, and dissolve them into single units.
However, this often results in clusters with far more than 100 BHH: too much for block leaders to cover. So, we must spatially aggregate these blocks up to a point where total black population hits 100 HH.
The workflow is as follows:
import pandas as pd
import geopandas as gpd
%matplotlib inline
#Get block data
fcBlocksAll = gpd.read_file('./data/WAKE/BlockMece.shp')
Subset the blocks, keeping only those with fewer than 50 black households. This will create clusters of blocks separated by other blocks with
#Subset blocks with fewer than 50 black households
fcBlocksSubset = fcBlocksAll[(fcBlocksAll.BlackHH < 50) & (fcBlocksAll.BlackHH > 10)].reset_index()
#Dissolve adjancet blocks
fcClusters = gpd.GeoDataFrame(geometry = list(fcBlocksSubset.unary_union))
#Add a static ID field to the geodataframe
fcClusters['ID'] = fcClusters.index
#Copy over crs to new file
fcClusters.crs = fcBlocksSubset.crs
Assign new columns to the block features. ID
will be the cluster to which the block belongs, and OrgID
will be it's assigned organizaitional unit ID.
#Spatially join the dissolved ID to the subset blocks layer
fcBlockSubset2 = gpd.sjoin(fcBlocksSubset,fcClusters,how='left',op='within').drop("index_right",axis=1)
#Initialize the field to contain new organizational unit IDs
fcBlockSubset2['OrgID'] = 0
Now we compute the total number of BHH found within each cluster. From that we can identify which clusters are good as is (those with between 50 and 100 aggregate BHH), those which remain untenable (fewer than 50 aggregate BHH), and those that need to be divided up (more than 100 aggregate BHH).
#Compute total BHH for the dissolved blocks and add as attribute to bloc
sumHH = fcBlockSubset2.groupby('ID').agg({'BlackHH':'sum'})
#Join total aggregate BHH to the cluster geoframe
fcClusters2=pd.merge(fcClusters,sumHH,left_index=True,right_index=True)
#Save clusters with between 50 and 100 HH to a a new geoframe
fcClusters_keep1 = fcClusters2[(fcClusters2.BlackHH > 50) & (fcClusters2.BlackHH <= 100)].reset_index()
fcClusters_keep1.shape[0]
26
Now to select new blocks with > 100 HH and break them up.
#Find IDs of dissolved blocks with HH > 100
fcTooBig = fcClusters2.query('BlackHH > 100')
fcTooBig.shape[0]
30
#We'll iterate through each
for idx in fcTooBig.ID:
#Get the cluster ID
clusterID = fcTooBig.loc[idx,"ID"]
#Get the blocks with that ID
fcCBlocks = fcBlockSubset2[(fcBlockSubset2.ID == clusterID) & (fcBlockSubset2.OrgID == 0)].reset_index()
#Add the X coordinate as a column
fcCBlocks['X'] = fcCBlocks.geometry.centroid.x
#Start with the feat with the min X
fcNbrs = fcCBlocks[fcCBlocks.X == fcCBlocks.X.max()]
#Get the number of aggregate BGG in the selection
BHH = fcNbrs.BlackHH.sum()
geom = fcNbrs.geometry.unary_union
#Increase the selection
iterX = 0 #Catch to avoid infinite loops
while BHH < 100 and iterX < 100:
#Find the blocks that touch
fcNbrs = fcBlockSubset2[(fcBlockSubset2.intersects(geom)) & #Select blocks that are adjacent
(fcBlockSubset2.BlackHH < 50) & #...that have fewer than 50 BHH
(fcBlockSubset2.BlackHH > 10) #...but have at least 10 BHH
]
BHH = fcNbrs.BlackHH.sum()
geom = fcNbrs.geometry.unary_union
iterX += 1
#Select blocks and assign its OrgID
fcBlockSubset2.loc[fcBlockSubset2.intersects(geom),'OrgID'] = idx+1000
fcClusters_keep2 = fcBlockSubset2.query('OrgID > 100').dissolve(by='OrgID',aggfunc={'BlackHH':'sum'})
fcClusters_keep2['ID']=fcClusters_keep2.index
fcAll = fcClusters_keep1.append(fcClusters_keep2,ignore_index=True,sort=False).reset_index()
fcAll.to_file("./scratch/OrgUnitsX.shp")
fcOrig = fcBlocksAll.loc[fcBlocksAll.BlackHH > 50,['BlackHH','geometry']]
fcOrig['ID']
fcOrig.head()
BlackHH | geometry | |
---|---|---|
0 | 70.072829 | POLYGON ((-78.33744 35.821984, -78.337457 35.8... |
4 | 52.784615 | POLYGON ((-78.572563 35.811895, -78.572035 35.... |
5 | 62.686567 | POLYGON ((-78.65722199999999 35.782297, -78.65... |
9 | 137.087805 | POLYGON ((-78.60934899999999 35.75329, -78.609... |
10 | 58.722581 | POLYGON ((-78.65531299999999 35.87068, -78.654... |