Reputation: 5
I am using shapely to find if a LineString and Polygon overlap.
I have defined the polygon based on the map area that I need to find if points are overlapping in:
polygon_map = Polygon([(-126.03561599522513,60.08405276856493),(-65.91842849522513,60.08405276856493),(-65.91842849522513,76.84958092750016),(-126.03561599522513,76.84958092750016),(-126.03561599522513,60.08405276856493)])
For the LineString, I have a long list of coordinates in two columns, x and y. I have taken the maximum and minimum coordinates of x and y to generate a LineString.
line = LineString([(32.823101,87.988993),(-153.01468,30.001368)])
When I plot these on a map, they overlap (as expected)
m = folium.Map([61.08405, -66.918], zoom_start=3, tiles='cartodbpositron')
folium.GeoJson(line).add_to(m)
folium.GeoJson(polygon_map).add_to(m)
folium.LatLngPopup().add_to(m)
m
[Image of map created showing intersecting polygon and linestring]
However, when I do:
line.overlaps(polygon_map)
It returns false, and I can't work out why.
I have simplified the LineString to only include the minimum and maximum coordinates as I have hundreds of coordinates in my original dataframe and I'm worried it will take too long to loop through each set of coordinates. I haven't used Shapely before so I'm not sure if this is why it isn't working.
Upvotes: 0
Views: 909
Reputation: 31226
This is all down to geographic projections. As pure cartesian geometry without accounting for curvature of the earth they do not overlap. (See below image). shapely has no knowledge of geographic projections, it is pure cartesian geometry. As cartesian geometric objects this polygon and LineString do not overlap.
Only after setting the CRS does folium show these geometries overlapping.
import geopandas as gpd
from shapely.geometry import Polygon, LineString
polygon_map = Polygon([(-126.03561599522513,60.08405276856493),(-65.91842849522513,60.08405276856493),(-65.91842849522513,76.84958092750016),(-126.03561599522513,76.84958092750016),(-126.03561599522513,60.08405276856493)]) # fmt: skip
line = LineString([(32.823101, 87.988993), (-153.01468, 30.001368)])
gdf = gpd.GeoDataFrame(geometry=[polygon_map, line])
gdf.explore(height=300, width=300)
# gdf.set_crs("epsg:4386").explore(height=300, width=300)
Upvotes: 1