From 7551440ac06282574e9164154d8e6457b1c65fdb Mon Sep 17 00:00:00 2001 From: Erik Shagdar Date: Mon, 7 Jan 2019 10:10:02 -0500 Subject: [PATCH] added existing geo functions to repo. Updated GeoDistance to return 0 if the specified points are the same. --- CustomFunctions/GeoBoxBounds.fmfn | 104 ++++++++++++++++++++++++++++++ CustomFunctions/GeoDistance.fmfn | 72 +++++++++++++++++++++ 2 files changed, 176 insertions(+) create mode 100644 CustomFunctions/GeoBoxBounds.fmfn create mode 100644 CustomFunctions/GeoDistance.fmfn diff --git a/CustomFunctions/GeoBoxBounds.fmfn b/CustomFunctions/GeoBoxBounds.fmfn new file mode 100644 index 0000000..bfffeee --- /dev/null +++ b/CustomFunctions/GeoBoxBounds.fmfn @@ -0,0 +1,104 @@ +// GeoBoxBounds ( latitude ; longitude ; radiusInMeters ) +// version 2012-09-27, Jeremy Bante, http://scr.im/jbante + +/** + * ===================================== + * GeoBoxBounds ( latitude ; longitude ; radiusInMeters ) + * + * PURPOSE: + * Calculates the approximate coordinate bounds for an approximately square region + * around a point on the surface of the Earth. This can be useful for finding data + * within a certain distance of a reference point. + * + * RETURNS: + * A return-delimited list of 4 numbers: + * 1. The south latitude boundary for the calculated region + * 2. The north latitude boundary for the calculated region + * 3. The west longitude boundary for the calculated region + * 4. The east longitude boundary for the calculated region + * + * PARAMETERS: + * latitude: Degrees between -90 and 90 for center of geobox + * longitude: Degrees between -180 and 180 for center of geobox + * radiusInMeters: The distance in meters from the center of the region to the + * closest points on the edges of the square bounding region. + * + * DEPENDENCIES: none + * + * EXAMPLE: Find records within a search radius + * Set Variable [$geoBox; Value:GeoBoxBounds ( $latitude ; $longitude ; $radius)] + * Enter Find Mode [] + * Set Field [Table::latitude; Value:GetValue ( $geoBox ; 1 ) & "..." & GetValue ( $geoBox ; 2 )] + * Set Field [Table::longitude; Value:GetValue ( $geoBox ; 3 ) & "..." & GetValue ( $geoBox ; 4 )] + * Perform Find [] + * + * NOTES: + * Use FileMaker's GetValue function to retrieve each value from the result. + * + * When the geobox includes a pole, the function will return the appropriate values + * to perform a find for coordinates enclosed by the lower latitude. For example: + * [90, 88, -180, 180] for a region around the north pole. + * + * When the geobox crosses the antimeridian (180º/-180º longitude), the function will + * return a west boundary coordinate greater than the east boundary coordinate. In + * this situation, a coordinate search within the geobox should be broken into + * separate searches: one search from the west bound to 180º, and another from -180º + * to the east bound. + * + * The bounding region's approximation to a square shape is less accurate as + * radiusInMeters grows large and as the region approaches either pole. However, the + * region will always contain the circular region of the same radius. + * + * Square bounding regions are more computationally convenient for finding data. + * Developers needing an exact circular region can use a square region to constrain + * the available data to a smaller found set, then use the GeoDistance function to + * omit points between the circular and square neighborhoods of the reference point. + * + * HISTORY: + * MODIFIED on 2012-09-27 by Jeremy Bante to account for + * geoboxes that may span a pole or meridian. + * CREATED on 2012-09-26 by Jeremy Bante . + * + * REFERENCES: + * Radius of Earth: http://www.wolframalpha.com/input/?i=radius+of+earth + * ===================================== + */ + +Let ( [ + ~earthRadius = 6367500; // meters + + // calculate latitude bounds + ~latitudeOffset = Degrees ( radiusInMeters / ~earthRadius ); + ~north = latitude + ~latitudeOffset; + ~south = latitude - ~latitudeOffset; + ~wrapNorth = ~north > 90; + ~wrapSouth = ~south < -90; + ~northBound = + Case ( + ~wrapNorth ; 90; + ~wrapSouth ; Max ( ~north ; -180 - ~south ); + ~north + ); + ~southBound = + Case ( + ~wrapSouth ; -90; + ~wrapNorth ; Min ( ~south ; 180 - ~north ); + ~south + ); + + // calculate longitude bounds + ~longitudeOffset = ~latitudeOffset * Cos ( Radians ( latitude ) ); + ~eastBound = + If ( ~wrapNorth or ~wrapSouth ; 180 ; /* Else */ longitude + ~longitudeOffset ); + ~eastBound = If ( ~eastBound > 180 ; ~eastBound - 360 ; /* Else */ ~eastBound ); + ~westBound = + If ( ~wrapNorth or ~wrapSouth ; -180 ; /* Else */ longitude - ~longitudeOffset ); + ~westBound = If ( ~westBound < -180 ; ~westBound + 360 ; /* Else */ ~westBound ) +]; + List ( + ~southBound; + ~northBound; + ~westBound; + ~eastBound + ) +) diff --git a/CustomFunctions/GeoDistance.fmfn b/CustomFunctions/GeoDistance.fmfn new file mode 100644 index 0000000..eef92ef --- /dev/null +++ b/CustomFunctions/GeoDistance.fmfn @@ -0,0 +1,72 @@ +// GeoDistance ( latitude1 ; longitude1 ; latitude2 ; longitude2 ) +// version 2019-01-04, Jeremy Bante, http://scr.im/jbante + +/******************************************************************************* + * LocationDistance ( latitude1 ; longitude1 ; latitude2 ; longitude2 ) + * Calculates the distance between two points on the surface of the Earth. This + * function assumes that coordinates are in the WGS84 coordinate system. + * + * @parameter latitude1: Latitude of the first point + * @parameter longitude1: Longitude of the first point + * @parameter latitude2: Latitude of the second point + * @parameter longitude2: Longitude of the second point + * + * @return The distance between the two points, in meters. + * + * @history 2012-09-26 - Jeremy Bante - Created + * @history 2014-09-15 - Jeremy Bante - Using the + * Andoyer-Lambert-Thomas formula approximating distances on an ellipsoid + * instead of a sphere to improve accuracy. + * @history 2017-03-16 - Erik Shagdar - updated function to use the + * specified parameters. + * @history 2019-01-04 - Erik Shagdar - Return 0 if the same point is used. + * + * @see http://www.dtic.mil/dtic/tr/fulltext/u2/703541.pdf + * @see http://en.wikipedia.org/wiki/World_Geodetic_System + ******************************************************************************/ + +If ( latitude1 = latitude2 and longitude1 = longitude2 ; + 0 ; + Let ( [ + // WGS84 + _equatorialRadius = 6378137 ; // meters + _flattening = .0033528106647475 ; // 1 / 298.257223563 + + _θ1 = + Atan ( + ( 1 - _flattening ) * Tan ( Radians ( latitude1 ) ) + ) ; + _θ2 = + Atan ( ( 1 - _flattening ) * Tan ( Radians ( latitude2 ) ) ) ; + _θm = ( _θ1 + _θ2 ) / 2 ; + _sinθm = Sin ( _θm ) ; + _cosθm = Cos ( _θm ) ; + _Δθm = ( _θ2 - _θ1 ) / 2 ; + _sinΔθm = Sin ( _Δθm ) ; + _cosΔθm = Cos ( _Δθm ) ; + _Δλ = Radians ( longitude2 - longitude1 ) ; + _h = _cosΔθm ^ 2 - _sinθm ^ 2 ; + _l = _sinΔθm ^ 2 + _h * Sin ( _Δλ / 2 ) ^ 2 ; + _cos_d = 1 - 2 * _l ; + _d = Acos ( _cos_d ) ; + _sin_d = Sin ( _d ) ; + _uu = 2 * _sinθm ^ 2 * _cosΔθm ^ 2 / ( 1 - _l ) ; + _vv = 2 * _sinΔθm ^ 2 * _cosθm ^ 2 / _l ; + _x = _uu + _vv ; + _t = _d / _sin_d ; + _e = 2 * _cos_d ; + _y = _uu - _vv ; + _dd = 4 * _t ^ 2 ; + _b = 2 * _dd ; + _a = _dd * _e ; + _c = _t - .5 * ( _a - _e ) ; + _n1 = _x * ( _a + _c * _x ) ; + _n2 = _y * ( _b + _e * _y ) ; + _n3 = _dd * _x * _y ; + _δ1d = .25 * _flattening * ( _t * _x - _y ) ; + _δ2d = ( _n1 - _n2 + _n3 ) * _flattening ^ 2 / 64 ; + _distance = _equatorialRadius * _sin_d * ( _t - _δ1d + _δ2d ) + ] ; + _distance + ) +)