12

I am trying to use Shapely's within function to do a 'spatial join' of a LineString and a Point file (FYI, the point file was generated using the interpolate function on the LineString). Problem is - nothing is being returned.

# this condition is never satisfied
if point.within(line):
    # here I write stuff to a file

where:

point = POINT (-9763788.9782693591000000 5488878.3678984242000000)
line = LINESTRING (-9765787.998118492 5488940.974948905, -9748582.801636808 5488402.127570709)

What am I missing?

Georgy
  • 12,464
  • 7
  • 65
  • 73
user14696
  • 657
  • 2
  • 10
  • 30

1 Answers1

24

There are floating point precision errors when finding a point on a line. Use the distance with an appropriate threshold instead.

from shapely.geometry import Point, LineString

line = LineString([(-9765787.9981184918, 5488940.9749489054), (-9748582.8016368076, 5488402.1275707092)])
point = Point(-9763788.9782693591, 5488878.3678984242)

line.within(point)  # False
line.distance(point)  # 7.765244949417793e-11
line.distance(point) < 1e-8  # True
Mike T
  • 41,085
  • 18
  • 152
  • 203
  • 1
    Hi. How can we choose the appropriate threshold on this? I am facing the same issue, it worked for me only when i choose 1e-4. On what basis we should choose that? I am using wgs84 lat long points – ds_user Aug 14 '17 at 03:48
  • I mean for someline, line.distance(point) returns true when i use <1e-4, and for some its 1e-3. Does that mean, I have to calculate the distance between the point and all line in my dataset and take the one with minimal distance? Is there an alternate efficient approach to this? – ds_user Aug 14 '17 at 03:59
  • @ds_user it really depends on how you obtained the data, or what resolution it represents (e.g. see [decimal degrees](https://en.wikipedia.org/wiki/Decimal_degrees)). You'd need to examine the unprojected data in a GIS to get an idea in your case. – Mike T Aug 14 '17 at 05:40
  • Actually I am trying to find the road name of an address. Two dataset I have, one with geopoints for address and another one with linestring for road. I need to get the road name based on which linestring is closest to my address point. I have a question on GIS stack exchange regarding this. – ds_user Aug 14 '17 at 05:43
  • https://gis.stackexchange.com/questions/251938/shapely-python-spatial-join-point-with-line – ds_user Aug 14 '17 at 05:50
  • 4
    This shouldn't be the case. `shapely` should handle the precision errors itself and not have users worry about this. python is supposed to be a high level language which abstracts machine precision from the user. Furthermore this renders the within and contain methods quite useless for linestrings and points in `shapely` – user32882 Sep 29 '18 at 08:18
  • Thanks for the wonderful solution. However, one issue that sprouts from here is: how to find whether, *e.g.* the `Point` is "on the road"-- *i.e.* within 6 m of the `LineString` which should denote the center-line of a 12 m wide road. – Partha D. Apr 20 '20 at 12:56
  • @ParthaD. sounds simple, just find `line.distance(point) < 6` – Mike T Apr 21 '20 at 03:16
  • @MikeT Seems my question was incomplete. My `point` and `line`s are in latitude, longitude, *e.g.* `line = LineString(-79.4176698 43.7005503, -79.4175783 43.7005733)` and `point = Point(-79.415847 43.70304155)` gives `line.distance(point)` as `0.0030149`. Later I found out that `from shapely.ops import nearest_points` and `shct = nearest_points(line, point)[0]` coupled with `1000*haversine(geo_dist(np.array(shct), np.array(loca))` using and the [haversine distance](https://stackoverflow.com/a/4913653/8854526) gives 199.00189`-- the required distance in metres :) – Partha D. Apr 21 '20 at 04:05
  • @ParthaD. that's a new question, but you generally can't reliably convert degrees to metres. – Mike T Apr 21 '20 at 08:36