1

I have this situation in which I want to detect elements that ought to be contained on the boundary of a geometrical feature, but because of a variety of reason these point-wise objects can be "seen" within or outside the geometry. The within part is not an issue since I wish to use "contains" as a method. So, what I do is that I buffer the geometries I have in order to "catch" elements that are just outside the boundary. But this is kind of an issue as well. Indeed, buffering might imply that an element belongs to two geometries when they actually can only be in one.

Here is an example (albeit very naïve, it still illustrate the problem).

import geopandas
import pandas as pd
%matplotlib inline
import matplotlib.pyplot as plt
from shapely.geometry import Point
from geopandas import datasets, GeoDataFrame, read_file
from geopandas.tools import overlay
from shapely.geometry import Polygon, LineString, Point
s = geopandas.GeoSeries(
    [
        Point(1.05, 1.1),
        Point(0.1, 1.1),
    ],
)
s2 = geopandas.GeoSeries(
    [
        Polygon([(0, 0), (1, 1), (0, 1)]),
        Polygon([(1.25, 1), (2, 2), (1.25, 2)]),
        
    ],
    index=range(1, 3),
)

envgdf1 = geopandas.GeoDataFrame(geometry=gpd.GeoSeries(s2))
envgdf2 = geopandas.GeoDataFrame(geometry=gpd.GeoSeries(s))

envgdf1 = envgdf1.rename_geometry('Object')
envgdf2 = envgdf2.rename_geometry('Point')

enter image description here

Picture the triangle as my objects and the black points as these object that actually ought to be on either one of these triangle. As I explained above, I buffer the triangle in order to ingulf these points:

envgdf1['buffered_object'] = envgdf1['Object'].buffer(0.3)
envgdf1 = envgdf1.set_geometry('buffered_object')

Which results in

enter image description here

Seen as a dataframe:

df = gpd.sjoin(envgdf1,envgdf2, how="inner", op='contains')


 Object  \
1  POLYGON ((0.00000 0.00000, 1.00000 1.00000, 0....   
1  POLYGON ((0.00000 0.00000, 1.00000 1.00000, 0....   
2  POLYGON ((1.25000 1.00000, 2.00000 2.00000, 1....   

                                     buffered_object  index_right  
1  POLYGON ((-0.30000 0.00000, -0.30000 1.00000, ...            1  
1  POLYGON ((-0.30000 0.00000, -0.30000 1.00000, ...            0  
2  POLYGON ((0.95000 1.00000, 0.95000 2.00000, 0....            0  

you can notice that there are two rows for the same object index_right = 0, meaning that the point is in both polygons.

Is there a way to deal with this in a good way? That is:

  1. Ignoring anything where the inclusion is ambiguous
  2. Even better: "If the point is closer to one of the geometries, then it belong to that one"
  3. Any other method?
1
  • @BERA Yes, I want to decide with polygon the point actually should belong to. The aim of buffering is to include the point in some polygon. The problem is that doing so it can belong to any of the two polygon. I want to be able to decide if it belong to the top-right one or botton-left one. I changed the picture slightly. In it the point between the triangles should belong to the bottom right triangle, but buffering makes it impossible to make that call. Oct 12, 2021 at 11:48

1 Answer 1

1

So, I found a way to answer my question. There might be something better than that and I would be grateful for an answer that makes mine seem cumbersome.

The idea is to create a buffer zone around the polygons in such a way that points will be absorded by the buffered poygon, then sjoin BUT keeping the points. A usual sjoin will disregard the geometry which one wants contained. So here we go:

Instead of

df = gpd.sjoin(envgdf1,envgdf2, how="inner", op='contains')

we'll ensure keeping the points by using savedgeom:

envgdf2['savedgeom'] = envgdf2.geometry
df = gpd.sjoin(envgdf1,envgdf2, how="inner", op='contains')

which results in:

 Object  \
1  POLYGON ((0.00000 0.00000, 1.00000 1.00000, 0....   
1  POLYGON ((0.00000 0.00000, 1.00000 1.00000, 0....   
2  POLYGON ((1.25000 1.00000, 2.00000 2.00000, 1....   

                                     buffered_object  index_right  \
1  POLYGON ((-0.30000 0.00000, -0.30000 1.00000, ...            1   
1  POLYGON ((-0.30000 0.00000, -0.30000 1.00000, ...            0   
2  POLYGON ((0.95000 1.00000, 0.95000 2.00000, 0....            0   

                 savedgeom  
1  POINT (0.10000 1.10000)  
1  POINT (1.10000 1.10000)  
2  POINT (1.10000 1.10000)  

Now, we can determine the distance of each point to the different polygons:

df.crs = "EPSG:4326"
df['dist'] = gpd.GeoSeries(df['Object']).distance(gpd.GeoSeries(df['savedgeom']))

which gives:

Object  \
1  POLYGON ((0.00000 0.00000, 1.00000 1.00000, 0....   
1  POLYGON ((0.00000 0.00000, 1.00000 1.00000, 0....   
2  POLYGON ((1.25000 1.00000, 2.00000 2.00000, 1....   

                                     buffered_object  index_right  \
1  POLYGON ((-0.30000 0.00000, -0.30000 1.00000, ...            1   
1  POLYGON ((-0.30000 0.00000, -0.30000 1.00000, ...            0   
2  POLYGON ((0.95000 1.00000, 0.95000 2.00000, 0....            0   

                 savedgeom      dist  
1  POINT (0.10000 1.10000)  0.100000  
1  POINT (1.10000 1.10000)  0.141421  
2  POINT (1.10000 1.10000)  0.150000  

Since the goal was to match a point to the closest polygon we can simply do

df = df.sort_values("dist", ascending=True).groupby(["index_right"]).first().reset_index()

which gives:

 index_right                                             Object  \
0            0  POLYGON ((0.00000 0.00000, 1.00000 1.00000, 0....   
1            1  POLYGON ((0.00000 0.00000, 1.00000 1.00000, 0....   

                                     buffered_object                savedgeom  \
0  POLYGON ((-0.30000 0.00000, -0.30000 1.00000, ...  POINT (1.10000 1.10000)   
1  POLYGON ((-0.30000 0.00000, -0.30000 1.00000, ...  POINT (0.10000 1.10000)   

       dist  
0  0.141421  
1  0.100000  

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.