1

I have a GeoDataFrame in which the geometry column consists of polygons.

Zones = 

  id                                           geometry
0   A1  POLYGON ((119.82334 28.35047, 119.79008 28.350...
1   A2  POLYGON ((138.79141 26.48802, 141.92247 26.261...
2   A3  POLYGON ((65.42208 48.14785, 65.42208 48.14785...
3   A4  POLYGON ((152.14557 25.88444, 155.50298 26.186...

and the geometry column is actually a geometry-type obejct:

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 103 entries, 0 to 102
Data columns (total 2 columns):
 #   Column    Non-Null Count  Dtype   
---  ------    --------------  -----   
 0   id        103 non-null    object  
 1   geometry  103 non-null    geometry
dtypes: geometry(1), object(1)
memory usage: 1.7+ KB
None

Now, I have used a function, subdivide_polygon that divides these polygons into small square (I post the function at the end of my question not to clutter the text). The manner in which I apply the function is the following to the first geometry:

polygon = geostore_zones_geometry_polygon.geometry.loc[0]
new = subdivide_polygon(polygon,1, 1)

Thee polygon looks like this:

enter image description here

and the subdivided one like this:

enter image description here

The problem if that it fails for some of the polygons with the error:

AttributeError: No geometry data set yet (expected in column 'geometry'.)

My first thought was that there was something inherent to the polygons that might cause this, but those for which the method fails are similar to the ones for which it work.

My second thought was that even thought the geometry column has dtype geometry, some rows might not. So, I checked that by making every row into individual dataframe and printing the info:

for i in range(0,len(geostore_zones_geometry_polygon)):
    df =  geostore_zones_geometry_polygon.loc[:i,:]
    print(df.info())

All of them are geometry type objects!

So, what is wrong?

Data I post here the first three rows of my dataframe. The method works well on the first two rows, not on the third (the only thing that is different is the number of decimals in some coordinate for the row that the method fails on, but surely, that can't be it!):

Zones = 

  id                                           geometry
0   A1  POLYGON ((119.82334 28.350468, 119.79008 28.350468, 122.85067 28.084328, 122.85067 44.851055, 119.92314 44.717983, 119.82334 28.350468))
1   A2  POLYGON ((138.79141 26.488018, 141.92247 26.261677, 141.73386 44.897156, 138.79141 44.897156, 138.79141 26.488018))
2   A3  POLYGON ((65.42207999999999 48.14785, 65.42207999999999 48.14785, 46.635708 51.165745, 46.597984 47.657444, 68.4777 44.0737, 65.42207999999999 48.14785))

Method

import numpy as np
import geopandas
from shapely.ops import split
from shapely.geometry import MultiPolygon, Polygon, Point, MultiPoint

def get_rect_from_geom(G, r=1):
    
    coordinate_arrays = G.exterior.coords.xy
    coordinates = list(
        zip(
            [np.round(c, r) for c in coordinate_arrays[0]],
            [np.round(c, r) for c in coordinate_arrays[1]],
        )
    )
    Rectangles = []
    for c1 in coordinates:
        Coords1 = [a for a in coordinates if a != c1]
        for c2 in Coords1:
            Coords2 = [b for b in Coords1 if b != c2]
            x1, y1 = c1[0], c1[1]
            x2, y2 = c2[0], c2[1]
            K1 = [k for k in Coords2 if k == (x1, y2)]
            K2 = [k for k in Coords2 if k == (x2, y1)]
            if (len(K1) > 0) & (len(K2) > 0):
                rect = [list(c1), list(K1[0]), list(c2), list(K2[0])]
                Rectangles.append(rect)
    return Rectangles


def get_squares_from_rect(rect, side_length):
    
    y_list = [r[1] for r in rect]
    x_list = [r[0] for r in rect]
    y1 = min(y_list)
    y2 = max(y_list)
    x1 = min(x_list)
    x2 = max(x_list)
    width = x2 - x1
    height = y2 - y1

    xcells, ycells = int(np.round(width / side_length)), int(
        np.round(height / side_length)
    )

    yindices = np.linspace(y1, y2, ycells + 1)
    xindices = np.linspace(x1, x2, xcells + 1)
    horizontal_splitters = [
        LineString([(x, yindices[0]), (x, yindices[-1])]) for x in xindices
    ]
    vertical_splitters = [
        LineString([(xindices[0], y), (xindices[-1], y)]) for y in yindices
    ]
    result = Polygon(rect)
    for splitter in vertical_splitters:
        result = MultiPolygon(split(result, splitter))
    for splitter in horizontal_splitters:
        result = MultiPolygon(split(result, splitter))
    square_polygons = list(result)
    return [np.stack(SQPOLY.exterior.coords.xy, axis=1) for SQPOLY in square_polygons]


def round_centroid(g, r=1):
   
    C = g.centroid.coords.xy
    return (np.round(C[0][0], r), np.round(C[1][0], r))


def subdivide_polygon(g, side_length, r):
    
    #   Number of squares required.
    num_squares_reqd = g.area // (side_length ** 2)

    # Some of these combinations can be dropped to improve performance.
    Rectangles = []
    Rectangles.extend(get_rect_from_geom(g))
    Rectangles.extend(get_rect_from_geom(g.minimum_rotated_rectangle))
    Rectangles.extend(get_rect_from_geom(g.convex_hull))

    Squares = []
    for r in range(len(Rectangles)):
        rect = Rectangles[r]
        Squares.extend(get_squares_from_rect(rect, side_length=side_length))  ### added r
    SquarePolygons = [Polygon(square) for square in Squares]
    GDF = geopandas.GeoDataFrame(SquarePolygons).rename(columns={0: "geometry"})
    GDF.loc[:, "centroid"] = GDF.geometry.apply(round_centroid, r=r)
    GDF = GDF.drop_duplicates(subset=["centroid"])

    wgeoms = GDF[GDF.within(g)].geometry.values
    igeoms = GDF[GDF.intersects(g)].geometry.values

    w = abs(num_squares_reqd - len(wgeoms))
    i = abs(num_squares_reqd - len(igeoms))

    print(w, i)
    if w <= i:
        return wgeoms
    else:
        return igeoms

Full traceback

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-225-b60f84617202> in <module>
----> 1 new = subdivide_polygon(polygon,1, 1)

<ipython-input-176-8a2e9714702b> in subdivide_polygon(g, side_length, r)
     84     SquarePolygons = [Polygon(square) for square in Squares]
     85     GDF = geopandas.GeoDataFrame(SquarePolygons).rename(columns={0: "geometry"})
---> 86     GDF.loc[:, "centroid"] = GDF.geometry.apply(round_centroid, r=r)
     87     GDF = GDF.drop_duplicates(subset=["centroid"])
     88 

~\Anaconda3\envs\conda-qgis\lib\site-packages\pandas\core\generic.py in __getattr__(self, name)
   5463             if self._info_axis._can_hold_identifiers_and_holds_name(name):
   5464                 return self[name]
-> 5465             return object.__getattribute__(self, name)
   5466 
   5467     def __setattr__(self, name: str, value) -> None:

~\Anaconda3\envs\conda-qgis\lib\site-packages\geopandas\geodataframe.py in _get_geometry(self)
    184     def _get_geometry(self):
    185         if self._geometry_column_name not in self:
--> 186             raise AttributeError(
    187                 "No geometry data set yet (expected in"
    188                 " column '%s'.)" % self._geometry_column_name

AttributeError: No geometry data set yet (expected in column 'geometry'.)
8
  • 1
    Possible that its the polygon itself that has a problem. Is there another function that you could call which would validate the geometry as a valid polygon geometry?
    – nr_aus
    Jun 14, 2021 at 6:58
  • @nr_aus Yes I have checked. point.within(polygon) works just fine. It would if polygonwasn't one. Jun 14, 2021 at 7:40
  • Debug your script because with your solution, the problem comes from get_rect_from_geom(G, r=1) -> Rectangle is empty for geometry A3
    – gene
    Jun 14, 2021 at 8:55
  • 1
    note that A3's 1st coordinate is repeated. With 14 decimals, it is possible that floating point arithmetic fails to identify equality between them, and so your function may not behave as intended. So you can try to 1) remove the duplicated point and 2) keep a meaningful precision
    – JGH
    Jun 14, 2021 at 13:13
  • 1

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.