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

mapviewer: Replace repeated calls to IntPower(2,..) by table lookup function ZoomLevel. Add elemental unit tests.

Browse Source 
git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@8790 8e941d3f-bd1b-0410-a28a-d453659cc2b4
This commit is contained in:
wp_xxyyzz
2023-04-17 15:42:06 +00:00
parent a530b02b06
commit 7a09190330
6 changed files with 374 additions and 26 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -171,27 +171,53 @@ function CalcGeoDistance(Lat1, Lon1, Lat2, Lon2: double;
function DMSToDeg(Deg, Min: Word; Sec: Double): Double;
function GPSToDMS(Angle: Double): string;
function GPSToDMS(Angle: Double; AFormatSettings: TFormatSettings): string;
function LatToStr(ALatitude: Double; DMS: Boolean): String;
function LatToStr(ALatitude: Double; DMS: Boolean; AFormatSettings: TFormatSettings): String;
function LonToStr(ALongitude: Double; DMS: Boolean): String;
function LonToStr(ALongitude: Double; DMS: Boolean; AFormatSettings: TFormatSettings): String;
function TryStrToGps(const AValue: String; out ADeg: Double): Boolean;
procedure SplitGps(AValue: Double; out ADegs, AMins, ASecs: Double);
function ZoomFactor(AZoomLevel: Integer): Int64;
var
HERE_AppID: String = '';
HERE_AppCode: String = '';
OpenWeatherMap_ApiKey: String = '';
ThunderForest_ApiKey: String = '';
DMS_Decimals: Integer = 1;
implementation
uses
Forms, laz2_xmlread, laz2_xmlwrite, laz2_dom, TypInfo,
Forms, TypInfo, laz2_xmlread, laz2_xmlwrite, laz2_dom,
mvJobs, mvGpsObj;
const
_K = 1024;
_M = _K*_K;
_G = _K*_M;
ZOOM_FACTOR: array[0..32] of Int64 = (
1, 2, 4, 8, 16, 32, 64, 128, 256, 512, // 0..9
_K, 2*_K, 4*_K, 8*_K, 16*_K, 32*_K, 64*_K, 128*_K, 256*_K, 512*_K, // 10..19
_M, 2*_M, 4*_M, 8*_M, 16*_M, 32*_M, 64*_M, 128*_M, 256*_M, 512*_M, // 20..29
_G, 2*_G, 4*_G // 31..32 // 30..32
);
function ZoomFactor(AZoomLevel: Integer): Int64;
begin
if (AZoomLevel >= Low(AZoomLevel)) and (AZoomLevel <= High(AZoomLevel)) then
Result := ZOOM_FACTOR[AZoomLevel]
else
Result := round(IntPower(2, AZoomLevel));
end;
type
{ TLaunchDownloadJob }
@@ -614,8 +640,9 @@ const
MAX_LATITUDE = 85.05112878;
MIN_LONGITUDE = -180;
MAX_LONGITUDE = 180;
TWO_PI = 2.0 * pi;
var
px, py: Extended;
factor, px, py: Extended;
pt: TRealPoint;
begin
// https://epsg.io/3857
@@ -624,8 +651,9 @@ begin
pt.Lat := Math.EnsureRange(ALonLat.Lat, MIN_LATITUDE, MAX_LATITUDE);
pt.Lon := Math.EnsureRange(ALonLat.Lon, MIN_LONGITUDE, MAX_LONGITUDE);
px := ( TILE_SIZE / (2 * pi)) * ( IntPower (2, AWin.Zoom) ) * (pt.LonRad + pi);
py := ( TILE_SIZE / (2 * pi)) * ( IntPower (2, AWin.Zoom) ) * (pi - ln( tan(pi/4 + pt.LatRad/2) ));
factor := TILE_SIZE / TWO_PI * ZoomFactor(AWin.Zoom);
px := factor * (pt.LonRad + pi);
py := factor * (pi - ln( tan(pi/4 + pt.LatRad/2) ));
Result.x := Round(px);
Result.y := Round(py);
@@ -644,7 +672,7 @@ var
pt: TRealPoint;
cfmpx, cfmpm: Extended;
Z: Integer;
two_power_Z: Extended; // 2**Z
zoomfac: Extended; // 2**Z
begin
// https://epsg.io/3395
// https://pubs.usgs.gov/pp/1395/report.pdf, page 44
@@ -652,14 +680,14 @@ begin
pt.Lon := Math.EnsureRange(ALonLat.Lon, MIN_LONGITUDE, MAX_LONGITUDE);
Z := 23 - AWin.Zoom;
two_power_Z := IntPower(2, Z);
zoomfac := ZoomFactor(Z);
cfmpx := IntPower(2, 31);
cfmpm := cfmpx / EARTH_CIRCUMFERENCE;
px := (EARTH_CIRCUMFERENCE/2 + EARTH_EQUATORIAL_RADIUS * pt.LonRad) * cfmpm / two_power_Z;
px := (EARTH_CIRCUMFERENCE/2 + EARTH_EQUATORIAL_RADIUS * pt.LonRad) * cfmpm / zoomfac;
sny := EARTH_ECCENTRICITY * sin(pt.LatRad);
lny := tan(pi/4 + pt.LatRad/2) * power((1-sny)/(1+sny), EARTH_ECCENTRICITY/2);
py := (EARTH_CIRCUMFERENCE/2 - EARTH_EQUATORIAL_RADIUS * ln(lny)) * cfmpm / two_power_Z;
py := (EARTH_CIRCUMFERENCE/2 - EARTH_EQUATORIAL_RADIUS * ln(lny)) * cfmpm / zoomfac;
Result.x := Round(px);
Result.y := Round(py);
@@ -683,7 +711,7 @@ var
iMapWidth: Int64;
mPoint : TPoint;
begin
iMapWidth := Round(IntPower(2, AWin.Zoom)) * TILE_SIZE;
iMapWidth := round(ZoomFactor(AWin.Zoom)) * TILE_SIZE;
mPoint.X := (APoint.X - AWin.X) mod iMapWidth;
while mPoint.X < 0 do
@@ -706,7 +734,7 @@ const
MIN_LONGITUDE = -180;
MAX_LONGITUDE = 180;
var
two_power_zoom: Extended; // 2**zoom
zoomfac: Extended;
begin
// https://epsg.io/3857
// https://pubs.usgs.gov/pp/1395/report.pdf, page 41
@@ -714,9 +742,9 @@ begin
// note: coth: ** for better readability, but breaking OmniPascal in VSCode
// Result.LonRad := ( APoints.X / (( TILE_SIZE / (2*pi)) * 2**Zoom) ) - pi;
// Result.LatRad := arctan( sinh(pi - (APoints.Y/TILE_SIZE) / 2**Zoom * pi*2) );
two_power_Zoom := IntPower(2, Zoom);
Result.LonRad := ( APoint.X / (( TILE_SIZE / (2*pi)) * two_power_Zoom) ) - pi;
Result.LatRad := arctan( sinh(pi - (APoint.Y/TILE_SIZE) / two_power_Zoom * pi*2) );
zoomFac := ZoomFactor(Zoom);
Result.LonRad := ( APoint.X / (( TILE_SIZE / (2*pi)) * zoomFac) ) - pi;
Result.LatRad := arctan( sinh(pi - (APoint.Y/TILE_SIZE) / zoomFac * pi*2) );
Result.Lat := Math.EnsureRange(Result.Lat, MIN_LATITUDE, MAX_LATITUDE);
Result.Lon := Math.EnsureRange(Result.Lon, MIN_LONGITUDE, MAX_LONGITUDE);
@@ -748,19 +776,19 @@ var
Cpm: Extended;
Z: Integer;
t, phi: Extended;
two_power_Z: Extended; // 2**Z
zoomFac: Extended; // 2**Z
i: Integer;
begin
// https://epsg.io/3395
// https://pubs.usgs.gov/pp/1395/report.pdf, page 44
Z := 23 - Zoom;
two_power_Z := IntPower(2, Z);
WorldSize := Round(IntPower(2, 31));
zoomFac := ZoomFactor(Z);
WorldSize := ZoomFactor(31);
Cpm := WorldSize / EARTH_CIRCUMFERENCE;
LonRad := (APoint.x / (Cpm/two_power_Z) - EARTH_CIRCUMFERENCE/2) / EARTH_EQUATORIAL_RADIUS;
LatRad := (APoint.y / (Cpm/two_power_Z) - EARTH_CIRCUMFERENCE/2);
LonRad := (APoint.x / (Cpm/zoomFac) - EARTH_CIRCUMFERENCE/2) / EARTH_EQUATORIAL_RADIUS;
LatRad := (APoint.y / (Cpm/zoomFac) - EARTH_CIRCUMFERENCE/2);
t := pi/2 - 2*arctan(exp(-LatRad/EARTH_EQUATORIAL_RADIUS));
@@ -1330,6 +1358,8 @@ begin
ADegs := trunc(AValue) + 1;
end else
ADegs := trunc(AValue);
if AValue < 0 then
ADegs := -ADegs;
end;
procedure SplitGps(AValue: Double; out ADegs, AMins, ASecs: Double);
@@ -1347,35 +1377,52 @@ begin
ADegs := ADegs + 1;
end;
end;
if AValue < 0 then
ADegs := -ADegs;
end;
function GPSToDMS(Angle: Double): string;
begin
Result := GPSToDMS(Angle, DefaultFormatSettings);
end;
function GPSToDMS(Angle: Double; AFormatSettings: TFormatSettings): string;
var
deg, min, sec: Double;
begin
SplitGPS(Angle, deg, min, sec);
Result := Format('%.0f° %.0f'' %.1f"', [deg, min, sec]);
Result := Format('%.0f° %.0f'' %.*f"', [deg, min, DMS_Decimals, sec], AFormatSettings);
end;
function LatToStr(ALatitude: Double; DMS: Boolean): String;
begin
Result := LatToStr(ALatitude, DMS, DefaultFormatSettings);
end;
function LatToStr(ALatitude: Double; DMS: Boolean; AFormatSettings: TFormatSettings): String;
begin
if DMS then
Result := GPSToDMS(abs(ALatitude))
Result := GPSToDMS(abs(ALatitude), AFormatSettings)
else
Result := Format('%.6f°',[abs(ALatitude)]);
Result := Format('%.6f°',[abs(ALatitude)], AFormatSettings);
if ALatitude > 0 then
Result := Result + ' N'
else
if ALatitude < 0 then
Result := Result + 'E';
Result := Result + ' S';
end;
function LonToStr(ALongitude: Double; DMS: Boolean): String;
begin
Result := LonToStr(ALongitude, DMS, DefaultFormatSettings);
end;
function LonToStr(ALongitude: Double; DMS: Boolean; AFormatSettings: TFormatSettings): String;
begin
if DMS then
Result := GPSToDMS(abs(ALongitude))
Result := GPSToDMS(abs(ALongitude), AFormatSettings)
else
Result := Format('%.6f°', [abs(ALongitude)]);
Result := Format('%.6f°', [abs(ALongitude)], AFormatSettings);
if ALongitude > 0 then
Result := Result + ' E'
else if ALongitude < 0 then

1
components/lazmapviewer/source/mvtypes.pas
View File

@@ -49,6 +49,7 @@ Type
BottomRight : TRealPoint;
end;
implementation
function TRealPoint.GetLonRad: Extended;