0

Since I'm not a mathematician I came here to ask what the most efficient way is to check if latitude and longitude coordinates are inside a range (for example 50 meters) of multiple latitude and longitude points (polygon).

I have a list of these longitude and latitude points:

[6.38537265,51.87721088],[6.38542453,51.87737201],[6.38523252,51.87739419],[6.38477205,51.87745015],[6.38426164,51.87751088],[6.38391099,51.87755068],[6.38386033,51.87738808],[6.38380232,51.87720004],[6.38376297,51.87708017],[6.38375183,51.87704018],[6.38373055,51.8769829],[6.38390723,51.87695904],[6.38389144,51.87691388],[6.38403002,51.87690001],[6.38440124,51.8768538],[6.38493939,51.87678787],[6.38522535,51.87675316],[6.38529885,51.87697928],[6.38537265,51.87721088]

And would like to check if these coordinates are in the range:

51.877368,6.383818

Here's a sketch of my question to visualize what I mean. Arrows represent a given latitude and longitude combination.

Click here for the image of my sketch

I was planning to use the haversine formula but I do not know how I could use that in a list of multiple points.

Thanks in advance!

SaltyPotato
  • 3
  • So does "in range of a point" mean that the distance to a point is smaller than some fixed value? I don't see what this has to do with polygons ... – Matti P. Dec 05 '18 at 14:08
  • @MattiP. Yes, well I'm trying to check if the given latitude and longitude is within the fixed value (50 meters) of the polygon (made of multiple latitude and longitude coordinates). – SaltyPotato Dec 05 '18 at 14:14
  • I want to ask you to give an answer to a sample problem so that I can know what you're trying to do. Suppose the input polygon is a triangle, all three vertices are on the equator, spaced 120 degrees apart. Is the north pole within this polygon? Is the south pole within it? – John Hughes Dec 05 '18 at 14:20
  • I'm wondering whether you care about the interior of the polygon at all. Perhaps your polygon is more a "poly-line", a sequence of segments that doesn't necessarily close up, and you want to know whether some test-point is within a distance $d$ of any point on this polyline. Is that in fact your question? I can imagine this being useful in navigational software ("Did the boat ever pass within 50 meters of mark "RB2"?"). – John Hughes Dec 05 '18 at 14:21
  • By the way, you're likely to get several more comments trying to clarify your question -- that may be the toughest part of the process. I can say with some confidence that if your real question is any of the several that I can imagine, then it won't be too hard to solve it. But I'm not willing to solve every possible question that might be the one you're asking (and nor are others), hence these comments asking for clarifications. – John Hughes Dec 05 '18 at 14:24
  • To follow up on my "poly-line" question, do you care about the test point being within distance $d$ of the whole poly-line (i.e., including edges) or just within 50 meters of any vertex of the polyline? I.e., are we doing a point-to-finite-set-of-points distance computation, or point-to-sequence-of-segments computation? – John Hughes Dec 05 '18 at 14:26
  • @JohnHughes The given latitude and longitude points always close up, they are data retrieved from a cadastral database and represent a building. The use case will be to check if a user is in the range of about 50 metres or inside the "polygon" of the building. So your use case of "("Did the boat ever pass within 50 meters of mark "RB2"?")" is kind of in the right direction I'm aiming for. – SaltyPotato Dec 05 '18 at 14:31
  • @JohnHughes The point must be within the distance d of the whole "polygon" including the edges. So a point-to-finite-set-of-points distance computation. – SaltyPotato Dec 05 '18 at 14:34
  • I think you've misunderstood me. What if the building is a square, 500 feet on each side? Checking the distance from the user to the 4 corners (dist-from-point-to-finite-set-of-points) won't do what you want. Checking distance to edges won't either (user could be in the middle of the square, 250 feet from each edge!). – John Hughes Dec 05 '18 at 15:05
  • Let me try to write your question in a way I think captures everything you've said: Given the sequence $p_1, p_2, \ldots, p_n$ of vertices of a spherical polygon on the unit sphere, where the "interior" of the polygon is to our left if we traverse the vertices in order, and where no two vertices are more than, say, 90 degrees apart on the sphere, and given a test point $q$ on the sphere, determine whether $q$ is either inside the polygon defined by the points $p_1, \ldots$, or whether the distance from $q$ to any of the edges of the polygon is less than or equal to some fixed distance $d$. – John Hughes Dec 05 '18 at 15:09
  • @JohnHughes I think we're on thinking on the same line now. Could you please check my initial question again? I've updated it with an image to visualize what I mean, to make sure we are both talking about the same thing. – SaltyPotato Dec 05 '18 at 15:24
  • Shouldn't the upper left corner of your red shape be replaced by a circle-arc, consisting of points that are 50m from the corner of the building? (And similarly for 4 other corners, but not the lower right "indented" corner). – John Hughes Dec 05 '18 at 15:38
  • @JohnHughes Yes, however I couldn't find that functionality in the sketch software I used haha. But you're correct. – SaltyPotato Dec 05 '18 at 15:43
  • Since the range of your points is so small, there's no point in using the haversine formula unless you need sub-millimeter precision. Just put the points on a plane. – Connor Harris Dec 05 '18 at 17:33

1 Answers1

0

This is the start of a general answer, but not yet complete. Still, I'm going to save it and continue work later, in part because I believe that the careful reformulation of the problem may well be the largest contribution I can make.

Presumably you're doing this in software. So I'm gonna write pseudocode in a matlab-like language. It's not pretty, but it'll do. 

Inputs: 
  r: The radius of the spherical approximation of earth, in some units
  d: the distance that the test point needs to be within to produce "true"; same units. 
  p: an n x 2 array of points on the sphere, in lat/lon coordinates, with 
     90 degrees being the north pole, -90 being the south pole, etc. 
  q: a 1 x 2 array consisting of the lat and lon of the "user" or "test point"

Assumptions: 
  1. The earth is spherical enough, and the terrain flat enough, that lat-lon 
  distance between points is a good approximation of absolute distance. A 
  failure-case for this would be at half-dome in Yosemite, where two people could be 
  100 yards apart when radially  projected onto the spherical approximation of the
  earth, but many hundreds of yards apart in reality, one at the top of half-
  dome, the other standing on the ground below. 

  2. The polygon defined by the points in the array p is defined so that 
  traversing the polygon edges in the order p1, p2, p3, ..., the interior
  of the polygon is always to our left. There are no duplicated points,
  and the polygon forms a simple closed curve on the surface (i.e., there are
  no self-intersections). 

  3. No two points of the polygon are more than 90 degrees distant on the 
  sphere (I'm not certain I'll use this, but just in case...); the entire 
  polygon is contained within some hemisphere (also not sure i need this,
  but...)

  4. The test-point is within a distance of 90 degrees of all the polygon 
  vertices. 

  5. If the input points are written in degrees, then the cosine and sine and 
  other trig functions work with degrees; if the inputs are in radians, then 
  trig functions work with radians. 

Output: TRUE if either one of the following two conditions holds:
    A. the test point q is within the polygon defined by the points p OR
    B. q is within spherical distance d of one of the polygon edges
  and FALSE otherwise.

// ---- A few helper functions ----

// convert lat-lon coordinates into an xyz triple representing the 
// corresponding point on the unit sphere. 
function [x, y, z] = rect_from_polar(lat, lon)
   y = sin(lat); 
   x = cos(lat) * cos(lon); 
   z = cos(lat) * sin(lon); 

function [lat, lon] = polar_from_rect(x, y, z)
  lat = arcsin(y); 
  if (y == 1) or (y == -1) {
     lon = 0; 
  }
  else {
     lon = atan2(z, x);
  }
John Hughes
  • 93,729