What I am researching now is how to extend that to "range ellipses". A description here http://desktop.arcgis.com/en/arcmap/10.3/tools/data-management-toolbox/table-to-ellipse.htm. If we take this geodetic ellipse on the surface of the earth and generalize it to have a constant radius then we will have the range ring. The azimuth I think I can do by doing what I do for range circles but then I should use
For az = 0 To 360 Step N
azCalc = az + azTilt
Use azCalc and a length L to calculate end point
Next az
What I am not able to do is come up with an equation on how to compute the length L at an azCalc. Would it be based on simple ellipse formulas? or do geodesic calculations affect it? My instinct says yes the geodesics will come into play.
Any help or guidance is appreciated.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Whether or not it is sensible to generalize this to a spheroid is up to you. But certainly this definition reduces to a geodesic circle when a = b.
My own preference would be to define the ellipse in terms of its foci: the sum of the distances to the two foci is constant for each point on the ellipse. In general, this will give you a different curve compared to the polar representation. However you will recover the geodesic circle when the two foci are coincident.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Here is where I am with this challenge. I assume all my coordinates are in WGS84 and the input is ptStart and ptEnd. The I proceed as follows in .NET using the wrapper classes:
'assume a WGS84 datum
Dim geod As Geodesic = New Geodesic 'WGS84
Dim s12 As Double, azi1 As Double, azi2 As Double
Dim a12 As Double = geod.Inverse(ptStart.Y, ptStart.X, ptEnd.Y, ptEnd.X, s12, azi1, azi2)
Dim geoLine As GeodesicLine = New GeodesicLine(geod, ptStart.Y, ptStart.X, azi1, Mask.ALL)
'calculate geodesic distance half-way between the two points in meters. this will eventually be our semi-major axis.
Dim halfDistMeters As Double = 0.5 * s12
'compute the half-way point along the geodesic
Dim xLonMid As Double, yLatMid As Double
geoLine.Position(halfDistMeters, yLatMid, xLonMid)
Dim ptMid As Point = New Point(xLonMid, yLatMid)
'calculate the azimuth between the mid point to destination in degrees
Dim azMidToPoa As Double = GeodesicAzimuthInDegrees(ptMid, ptEnd)
'for great circle the flattening constant is zero. The flattening constant is f=(a-b)/a. if the semi-major axis, a, and the semi-minor axis, b,
'are equal (as they would be in a perfect sphere) then that becomes f=(a-a)/a=0/a=0
Dim a As Double = halfDistMeters
Dim f As Double = 0.6 'large to test the flattening of the ellipse. change to desired factor.
Dim b As Double = a * (1 - f)
'we will hard code the angle sweep to be from 0 to 360 by 2 degrees. in the future we can probably make this user configurable.
Dim radiusInMeters As Double, yLatPerim As Double, xLonPerim As Double
For az As Double = 0 To 360 Step 2
'radiusInMeters = GetDistAtTheta(az, a, b, azPodToPoa) 'this does not work
radiusInMeters = GetDistAtTheta(az, a, b, azMidToPoa) 'this does work
geod.Direct(ptMid.Y, ptMid.X, az, radiusInMeters, yLatPerim, xLonPerim)
Dim ptPerim As Point = New Point(xLonPerim, yLatPerim)
Next az
Public Function GeodesicAzimuthInDegrees(ByVal pt1 As Point, ByVal pt2 As Point) As Double
'assume a WGS84 datum
Dim geod As Geodesic = New Geodesic ' WGS84
Dim s12gd As Double, azi1gd As Double, azi2gd As Double
Dim a12gd As Double = geod.Inverse(pt1.Y, pt1.X, pt2.Y, pt2.X, s12gd, azi1gd, azi2gd)
If azi1gd < 0 Then azi1gd += 360 'adjust for negative angle and always make positive
Return azi1gd
End Function
Private Function GetDistAtTheta(ByVal azInDeg As Double, ByVal a As Double, ByVal b As Double, ByVal angRotationInDeg As Double) As Double
Dim thetaInDeg As Double = AzimuthToPolar(azInDeg) - AzimuthToPolar(angRotationInDeg)
'In polar coordinates, with the origin at the center of the ellipse and with the angular coordinate theta measured from the major axis,
'the ellipse's equation is: r(theta) = a*b / sqrt(b^2 * cos(theta)^2 + a^2 * sin(theta)^2)
Dim thetaInRad As Double = Math.PI * thetaInDeg / 180.0
Dim r As Double = (a * b) / Math.Sqrt(b ^ 2 * Math.Cos(thetaInRad) ^ 2 + a ^ 2 * Math.Sin(thetaInRad) ^ 2)
Return r
End Function
'https://math.stackexchange.com/questions/926226/conversion-from-azimuth-to-counterclockwise-angle
Private Function AzimuthToPolar(ByVal azInDeg As Double) As Double
'For azimuth a, compute h=450−a; if h>360, return h−360; otherwise return h.
Dim h As Double = 450 - azInDeg
If h > 360 Then
Return h - 360
Else
Return h
End If
End Function
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
OK, this looks fine. However, I remind you that you've only implemented one possible definition of an "ellipse" on a general surface. Which is the best definition to use will depend on your application.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Completely understand. I took the specific problem of a range circle and tried to generalize it to a range ellipse. I tried to incorporate geodesics to come up with a suitable definition. For my purposes I am trying to minimize the points used in a path finding algorithm. The flatness of the ellipse will determine how "far out" we will search for connection points. Many thanks for your help.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Based on the discussion on range rings from this post https://sourceforge.net/p/geographiclib/discussion/1026621/thread/1df0a86b/ I have been able to recreate geodesic circles that are dead-on when compared to ArcGIS.
What I am researching now is how to extend that to "range ellipses". A description here http://desktop.arcgis.com/en/arcmap/10.3/tools/data-management-toolbox/table-to-ellipse.htm. If we take this geodetic ellipse on the surface of the earth and generalize it to have a constant radius then we will have the range ring. The azimuth I think I can do by doing what I do for range circles but then I should use
For az = 0 To 360 Step N
azCalc = az + azTilt
Use azCalc and a length L to calculate end point
Next az
What I am not able to do is come up with an equation on how to compute the length L at an azCalc. Would it be based on simple ellipse formulas? or do geodesic calculations affect it? My instinct says yes the geodesics will come into play.
Any help or guidance is appreciated.
For a plane ellipse, the polar representation relative to its center is
see https://en.wikipedia.org/wiki/Ellipse#Polar_form_relative_to_center
Whether or not it is sensible to generalize this to a spheroid is up to you. But certainly this definition reduces to a geodesic circle when a = b.
My own preference would be to define the ellipse in terms of its foci: the sum of the distances to the two foci is constant for each point on the ellipse. In general, this will give you a different curve compared to the polar representation. However you will recover the geodesic circle when the two foci are coincident.
OK let me give this some thought and i will see what I can do with it.
Here is where I am with this challenge. I assume all my coordinates are in WGS84 and the input is ptStart and ptEnd. The I proceed as follows in .NET using the wrapper classes:
Here are a couple of samples:
MIA to ORD
BNA to CUN
LAX to YYZ
Last edit: AGP 2019-08-09
OK, this looks fine. However, I remind you that you've only implemented one possible definition of an "ellipse" on a general surface. Which is the best definition to use will depend on your application.
Completely understand. I took the specific problem of a range circle and tried to generalize it to a range ellipse. I tried to incorporate geodesics to come up with a suitable definition. For my purposes I am trying to minimize the points used in a path finding algorithm. The flatness of the ellipse will determine how "far out" we will search for connection points. Many thanks for your help.