Python GeoPandas – For-Loop with Geometry Constructor

geopandaspythonshapely

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?

Best Answer

I don't know your data (duplicate of the question in reddit/learnPython), but from your script, it is not necessary to iterate over rows in the GeoDataFrame (see How to iterate over rows in a DataFrame in Pandas):

First select the data in a temporary GeoDataFrame (no need of if:else:pass)

temp = parcel_bldg_clip[parcel_bldg_clip['can_it_adu']  == 'Y']

Create the buffer in this temporary GeoDataFrame

a) use list_comprehension

pl = [polylabel(geom) for geom in temp.geometry]
buf = [geom.buffer(4, cap_style=3) for geom in pl]
temp['geometry'] = buf

b) or use .apply

from shapely.ops import polylabel
temp['geometry'] = temp.apply(lambda row:polylabel(row.geometry), axis=1)
temp['geometry'] = temp.geometry.buffer(4, cap_style=3)

Check if the hypothetical building polygon (temporary GeoDataFrame) is within the parcel (original GeoDataFrame) or rotate it (here with a for-loop but outside the GeoDataFrame, no .iterrows(), .itertuples() or .iteritems(), and modify the original GeoDataFrame

 for i,(geom1,geom2) in enumerate(zip(parcel_bldg_clip.geometry,temp.geometry)):
    if geom1.contains(geom2) :
        # modify the original GeoDataFrame
        parcel_bldg_clip.loc[i,['will_it_ad']] = 'Y'
    else:
        for angle in range(0,360,10):
            # rotate from 0 to 360° with steps of 10°
            geom2= gemom2.rotate(angle, origin='centroid') 
                if geom1.contains(geom2) :
                    parcel_bldg_clip.loc[i,['will_it_ad']] = 'Y'
Related Question