I am trying to create a tool to find tax parcels that have enough open space to build a new structure of a certain size.
I start with a GeoJSON of tax parcels and a GeoJSON of building footprints and create two separate geodataframes for them. I then intersected the parcels with the building footprints to get one polygon geodataframe that just covers open space using gpd.overlay - all good so far.
My trouble is constructing the actual function to test each of the parcels for sufficient size/space. In my head, it should go something like this:
- Check if the area of the parcel is bigger than the hypothetical building. If the parcel is not big enough, flag the parcel and exit
- If the area is sufficient, find the Pole of Inaccessibility for each parcel
- Construct a hypothetical building polygon with a square buffer centered on the parcel's pole of inaccessibility
- Check if the hypothetical building polygon is within the parcel. If it is, flag the parcel and exit
- If the hypothetical building is not fully within the parcel, rotate the polygon 10 degrees and recheck if it falls within the parcel. If not, keep rotating until a fit is found or a 360-degree rotation has occurred. Either way, flag appropriately and exit.
Here's the code I've written so far for this part (it's no good, but I figure it's a start):
# iterate over each parcel in the parcel_bldg_clip geodataframe and check if the value of parcel_area minus clip_area is greater than 65
for index, row in parcel_bldg_clip.iterrows():
if row['can_it_adu'] == 'Y': # if the parcel has a minimum of 65 sq meters of land area (previously calculated)
parcel_pole = polylabel(row['geometry']) # find the pole of inaccessibility of the parcel
parcel_pole_gdf = gpd.GeoDataFrame(geometry=[test_pole], crs=3857) # type: ignore # create a geodataframe from the pole
parcel_pole_sq_buffer = parcel_pole_gdf.buffer(4, cap_style=3) # 4 meter RADIUS buffer with square cap style to create a square with 8 meter sides (approx 26 feet). This is the hypothetical building footprint
if parcel_pole_sq_buffer.within(row[row['geometry']]): # test if the hypothetical building footprint is within the parcel
row['will_it_adu'] = 'Y' # if yes, then the parcel can have an building
else: # if the hypothetical building footprint is not within the parcel, we need to test if rotating the hypothetical building footprint will make it fit
rotation = 0 # initialize rotation counter to 0
while rotation < 360: # while the rotation counter is less than 360 degrees
parcel_pole_sq_buffer = parcel_pole_sq_buffer.rotate(10, origin='centroid') # rotate the hypothetical building footprint 10 degrees
if parcel_pole_sq_buffer.within(row['geometry']): # re-test if the hypothetical building footprint is within the parcel
row['will_it_adu'] = 'Y' # if yes, then the parcel can have an building
break # exit the while loop
rotation += 10 # increment the rotation counter by 10 degrees
else:
pass
The current error I get with this section of my code is:
TypeError: 'POLYGON ((-85...15))' is an invalid key
During handling of the above exception, another exception occurred:
InvalidIndexError: POLYGON ((-85...15))
One thing I keep coming across are recommendations not to use a for-loop with pandas - however, it seems to me that it's my only option because I'm working with data outside of the geodataframe that is referencing information inside the gdf? Any thoughts on that, as well as generally how I can get this working?