pascal/lazarus-ccr
1
0
Fork 0
Code Issues Releases Activity

LazMapViewer: Add unit test for distance calculation

Browse Source 
git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@8793 8e941d3f-bd1b-0410-a28a-d453659cc2b4
This commit is contained in:
wp_xxyyzz
2023-04-18 17:48:01 +00:00
parent fd7b0d8af9
commit 649286bb6b
2 changed files with 67 additions and 31 deletions
Download Patch File Download Diff File Expand all files Collapse all files

52
components/lazmapviewer/source/mvengine.pas
View File

@@ -166,6 +166,7 @@ type
function RealPoint(Lat, Lon: Double): TRealPoint;
function HaversineDist(Lat1, Lon1, Lat2, Lon2, Radius: Double): Double;
function CalcGeoDistance(Lat1, Lon1, Lat2, Lon2: double;
AUnits: TDistanceUnits = duKilometers): double;
@@ -734,7 +735,7 @@ const
MIN_LONGITUDE = -180;
MAX_LONGITUDE = 180;
var
zoomfac: Extended;
zoomfac: Int64;
begin
// https://epsg.io/3857
// https://pubs.usgs.gov/pp/1395/report.pdf, page 41
@@ -776,7 +777,7 @@ var
Cpm: Extended;
Z: Integer;
t, phi: Extended;
zoomFac: Extended; // 2**Z
zoomFac: Int64;
i: Integer;
begin
// https://epsg.io/3395
@@ -1549,42 +1550,37 @@ begin
ADeg := sgn * (abs(ADeg) + mins / 60 + secs / 3600);
end;
// https://stackoverflow.com/questions/73608975/pascal-delphi-11-formula-for-distance-in-meters-between-two-decimal-gps-point
function HaversineDist(Lat1, Lon1, Lat2, Lon2, Radius: Double): Double;
var
latFrom, latTo, lonDiff: Double;
dx, dy, dz: Double;
begin
lonDiff := DegToRad(Lon1 - Lon2);
latFrom := DegToRad(Lat1);
latTo := DegToRad(Lat2);
dz := sin(latFrom) - sin(latTo);
dx := cos(lonDiff) * cos(latFrom) - cos(latTo);
dy := sin(lonDiff) * cos(latFrom);
Result := arcsin(sqrt(sqr(dx) + sqr(dy) + sqr(dz)) / 2) * Radius * 2;
end;
{ Returns the direct distance (air-line) between two geo coordinates
If latitude NOT between -90°..+90° and longitude NOT between -180°..+180°
the function returns -1.
Usage: FindDistance(51.53323, -2.90130, 51.29442, -2.27275, duKilometers);
the function returns NaN.
Usage: CalcGeoDistance(51.53323, -2.90130, 51.29442, -2.27275, duKilometers);
}
function CalcGeoDistance(Lat1, Lon1, Lat2, Lon2: double;
AUnits: TDistanceUnits = duKilometers): double;
const
EPS = 1E-12;
var
d_radians: double; // distance in radians
lat1r, lon1r, lat2r, lon2r: double;
arg: Double;
begin
// Validate
if (Lat1 < -90.0) or (Lat1 > 90.0) then exit(NaN);
// if (Lon1 < -180.0) or (Lon1 > 180.0) then exit(NaN);
if (Lat2 < -90.0) or (Lat2 > 90.0) then exit(NaN);
// if (Lon2 < -180.0) or (Lon2 > 180.0) then exit(NaN);
// Turn lat and lon into radian measures
lat1r := (PI / 180.0) * Lat1;
lon1r := (PI / 180.0) * Lon1;
lat2r := (PI / 180.0) * Lat2;
lon2r := (PI / 180.0) * Lon2;
// calc
arg := sin(lat1r) * sin(lat2r) + cos(lat1r) * cos(lat2r) * cos(lon1r - lon2r);
if (arg < -1) or (arg > +1) then
exit(NaN);
if SameValue(abs(Lon1-Lon2), 360, EPS) and SameValue(abs(arg), 1.0, EPS) then
d_radians := PI * 2.0
else
d_radians := arccos(arg);
Result := EARTH_EQUATORIAL_RADIUS * d_radians;
Result := HaversineDist(Lat1, Lon1, Lat2, Lon2, EARTH_EQUATORIAL_RADIUS);
case AUnits of
duMeters: ;
duKilometers: Result := Result * 1E-3;