Detecting lines which pass close to points, lines or polygons

[mikemcc]

I'm looking for software to analyze a large number of polylines and
detect whether each one passes:
- either within a radius of x miles of a point which I specify
- or within x miles of another polyline which I specify
- or through a polygon which I specify.

I then need to display just the relevant polylines in the GE plugin
(because that's what my customers can use now).

I have all of the polylines, points and polygons defined in Google
Earth kml files. I also have the points defined in a MySql table and
have other MySql tables which refer to all the above kml files.

I have written various php scripts to open, analyze and manipulate all
this info in other ways and display results in Google Earth, so my
first instinct is to try to find the necessary geometric functions and
write the scripts myself.

However, I haven't found any yet (inc asking in the server-side forum)
and I'd happily buy some "industry-standard" software if necessary. As
I'm not a mathemetician and as GE is the only GIS software I have used
so far, I'm not sure:
- whether the above analysis is geometrically possible (I imagine it
is)
- what GIS software can do and how well it interfaces with GE.

So my question is: are there any software products which can allow me
to do the above analysis so I can present the results in Google
Earth ?

Thanks. Mike.

[Frankv]

Hi Mike,

I've done the first two in Visual Basic .Net.

I found all the information using Google searches...

Here's code for distance from a point to a line segment... repeat this
for each segment in the polyline. Note that you can't do this in lat/
lon -- you have to either do 3D calcs (as I do... assuming a spherical
Earth) or project to a 2D representation to calculate distance.

Class xyz
Public x, y, z As Double
Public Sub New(ByVal x As Double, ByVal y As Double, ByVal z As
Double)
Me.x = x
Me.y = y
Me.z = z
End Sub
Public Sub New(ByVal P As coord)
' P needs to be lat/long in degrees!
Me.New(toRad(P.longitude), toRad(P.latitude))
End Sub
Public Sub New(ByVal longitudeRad As Double, ByVal latitudeRad As
Double)
' P needs to be in Radians!
Dim theta As Double = longitudeRad
' spherical coordinate use "co-latitude", not "latitude"
' lattiude = [-90, 90] with 0 at equator
' co-latitude = [0, 180] with 0 at north pole
Dim phi As Double = Math.PI / 2.0 - latitudeRad
Me.x = Math.Cos(theta) * Math.Sin(phi)
Me.y = Math.Sin(theta) * Math.Sin(phi)
Me.z = Math.Cos(phi)
End Sub
End Class

Private Function sqrDistPP3D(ByVal v1 As xyz, ByVal v2 As xyz) As
Double
Return (v1.x - v2.x) * (v1.x - v2.x) + (v1.y - v2.y) * (v1.y -
v2.y) + (v1.z - v2.z) * (v1.z - v2.z)
End Function

Public Function distance_lineseg(ByVal A As coord, ByVal B As
coord, Optional ByRef near_pt As coord = Nothing) As Double
If A.Equals(B) Then
near_pt = A
Return Me.distance(A)
End If
If Me.Equals(A) Then
near_pt = A
Return 0
End If
If Me.Equals(B) Then
near_pt = B
Return 0
End If

' Convert form lat/lon to Cartesian coordinates.
Dim P1 As xyz = New xyz(Me)
Dim A1 As xyz = New xyz(A)
Dim B1 As xyz = New xyz(B)

Dim sqrP_LB As Double = sqrDistPP3D(P1, A1)
Dim sqrP_LE As Double = sqrDistPP3D(P1, B1)
Dim sqrLB_LE As Double = sqrDistPP3D(A1, B1)
Dim LB_LE As Double = Math.Sqrt(sqrLB_LE)
Dim I_LB As Double = (sqrP_LB + sqrLB_LE - sqrP_LE) / (2 *
LB_LE)
If (I_LB <= 0) Then
near_pt = A ' intersection point is before beginning line
segment
ElseIf (I_LB >= LB_LE) Then
near_pt = B ' intersection point is after end line segment
Else

' use this part if you want to know what the closest
point on the line is.
Dim u As Double = I_LB / LB_LE
near_pt = New coord(A1.x + u * (B1.x - A1.x), A1.y + u *
(B1.y - A1.y), A1.z + u * (B1.z - A1.z))
End If
Return Me.distance(near_pt)
End Function

How it works...
Call the end points of the segment A & B, and the given point P
Imagine a circle centred on P and tangential to the line (not
segment) extending through A & B. The radius of this circle is the
distance from the line to P.
If the intersection point I is between A & B, then the shortest
distance is the radius of the circle. Otherwise, the shortest distance
is the minimum of the lengths AP and BP.

So you need to find the intersection point I. You're looking for the
intersection of the line passing through AB and the line IP (starting
at P, direction at right angles to AB).

If you can understand the above, get back to me and I'll sort out the
distance-between-polylines code that I have somewhere.

Identifying intersecting polylines and polygons is relatively simple
-- it breaks down to checking each polyline segment in turn, to see
whether it intersects the polygon. You also need to handle the special
case where both ends are inside the polygon, but there's a standard
algorithm for identifying if a point is inside a polygon or not. But
this can get very complex, especially if you want the polygon to be a
complex shape (e.g. has an InnerRing and/or is an irregular shape and/
or segments that cross each other).

If you want to pay someone to write this code, I'm happy to do it. :)
I can write PHP, no problem.

Frank

[mikemcc]

Wow, thanks very much, Frank. Please give me time to absorb this.
It'll involve understanding the VB syntax too.

Meanwhile I had:
1. Worked out some of the maths terminology eg circle line
intersection
2. Googled that terminology to find a formula at
http://mathworld.wolfram.com/Circle-LineIntersection.html
3. Written most of that formula in php; I am confused by the +/- stuff
4. Prepared some questions for my daughter who promises to try to help
me shortly:

In the equations for the intersections x and y, what is the syntax
immediately after the plus/minus symbols meant to say ?

Why aren't the points of intersection defined as 2 coordinates each ?

What is the unit of length ? Nautical miles ?

Is this calculation valid at any latitude (where the lines of
longitude get closer together) ? I note (but don't yet understand)
what you say about 2D and 3D. In my PHP code I have reduced the coords
of the line to be relative to the centre of the circle. These calcs
don't have to be super-accurate - as long as I know how inaccurate
they are for eg a 20 mile radius at differnt parts of the world.

Is it the case that working out intersections x and y is irrelevant to
determining whether the line intersects ?

Other scenarios
---------------
What if the line is finite ? ie I'm analyzing using straight, finite
segments of polylines
What if the finite line ends within the circle ?
What if the finite line starts and ends within the circle ?

Next calculation
----------------
If there are two intersections or the line starts and/or ends inside
the circle, calculate the length of the line within the circle.

Once she educates me a bit about the above, I imagine I'll be back to
you, knowing more about what I need than I do now.

PS what time zone are you in ? I'm in the UK.

Thanks again. Mike.

[Frankv]

> 2. Googled that terminology to find a formula at http://mathworld.wolfram.com/Circle-LineIntersection.html

I recall that page... probably I used it in developing my code.

> 3. Written most of that formula in php; I am confused by the +/- stuff

There are always two results of a square root, with the same
magnitude, but positive or negative. (e.g. SQRT(4) can be either +2 or
-2). In geometry of circles, these typically correspond to the
intersections on opposite sides of the circle.

> In the equations for the intersections x and y, what is the syntax
> immediately after the plus/minus symbols meant to say ?

The function sgn*() is described further down the page.

The vertical bars mean "absolute value"

> Why aren't the points of intersection defined as 2 coordinates each ?

If you work with a tangent, then there is a single point of
intersection.

> What is the unit of length ? Nautical miles ?

I work in metres. But to some extent it's irrelevant... the important
thing is to use the same units throughout. The answer will be in terms
of whatever units you you for the various distances.

> Is this calculation valid at any latitude (where the lines of
> longitude get closer together) ? I note (but don't yet understand)
> what you say about 2D and 3D. In my PHP code I have reduced the coords
> of the line to be relative to the centre of the circle. These calcs
> don't have to be super-accurate - as long as I know how inaccurate
> they are for eg a 20 mile radius at differnt parts of the world.

You can't work in lat/lon except at the equator. I convert from lat/
lon to XYZ... points relative to the centre of the Earth, in metres.

> Is it the case that working out intersections x and y is irrelevant to
> determining whether the line intersects ?

x & y aren't intersections... they are components of one intersection.

Frank

[Frankv]

I ran out of time writing the previous post.

A Google for 'distance point line' (without the quotes) brings up lots
of hits

http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html is
quite technical.
http://en.wikipedia.org/wiki/Tangent_lines_to_circles has a good
description of circles and tangents.
http://www.worsleyschool.net/science/files/linepoint/distance.html has
5 different methods.... wish I'd found this before I wrote my
code :rolleyes:

Method 5 is kind of a cookbook approach, and you would think it was
easiest. But this solution is for a line, NOT a line segment. So in
fact Method 4 is better.

The following doesn't match my code in my previous post, and I haven't
checked it at all. So take it all with a grain of salt, and test
thoroughly.

Following Method 4, then:
Step 1:
You already have 2 points A (x1, y1) and B (x2, y2) defining a line
segment so you don't need to do step 1.

Step 2:
By Pythagoras, the distance between two points is sqrt( (x1-x2)^2 +
(y1-y2)^2 )

But this applies to planar, Cartesian geometry, not latitude/longitude
coordinates. So here's a simple, ROUGH, way to calculate distances
between points specified as lat/lon -- it gives answers in nautical
miles.

This assumes a locally flat Earth... it takes straight lines, rather
than following the curve of the Earth, so it will always underestimate
a little. The further apart the points, the greater the error. But it
should work OK for short distances, say less than 300nm or so. It also
assumes a spherical Earth, so I guess it will overestimate the North/
South component of the distance.

1 degree of latitude = 60nm (near enough, assuming a spherical Earth)
As you said, the length of a degree of longitude varies depending on
the latitude... at the equator, it is 60nm (by definition); at the
poles, it is 0
So 1 degree of longitude = 60 * (90-abs(lat))/90 nm

So Pythagoras's formula for distances on the Earth ROUGHLY =
sqrt( ((x1-x2)*60*((90-abs(y1+y2)/2)/90))^2 + ((y1-y2)*60)^2 )

where (x1,y1) is the long/lat of one point, (x2, y2) is the other.
Note that I take the average of the two latitudes (y1+y2)/2 to
calculate the length of a degree of longitude.

I use abs() to indicate absolute value, sqrt() to indicate square
root, and ^2 to indicate squared. I'm sure PHP has a function for the
first two, probably a power() function or similar too???

You'll need to use this numerous times to calculate distances between
several points, so write it as a function in PHP.

Step 3: as per the web page, but note that this assumes that the
intercept lies between A & B -- you should check for an oblique
triangle by comparing s with a, b, and c. If s is less than any of
these, then the triangle is oblique, and the calculation s * (s- a) *
(s - b) * (s - c) would result in a negative number, which you then
couldn't calculate the square root of. If s = a or b or c then the
calculation will be 0, which also won't let you get the right answer.
But all is not lost... in either of these cases, the nearest point to
P will be one of the ends of the segment (either A or B)... simply
compare the distance AP against BP... the minimum of these will be the
distance from P to the line segment.

There's also a special case of A,B, & P being colinear, in this case
the calculation of s will also be 0. You need to check whether P lies
between A & B -- if it's between (just compare (say) the x coordinate
of P to that of A & B), then distance is 0, otherwise the distance is
the minimum of the distances AP and BP.

Step 4:
As per the web site.

[mikemcc]

Thanks very much Frankv and I apologize for not spotting and
responding to your major contribution until just now. I must absorb it
as and when I have a moment - I'm trying to develop the application
while doing all my normal work. Much appreciated. Mike.