intersection not returning true on polygons that have 1 point in common

[davidanderson-1701]

I am processing a largish spatial dataset wherein group together various polygons based on attributes, then do a spatial union. For the case where the polygons touch at just point, the union is creating two polygons and not one.
Checking the data for individual cases, the polygons neither intersect nor are disjoint.
Here is the data in wkt
array([ 'POLYGON ((503988.2057999996500000 5214800.2800999992000000, 503838.2706000004000000 5214797.4101000000000000, 503835.3997999997800000 5214947.3000000007000000, 503985.3349000001300000 5214950.1699999999000000, 503988.2057999996500000 5214800.28009999

array([ 'POLYGON ((504141.0118000004400000 5214653.2601999994000000, 503991.0766000002600000 5214650.3902000003000000, 503988.2057999996500000 5214800.2800999992000000, 504138.1409000000000000 5214803.1501000002000000, 504141.0118000004400000 5214653.2601999994000000))'], dtype=object)

The common point is :
503988.2057999996500000 5214800.2800999992000000

Here are the commands that I tried:
spatial_data[spatial_data['SIMPLE_ID']==4309842].intersects(spatial_data[spatial_data['SIMPLE_ID']==4305747])
Out[20]:
77007 False
77681 False

spatial_data[spatial_data['SIMPLE_ID']==4309842].disjoint(spatial_data[spatial_data['SIMPLE_ID']==4305747])
Out[21]:
77007 False
77681 False

spatial_data[spatial_data['SIMPLE_ID']==4309842].touches(spatial_data[spatial_data['SIMPLE_ID']==4305747])
Out[22]:
77007 False
77681 False

Loading the data into Spatialite to test (since it uses the same geos library) does result in two polygons intersecting.

[jorisvandenbossche]

@davidanderson-1701 Would it be possible to provide a reproducible, self-contained example? (runnable code). Because it is difficult to say what is going on without knowing exactly what is inside spatial_data.

Anyway, two remarks:

• Polygons that share only one point in any case do 'intersect', given the following simple shapely example (the geopandas intersects methods just dispatches to this under the hood):

  In [5]: p1 = shapely.geometry.Polygon([(0,0), (1,0), (1,1), (0,1)])
  
  In [8]: p2 = shapely.geometry.Polygon([(1,1), (2,1), (2,2), (1,2)])
  
  In [10]: p1.intersects(p2)
  Out[10]: True
  

• When providing a GeoDataFrame/Series to GeoDataFrame.intersects, it will first align the calling object and provided object (match rows with matching index). So if in the above eg spatial_data[spatial_data['SIMPLE_ID']==4305747] does select a single element (but so a single row GeoDataFrame), it will first align and so maybe (depending on the data and what you want to obtain) does not what you expect it to do.
  To illustrate with the same simple example:

  In [22]: s1 = geopandas.GeoSeries([p1])
  
  In [23]: s2 = geopandas.GeoSeries([p2])
  
  In [24]: s1.intersects(s2)
  Out[24]:
  0    True
  dtype: bool
  
  # giving it another index, s1 and s2 will now first align (and in this case return False
  # for each index, as there are no polygons that intersect and have matching index)
  In [25]: s2 = geopandas.GeoSeries([p2], index=[1])
  
  In [26]: s1.intersects(s2)
  Out[26]:
  0    False
  1    False
  dtype: bool
  

  In the last case, you can eg do

  In [29]: s1.intersects(s2.geometry.iloc[0])
  Out[29]:
  0    True
  dtype: bool
  

  if you have a one row series/dataframe to select the single geometry object (in that case, it will not align on the calling series/dataframe but return intersects with all rows of the calling object)

[davidanderson-1701]

I've attached a shapefile with some representative data.
Google Drive

The process I used is

`def slink_union(slink_list):
polygons = gpd.GeoSeries(spatial_data[spatial_data['slink'].isin(slink_list)].unary_union)
return polygons

spatial_data = gpd.read_file(r'F:\workspace\lb_test.shp')
slink_list = [4212048,4216113,4220178]
spatial_data['slink']=spatial_data['SIMPLE_ID']
poly=slink_union(slink_list)
poly_ex=poly.explode()
`

The three polygons are arranged like this:
A
B
C

So one corner touching and one line touching.

The code produces two polygons. One that is a single and one that is the
two line touching polygons merged into one polygon.

Really the issue is that unary_union is not seeing those polygons as intersecting and thus not combining them.
Cross testing in Spatialite shows the same issue. I would guess that this is a GEOS issue since the two apps (geopandas and spatialite) both use that library.

On Mon, May 16, 2016 at 1:05 PM, Joris Van den Bossche <
notifications@github.com> wrote:

@davidanderson-1701 https://github.com/davidanderson-1701 Would it be
possible to provide a reproducible, self-contained example? (runnable
code). Because it is difficult to say what is going on without knowing
exactly what is inside spatial_data.

Anyway, two remarks:

Polygons that share only one point in any case do 'intersect', given
the following simple shapely example (the geopandas intersects methods
just dispatches to this under the hood):

In [5]: p1 = shapely.geometry.Polygon([(0,0), (1,0), (1,1), (0,1)])

In [8]: p2 = shapely.geometry.Polygon([(1,1), (2,1), (2,2), (1,2)])

In [10]: p1.intersects(p2)
Out[10]: True

When providing a GeoDataFrame/Series to GeoDataFrame.intersects, it
will first align the calling object and provided object (match rows with
matching index). So if in the above eg
spatial_data[spatial_data['SIMPLE_ID']==4305747] does select a single
element (but so a single row GeoDataFrame), it will first align and so
maybe (depending on the data and what you want to obtain) does not what you
expect it to do.
To illustrate with the same simple example:

In [22]: s1 = geopandas.GeoSeries([p1])

In [23]: s2 = geopandas.GeoSeries([p2])

In [24]: s1.intersects(s2)
Out[24]:
0 True
dtype: bool

giving it another index, s1 and s2 will now first align (and in this case return False

for each index, as there are no polygons that intersect and have matching index)

In [25]: s2 = geopandas.GeoSeries([p2], index=[1])

In [26]: s1.intersects(s2)
Out[26]:
0 False
1 False
dtype: bool

In the last case, you can eg do

In [29]: s1.intersects(s2.geometry.iloc[0])
Out[29]:
0 True
dtype: bool

if you have a one row series/dataframe to select the single geometry
object (in that case, it will not align on the calling series/dataframe but
return intersects with all rows of the calling object)

—
You are receiving this because you were mentioned.
Reply to this email directly or view it on GitHub
#317 (comment)