Find which Shapely Polygon contains an intersection

polygonshapely

Intersection between a line and a Polygon is not inside the Polygon(this). Polygon does not have any interiors. My intention is to find which Polygon contains the intersection, if there are multiple(10) polygons overlapping. I tried projecting on to the polygon, however, it does not work. Is there any way to handle this?

shapely version: 1.8.5.post1 | descartes version: 1.1.0

#!/usr/bin/python3
import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Point, Polygon, LineString
from descartes import PolygonPatch

fig, ax = plt.subplots()

polygon = Polygon([[0.07, 0.0], [0.13, 0.0], [0.13, 0.04], [0.07, 0.04]])
patch = PolygonPatch(polygon, edgecolor='blue', facecolor='none', alpha=1.0)
ax.add_patch(patch)

line = LineString([(0.10,0.03), (0.14,0.03)])
ax.plot(*line.xy, color='black')

print("Intersection points <before> projection, Polygon.contains?")
interesection_points = np.array(polygon.boundary.intersection(line))
ax.plot(interesection_points[0], interesection_points[1], color='blue', marker='h', linestyle='', markersize=15, fillstyle='none', label='intersection before projection')
print (interesection_points, "contains:", polygon.contains(Point(interesection_points)))

print("Intersection points <after> projection, Polygon.contains?")
interesection_points_after_projection_case = polygon.boundary.intersection(line)
project_points = []
if not polygon.contains(interesection_points_after_projection_case):
    nearest_point = np.array(polygon.boundary.interpolate(polygon.boundary.project(interesection_points_after_projection_case)))
    project_points.append(nearest_point)
    print(nearest_point, "contains:", polygon.contains(Point(nearest_point)))

project_points = np.asarray(project_points)
ax.plot(project_points[:,0], project_points[:,1], color='green', marker='^', linestyle='', markersize=15, fillstyle='none', label='intersection after projection')

ax.legend()
plt.show()

Intersection points <before> projection, Polygon.contains?
[0.13 0.03] contains: False

Intersection points <after> projection, Polygon.contains?
[0.13 0.03] contains: False

Best Answer

The definition of the "contains" predicate is as follows (source):

Returns True if no points of other lie in the exterior of the object and at least one point of the interior of other lies in the interior of object.

At least one point needs to be in the interior of the other is the issue here. The point you are comparing is on the boundary of the polygon, not in the interior of the polygon.

So, for this example if the point being on the boundary is already OK you could use e.g.

intersects: point is in the interior or on the boundary of the polygon
touches: point is on the boundary of the polygon

Related Solutions

[GIS] Shapely unable to tell if polygon contains point

If we examine your polygon:

polygon = shapefile_record['geometry']
print polygon.bounds
(77.84476181915733, 30.711096140487314, 78.59476181915738, 31.28199614048725)

From Shapely manual, object.bounds:

Returns a (minx, miny, maxx, maxy) tuple (float values) that bounds the object.

Here minx = 77.84476181915733, miny = 30.711096140487314

= here, min longitude = 77.844... and min latitude = 30.711...

Setting your point as Point(30.9787, 78.8900) or Point(latitude, longitude), you reverse the coordinates:

your_point = Point(30.9787, 78.8900) # = Point(y, x)
xy_point = Point(78.8900,30.9787) # = Point(x,y)
polygon.contains(xy_point)
False

Which is correct:

enter image description here

So, the script with your example:

from shapely.geometry import Point, shape
point = Point(78.8900,30.9787)
fc = fiona.open("AC_All_Final.shp")
print fc.schema
{'geometry': 'Polygon', 'properties': OrderedDict([(u'DIST_NAME', 'str:254'), (u'Assem_No', 'int:4'), (u'AREA', 'float:18.3'), (u'PERIMETER', 'float:18.3'), (u'INDIAASSEM', 'float:11'), (u'INDIAASS', 'float:11'), (u'NO', 'float:11'), (u'ST_CODE', 'str:254'), (u'AC_NAME', 'str:254'), (u'AC_TYPE', 'str:254'), (u'PC_NO', 'float:11'), (u'AC_NO', 'float:19.11'), (u'AC_HNAME', 'str:254'), (u'PARTY', 'str:254'), (u'CODE_NO', 'str:254'), (u'Longitude', 'float:19.17'), (u'Latitude', 'float:19.17'), (u'State', 'str:50'), (u'PC_NAME', 'str:25'), (u'PC_TYPE', 'str:3'), (u'PC_CODE', 'str:10'), (u'PC_NO_1', 'int:4')])}
from feature in fc:
    if shape(feature['geometry']).contains(point):
          print feature['properties']['DIST_NAME'], feature['properties']['Longitude'], feature['properties']['Latitude']

Uttarkashi 78.8900584624 30.9787742252

Result:

enter image description here

[GIS] Which method to find intersection between lines with polygon

I tried to take your question becuase it is an interesting one to be solved with an SLQ statement.

Reading your tittle the case was either a linestring & polygon or a linestring & multipolygon type on the picture below (Source: Introduction to PostGIS):

But it seems you are mixing concepts. I GIS there is a limited number of geometry types compatred with the CAD paradigma. In GIS there are fewer vector types, the main ones are: points, lines and polygons (or nodes, edges and faces)

In GIS a point is an entity that had zero dimension (0D), it has only a location, a line has 1-D and it is made out of a set of points and a polygon is a 2D entity maede aout of a set of lines.

This picture will show the main three types of geometry entities (Source: A Gentle Introduction to GIS):

There are more types, but main differences with the CAD paradigma are in the confusion between lines, polylines and polygons. A polyline in CAD is just a line type in GIS, and a closed polyline in CAD could be read as line or a polygon type depending on the success of the importing tool.

Then we have the difference of the layer paradigma. CAD applications might store any kind of entity types in a layer, while GIS applications store only one type per layer.

Back to your question, it does not matter what the layer is called, it can store only one kind of data, and that can be read for the icon.

Line and polygon has only line type entities and Muffe and Arrow polygons.

My main suggestion now to you now is to read about GIS to get used to this new GIS paradigma about entity times and storage organization (layers). I highly recomend A Gentle Introduction to GIS Link to .pdf. Then, I suggest you to choose another file format than .dxf for your work in GIS, dxf does not fit very well with any gis analysis.

Best Answer

Related Solutions

[GIS] Shapely unable to tell if polygon contains point

[GIS] Which method to find intersection between lines with polygon

Related Question