LazMapViewer: Introduce map projections (EPSG3857, EPSG3395). Refactor calculation of projection. Issue #38279, patch by regs.

git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@7948 8e941d3f-bd1b-0410-a28a-d453659cc2b4

components/lazmapviewer/source/mvengine.pas

@@ -20,19 +20,17 @@ unit mvEngine;
 interface
 
 uses
-  Classes, SysUtils, IntfGraphics, Controls,
+  Classes, SysUtils, IntfGraphics, Controls, Math,
   mvTypes, mvJobQueue, mvMapProvider, mvDownloadEngine, mvCache, mvDragObj;
 
 const
-  EARTH_RADIUS = 6378137;
+  EARTH_EQUATORIAL_RADIUS = 6378137;
-  MIN_LATITUDE = -85.05112878;
+  EARTH_POLAR_RADIUS = 6356752.3142;
-  MAX_LATITUDE = 85.05112878;
+  EARTH_CIRCUMFERENCE = 2 * pi * EARTH_EQUATORIAL_RADIUS;
-  MIN_LONGITUDE = -180;
+  EARTH_ECCENTRICITY = sqrt(1-( exp(2*ln(EARTH_POLAR_RADIUS)) / exp(2*ln(EARTH_EQUATORIAL_RADIUS))));
-  MAX_LONGITUDE = 180;
+    // power() doesn't work in interface
-  SHIFT = 2 * pi * EARTH_RADIUS / 2.0;
 
-
+type
-Type
   TDrawTileEvent = Procedure (const TileId: TTileId; X,Y: integer;
     TileImg: TLazIntfImage) of object;
 
@@ -89,10 +87,14 @@ Type
       procedure SetUseThreads(AValue: Boolean);
       procedure SetWidth(AValue: integer);
       procedure SetZoom(AValue: integer);
-      function LonLatToMapWin(const aWin: TMapWindow; aPt: TRealPoint): TPoint;
+      function DegreesToMapPixels(const AWin: TMapWindow; ALonLat: TRealPoint): TPoint;
-      Function MapWinToLonLat(const aWin: TMapWindow; aPt : TPoint) : TRealPoint;
+      function MapPixelsToDegrees(const AWin: TMapWindow; APoint: TPoint): TRealPoint;
-      Procedure CalculateWin(var aWin: TMapWindow);
+      function PixelsToDegreesEPSG3395(APoint: TPoint; Zoom: Integer): TRealPoint;
-      Procedure Redraw(const aWin: TmapWindow);
+      function PixelsToDegreesEPSG3857(APoint: TPoint; Zoom: Integer): TRealPoint;
+      procedure CalculateWin(var AWin: TMapWindow);
+      function DegreesToPixelsEPSG3395(const AWin: TMapWindow; ALonLat: TRealPoint): TPoint;
+      function DegreesToPixelsEPSG3857(const AWin: TMapWindow; ALonLat: TRealPoint): TPoint;
+      procedure Redraw(const aWin: TMapWindow);
       function CalculateVisibleTiles(const aWin: TMapWindow) : TArea;
       function IsCurrentWin(const aWin: TMapWindow) : boolean;
     protected
@@ -106,15 +108,15 @@ Type
       constructor Create(aOwner: TComponent); override;
       destructor Destroy; override;
 
-      function AddMapProvider(OpeName: String; Url: String;
+      function AddMapProvider(OpeName: String; ProjectionType: TProjectionType; Url: String;
         MinZoom, MaxZoom, NbSvr: integer; GetSvrStr: TGetSvrStr = nil;
         GetXStr: TGetValStr = nil; GetYStr: TGetValStr = nil;
         GetZStr: TGetValStr = nil): TMapProvider;
       procedure CancelCurrentDrawing;
       procedure ClearMapProviders;
       procedure GetMapProviders(AList: TStrings);
-      function LonLatToScreen(aPt: TRealPoint): TPoint;
+      function LonLatToScreen(ALonLat: TRealPoint): TPoint;
-      function LonLatToWorldScreen(aPt: TRealPoint): TPoint;
+      function LonLatToWorldScreen(ALonLat: TRealPoint): TPoint;
       function ReadProvidersFromXML(AFileName: String; out AMsg: String): Boolean;
       procedure Redraw;
       Procedure RegisterProviders;
@@ -177,7 +179,7 @@ var
 implementation
 
 uses
-  Math, Forms, laz2_xmlread, laz2_xmlwrite, laz2_dom,
+  Forms, laz2_xmlread, laz2_xmlwrite, laz2_dom, TypInfo,
   mvJobs, mvGpsObj;
 
 type
@@ -354,9 +356,9 @@ begin
   inherited Destroy;
 end;
 
-function TMapViewerEngine.AddMapProvider(OpeName: String; Url: String;
+function TMapViewerEngine.AddMapProvider(OpeName: String; ProjectionType: TProjectionType;
-  MinZoom, MaxZoom, NbSvr: integer; GetSvrStr: TGetSvrStr; GetXStr: TGetValStr;
+  Url: String; MinZoom, MaxZoom, NbSvr: integer; GetSvrStr: TGetSvrStr;
-  GetYStr: TGetValStr; GetZStr: TGetValStr): TMapProvider;
+  GetXStr: TGetValStr; GetYStr: TGetValStr; GetZStr: TGetValStr): TMapProvider;
 var
   idx :integer;
 Begin
@@ -368,7 +370,7 @@ Begin
   end
   else
     Result := TMapProvider(lstProvider.Objects[idx]);
-  Result.AddUrl(Url, NbSvr, MinZoom, MaxZoom, GetSvrStr, GetXStr, GetYStr, GetZStr);
+  Result.AddUrl(Url, ProjectionType, NbSvr, MinZoom, MaxZoom, GetSvrStr, GetXStr, GetYStr, GetZStr);
 end;
 
 function TMapViewerEngine.CalculateVisibleTiles(const aWin: TMapWindow): TArea;
@@ -385,22 +387,18 @@ begin
   Result.Bottom := startY + MaxY;
 end;
 
-procedure TMapViewerEngine.CalculateWin(var aWin: TMapWindow);
+procedure TMapViewerEngine.CalculateWin(var AWin: TMapWindow);
 var
-  mx, my: Extended;
+  PixelLocation: TPoint; // review: coth: Should it use Int64?
-  res: Extended;
-  px, py: Int64;
 begin
-  mx := aWin.Center.Lon * SHIFT / 180.0;
+  case AWin.MapProvider.ProjectionType of
-  my := ln( tan((90 - aWin.Center.Lat) * pi / 360.0 )) / (pi / 180.0);
+    ptEPSG3857: PixelLocation := DegreesToPixelsEPSG3857(AWin, AWin.Center);
-  my := my * SHIFT / 180.0;
+    ptEPSG3395: PixelLocation := DegreesToPixelsEPSG3395(AWin, AWin.Center);
+    else PixelLocation := DegreesToPixelsEPSG3857(AWin, AWin.Center);
+  end;
 
-  res := (2 * pi * EARTH_RADIUS) / (TILE_SIZE * (1 shl aWin.Zoom));
+  AWin.X := Int64(AWin.Width div 2) - PixelLocation.x;
-  px := Round((mx + shift) / res);
+  AWin.Y := Int64(AWin.Height div 2) - PixelLocation.y;
-  py := Round((my + shift) / res);
-
-  aWin.X := aWin.Width div 2 - px;
-  aWin.Y := aWin.Height div 2 - py;
 end;
 
 procedure TMapViewerEngine.CancelCurrentDrawing;
@@ -572,88 +570,193 @@ begin
             (aTile.Y >= 0) and (aTile.Y <= tiles-1);
 end;
 
-function TMapViewerEngine.LonLatToMapWin(const aWin: TMapWindow;
+function TMapViewerEngine.DegreesToMapPixels(const AWin: TMapWindow;
-  aPt: TRealPoint): TPoint;
+  ALonLat: TRealPoint): TPoint;
 var
-  tiles: Int64;
+  pixelLocation: TPoint;
-  circumference: Int64;
-  res: Extended;
-  tmpX,tmpY : Double;
 begin
-  tiles := 1 shl aWin.Zoom;
+  case AWin.MapProvider.ProjectionType of
-  circumference := tiles * TILE_SIZE;
+    ptEPSG3395: pixelLocation := DegreesToPixelsEPSG3395(AWin, ALonLat);
-  tmpX := ((aPt.Lon+ 180.0)*circumference)/360.0;
+    ptEPSG3857: pixelLocation := DegreesToPixelsEPSG3857(AWin, ALonLat);
-
+    else pixelLocation := DegreesToPixelsEPSG3857(AWin, ALonLat);
-  res := (2 * pi * EARTH_RADIUS) / circumference;
+  end;
-
+  Result.X := pixelLocation.x + AWin.X;
-  tmpY := -aPt.Lat;
+  Result.Y := pixelLocation.y + AWin.Y;
-  tmpY := ln(tan((degToRad(tmpY) + pi / 2.0) / 2)) *180 / pi;
-  tmpY:= (((tmpY / 180.0) * SHIFT) + SHIFT) / res;
-
-  tmpX := tmpX + aWin.X;
-  tmpY := tmpY + aWin.Y;
-  Result.X := trunc(tmpX);
-  Result.Y := trunc(tmpY);
 end;
 
-function TMapViewerEngine.LonLatToScreen(aPt: TRealPoint): TPoint;
+// review: coth: Should it use Int64?
+function TMapViewerEngine.DegreesToPixelsEPSG3857(const AWin: TMapWindow;
+  ALonLat: TRealPoint): TPoint;
+const
+  MIN_LATITUDE = -85.05112878;
+  MAX_LATITUDE = 85.05112878;
+  MIN_LONGITUDE = -180;
+  MAX_LONGITUDE = 180;
+var
+  px, py: Extended;
+  pt: TRealPoint;
+begin
+  // https://epsg.io/3857
+  // https://pubs.usgs.gov/pp/1395/report.pdf, page 41
+  // https://en.wikipedia.org/wiki/Web_Mercator_projection
+  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) ));
+
+  Result.x := Round(px);
+  Result.y := Round(py);
+end;
+
+// review: coth: Should it use Int64?
+function TMapViewerEngine.DegreesToPixelsEPSG3395(const AWin: TMapWindow;
+  ALonLat: TRealPoint): TPoint;
+const
+  MIN_LATITUDE = -80;
+  MAX_LATITUDE = 84;
+  MIN_LONGITUDE = -180;
+  MAX_LONGITUDE = 180;
+var
+  px, py, lny, sny: Extended;
+  pt: TRealPoint;
+  cfmpx, cfmpm: Extended;
+  Z: Integer;
+  two_power_Z: Extended;  // 2**Z
+begin
+  // https://epsg.io/3395
+  // https://pubs.usgs.gov/pp/1395/report.pdf, page 44
+  pt.Lat := Math.EnsureRange(ALonLat.Lat, MIN_LATITUDE, MAX_LATITUDE);
+  pt.Lon := Math.EnsureRange(ALonLat.Lon, MIN_LONGITUDE, MAX_LONGITUDE);
+
+  Z := 23 - AWin.Zoom;
+  two_power_Z := IntPower(2, Z);
+  cfmpx := IntPower(2, 31);
+  cfmpm := cfmpx / EARTH_CIRCUMFERENCE;
+  px := (EARTH_CIRCUMFERENCE/2 + EARTH_EQUATORIAL_RADIUS * pt.LonRad) * cfmpm / two_power_Z;
+
+  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;
+
+  Result.x := Round(px);
+  Result.y := Round(py);
+end;
+
+
+function TMapViewerEngine.LonLatToScreen(ALonLat: TRealPoint): TPoint;
 Begin
-  Result := LonLatToMapWin(MapWin, aPt);
+  Result := DegreesToMapPixels(MapWin, ALonLat);
 end;
 
-function TMapViewerEngine.LonLatToWorldScreen(aPt: TRealPoint): TPoint;
+function TMapViewerEngine.LonLatToWorldScreen(ALonLat: TRealPoint): TPoint;
 begin
-  Result := LonLatToScreen(aPt);
+  Result := LonLatToScreen(ALonLat);
   Result.X := Result.X + MapWin.X;
   Result.Y := Result.Y + MapWin.Y;
 end;
 
-function TMapViewerEngine.MapWinToLonLat(const aWin: TMapWindow;
+function TMapViewerEngine.MapPixelsToDegrees(const AWin: TMapWindow;
-  aPt: TPoint): TRealPoint;
+  APoint: TPoint): TRealPoint;
 var
-  tiles: Int64;
+  iMapWidth: Int64;
-  circumference: Int64;
-  lat: Extended;
-  res: Extended;
   mPoint : TPoint;
 begin
-  tiles := 1 shl aWin.Zoom;
+  // review: coth: respective projection check. move to subfunctions? figure out what did i mean here...
-  circumference := tiles * TILE_SIZE;
+  iMapWidth := Round(IntPower(2, AWin.Zoom)) * TILE_SIZE;
 
-  mPoint.X := aPt.X - aWin.X;
+  mPoint.X := EnsureRange(APoint.X - AWin.X, 0, iMapWidth);
-  mPoint.Y := aPt.Y - aWin.Y;
+  mPoint.Y := EnsureRange(APoint.Y - AWin.Y, 0, iMapWidth);
 
-  if mPoint.X < 0 then
+  case aWin.MapProvider.ProjectionType of
-    mPoint.X := 0
+    ptEPSG3857: Result := PixelsToDegreesEPSG3857(mPoint, AWin.Zoom);
-  else
+    ptEPSG3395: Result := PixelsToDegreesEPSG3395(mPoint, AWin.Zoom);
-  if mPoint.X > circumference then
+    else Result := PixelsToDegreesEPSG3857(mPoint, AWin.Zoom);
-    mPoint.X := circumference;
+  end;
+end;
 
-  if mPoint.Y < 0 then
+function TMapViewerEngine.PixelsToDegreesEPSG3857(APoint: TPoint; Zoom: Integer): TRealPoint;
-    mPoint.Y := 0
+const
-  else
+  MIN_LATITUDE = -85.05112878;
-  if mPoint.Y > circumference then
+  MAX_LATITUDE = 85.05112878;
-    mPoint.Y := circumference;
+  MIN_LONGITUDE = -180;
+  MAX_LONGITUDE = 180;
+var
+  two_power_zoom: Extended;   // 2**zoom
+begin
+  // https://epsg.io/3857
+  // https://pubs.usgs.gov/pp/1395/report.pdf, page 41
 
-  Result.Lon := ((mPoint.X * 360.0) / circumference) - 180.0;
+  // 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) );
 
-  res := (2 * pi * EARTH_RADIUS) / circumference;
+  Result.Lat := Math.EnsureRange(Result.Lat, MIN_LATITUDE, MAX_LATITUDE);
-  lat := ((mPoint.Y * res - SHIFT) / SHIFT) * 180.0;
+  Result.Lon := Math.EnsureRange(Result.Lon, MIN_LONGITUDE, MAX_LONGITUDE);
+end;
 
-  lat := radtodeg (2 * arctan( exp( lat * pi / 180.0)) - pi / 2.0);
+Function TMapViewerEngine.PixelsToDegreesEPSG3395(APoint: TPoint; Zoom: Integer): TRealPoint;
-  Result.Lat := -lat;
 
-  if Result.Lat > MAX_LATITUDE then
+  function PhiIteration(y, phi: Extended): Extended;
-    Result.Lat := MAX_LATITUDE
+  var
-  else
+    t: Extended;
-  if Result.Lat < MIN_LATITUDE then
+    sin_phi: Extended;
-    Result.Lat := MIN_LATITUDE;
+    arg: Extended;
+  begin
+    t := exp(y/EARTH_EQUATORIAL_RADIUS);
+    sin_phi := sin(phi);
+    arg := (1 - EARTH_ECCENTRICITY * sin_phi) / (1 + EARTH_ECCENTRICITY * sin_phi);
+    Result := pi/2 - 2*arctan( t * Math.power(arg, EARTH_ECCENTRICITY/2) );
+  end;
 
-  if Result.Lon > MAX_LONGITUDE then
+const
-    Result.Lon := MAX_LONGITUDE
+  MIN_LATITUDE = -80;
-  else
+  MAX_LATITUDE = 84;
-  if Result.Lon < MIN_LONGITUDE then
+  MIN_LONGITUDE = -180;
-    Result.Lon := MIN_LONGITUDE;
+  MAX_LONGITUDE = 180;
+  EPS = 1e-8;
+var
+  LonRad, LatRad: Extended;
+  WorldSize: Int64;
+  Cpm: Extended;
+  Z: Integer;
+  t, phi:  Extended;
+  two_power_Z: 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));
+  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);
+
+  t := pi/2 - 2*arctan(exp(-LatRad/EARTH_EQUATORIAL_RADIUS));
+
+  i := 0;
+  repeat
+    phi := t;
+    t := PhiIteration(LatRad, phi);
+    inc(i);
+    if i>10 then
+      Break;
+      //raise Exception.Create('Phi iteration takes too long.');
+  until (abs(phi - t) < EPS);
+
+  LatRad := t;
+
+  Result.LonRad := LonRad;
+  Result.LatRad := LatRad;
+
+  Result.Lat := Math.EnsureRange(Result.Lat, MIN_LATITUDE, MAX_LATITUDE);
+  Result.Lon := Math.EnsureRange(Result.Lon, MIN_LONGITUDE, MAX_LONGITUDE);
 end;
 
 procedure TMapViewerEngine.MouseDown(Sender: TObject; Button: TMouseButton;
@@ -705,7 +808,7 @@ Begin
   old := TMemObj(Sender.LnkObj);
   aPt.X := old.FWin.Width DIV 2-Sender.OfsX;
   aPt.Y := old.FWin.Height DIV 2-Sender.OfsY;
-  nCenter := MapWinToLonLat(old.FWin,aPt);
+  nCenter := MapPixelsToDegrees(old.FWin,aPt);
   SetCenter(nCenter);
 end;
 
@@ -756,6 +859,7 @@ var
   doc: TXMLDocument = nil;
   node, layerNode: TDOMNode;
   providerName: String;
+  projectionType: TProjectionType;
   url: String;
   minZoom: Integer;
   maxZoom: Integer;
@@ -796,6 +900,8 @@ begin
         s := GetAttrValue(layerNode, 'serverCount');
         if s = '' then svrCount := 1
           else svrCount := StrToInt(s);
+        s := Concat('pt', GetAttrValue(layerNode, 'projection'));
+        projectionType := TProjectionType(GetEnumValue(TypeInfo(TProjectionType), s)); //-1 will default to ptEPSG3857
         svrProc := GetAttrValue(layerNode, 'serverProc');
         xProc := GetAttrValue(layerNode, 'xProc');
         yProc := GetAttrValue(layerNode, 'yProc');
@@ -806,7 +912,7 @@ begin
         ClearMapProviders;
         first := false;
       end;
-      AddMapProvider(providerName,
+      AddMapProvider(providerName, projectionType,
         url, minZoom, maxZoom, svrCount,
         GetSvrStr(svrProc), GetValStr(xProc), GetValStr(yProc), GetValStr(zProc)
       );
@@ -855,17 +961,30 @@ procedure TMapViewerEngine.RegisterProviders;
 var
   HERE1, HERE2: String;
 begin
-//  AddMapProvider('Aucun','',0,30, 0);  ???
 
-  AddMapProvider('Google Normal',
+  // dev links
-    'http://mt%serv%.google.com/vt/lyrs=m@145&v=w2.104&x=%x%&y=%y%&z=%z%',
+  //https://gis-lab.info/forum/viewtopic.php?f=19&t=19763
-    0, 19, 4, nil);
+  //https://a.tile.openstreetmap.org/16/51693/32520.png
-  AddMapProvider('Google Hybrid',
+  //https://vec01.maps.yandex.net/tiles?l=map&x=51693+570&y=32520&z=16&scale=1&lang=ru_RU
-    'http://mt%serv%.google.com/vt/lyrs=h@145&v=w2.104&x=%x%&y=%y%&z=%z%',
+  //https://www.linux.org.ru/forum/development/9038716
-    0, 19, 4, nil);
+  //https://wiki.openstreetmap.org/wiki/Tiles
-  AddMapProvider('Google Physical',
+  //https://pubs.usgs.gov/pp/1395/report.pdf
-    'http://mt%serv%.google.com/vt/lyrs=t@145&v=w2.104&x=%x%&y=%y%&z=%z%',
+  //https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Tile_numbers_to_lon..2Flat.
-    0, 19, 4, nil);
+  //https://mc.bbbike.org/mc/?num=2
+  //https://mc.bbbike.org/mc/?lon=37.62178&lat=55.740937&zoom=14&num=1&mt0=opentopomap&mt1=mapnik-german
+  //https://t.ssl.ak.dynamic.tiles.virtualearth.net/comp/ch/12031010103311?mkt=ru-RU&it=G,BX,RL&shading=hill&n=z&og=677&c4w=1&cstl=vb&src=h
+
+  {
+  // needs hybrid overlays
+  AddMapProvider('Google Hybrid', ptEPSG3857, 'http://mt%serv%.google.com/vt/lyrs=h@145&v=w2.104&x=%x%&y=%y%&z=%z%', 0, 19, 4, nil);
+  AddMapProvider('Yandex.Maps Hybrid', ptEPSG3395, 'https://vec0%serv%.maps.yandex.net/tiles?l=skl&x=%x%&y=%y%&z=%z%', 0, 19, 4, @GetSvrBase1, nil, nil, nil);
+  }
+
+  // Google
+  AddMapProvider('Google Maps', ptEPSG3857, 'http://mt%serv%.google.com/vt/lyrs=m@145&v=w2.104&x=%x%&y=%y%&z=%z%', 0, 19, 4, nil);
+  AddMapProvider('Google Sattelite', ptEPSG3857, 'http://khm%serv%.google.com/kh/v=863?x=%x%&y=%y%&z=%z%', 0, 19, 4, nil);
+  //AddMapProvider('Google Physical', ptEPSG3857, 'http://mt%serv%.google.com/vt/lyrs=t@145&v=w2.104&x=%x%&y=%y%&z=%z%', 0, 19, 4, nil); // review: doesn't work?
+
 
   {
   AddMapProvider('Google Hybrid','http://khm%d.google.com/kh/v=82&x=%x%&y=%y%&z=%z%&s=Ga',4);
@@ -879,7 +998,7 @@ begin
   //AddMapProvider('Yahoo Satellite','http://maps%serv%.yimg.com/ae/ximg?v=1.9&t=a&s=256&.intl=en&x=%d&y=%d&z=%d&r=1', 0,20,3,@GetYahooSvr, nil, @getYahooY, @GetYahooZ); //[Random(3)+1, X, YahooY(Y), Z+1]));
   //AddMapProvider('Yahoo Hybrid','http://maps%serv%.yimg.com/ae/ximg?v=1.9&t=a&s=256&.intl=en&x=%x%&y=%y%&z=%z%&r=1', 0,20,3,@GetYahooSvr, nil, @getYahooY, @GetYahooZ); //[Random(3)+1, X, YahooY(Y), Z+1]));
   //AddMapProvider('Yahoo Hybrid','http://maps%serv%.yimg.com/hx/tl?b=1&v=4.3&t=h&.intl=en&x=%x%&y=%y%&z=%z%&r=1'    , 0,20,3,@GetYahooSvr, nil, @getYahooY, @GetYahooZ); //[Random(3)+1, X, YahooY(Y), Z+1]));
-
+                               (*
   // opeName, Url, MinZoom, MaxZoom, NbSvr, GetSvrStr, GetXStr, GetYStr, GetZStr
   MapWin.MapProvider := AddMapProvider('OpenStreetMap Mapnik',
     'http://%serv%.tile.openstreetmap.org/%z%/%x%/%y%.png',
@@ -890,6 +1009,7 @@ begin
   AddMapProvider('Open Topo Map',
     'http://%serv%.tile.opentopomap.org/%z%/%x%/%y%.png',
     0, 19, 3, @GetSvrLetter);
+
   AddMapProvider('Virtual Earth Bing',
     'http://ecn.t%serv%.tiles.virtualearth.net/tiles/r%x%?g=671&mkt=en-us&lbl=l1&stl=h&shading=hill',
     1, 19, 8, nil, @GetStrQuadKey);
@@ -902,6 +1022,24 @@ begin
   AddMapProvider('Virtual Earth Hybrid',
     'http://h%serv%.ortho.tiles.virtualearth.net/tiles/h%x%.jpg?g=72&shading=hill',
     1, 19, 4, nil, @GetStrQuadKey);
+                            *)
+  // OpenStreetMap section
+  MapWin.MapProvider := AddMapProvider('OpenStreetMap Mapnik', ptEPSG3857, 'http://%serv%.tile.openstreetmap.org/%z%/%x%/%y%.png', 0, 19, 3, @GetSvrLetter);
+  AddMapProvider('OpenStreetMap Wikipedia', ptEPSG3857, 'https://maps.wikimedia.org/osm-intl/%z%/%x%/%y%.png', 0, 19, 3, @GetSvrLetter);
+  AddMapProvider('OpenStreetMap Sputnik', ptEPSG3857, 'https://%serv%.tilessputnik.ru/tiles/kmt2/%z%/%x%/%y%.png', 0, 19, 3, @GetSvrLetter);
+  AddMapProvider('OpenStreetMap.fr Hot', ptEPSG3857, 'https://%serv%.tile.openstreetmap.fr/hot/%z%/%x%/%y%.png', 0, 18, 3, @GetSvrLetter);
+  AddMapProvider('Open Topo Map', ptEPSG3857, 'http://%serv%.tile.opentopomap.org/%z%/%x%/%y%.png', 0, 19, 3, @GetSvrLetter);
+  AddMapProvider('OpenStreetMap.fr Cycle Map', ptEPSG3857, 'https://dev.%serv%.tile.openstreetmap.fr/cyclosm/%z%/%x%/%y%.png', 0, 18, 3, @GetSvrLetter);
+  // todo: requires an optional key
+  AddMapProvider('Open Cycle Map', ptEPSG3857, 'http://%serv%.tile.opencyclemap.org/cycle/%z%/%x%/%y%.png', 0, 18, 3, @GetSvrLetter);
+  AddMapProvider('OpenStreetMap Transport', ptEPSG3857, 'https://%serv%.tile.thunderforest.com/transport/%z%/%x%/%y%.png', 0, 18, 3, @GetSvrLetter);
+
+  // Bing
+  AddMapProvider('Virtual Earth Bing', ptEPSG3857, 'http://ecn.t%serv%.tiles.virtualearth.net/tiles/r%x%?g=671&mkt=en-us&lbl=l1&stl=h&shading=hill', 1, 19, 8, nil, @GetStrQuadKey);
+  // review: remove? fully identical to Bing Maps?
+  //AddMapProvider('Virtual Earth Road', ptEPSG3857, 'http://r%serv%.ortho.tiles.virtualearth.net/tiles/r%x%.png?g=72&shading=hill', 1, 19, 4, nil, @GetStrQuadKey);
+  AddMapProvider('Virtual Earth Aerial', ptEPSG3857, 'http://a%serv%.ortho.tiles.virtualearth.net/tiles/a%x%.jpg?g=72&shading=hill', 1, 19, 4, nil, @GetStrQuadKey);
+  AddMapProvider('Virtual Earth Hybrid', ptEPSG3857, 'http://h%serv%.ortho.tiles.virtualearth.net/tiles/h%x%.jpg?g=72&shading=hill', 1, 19, 4, nil, @GetStrQuadKey);
 
   if (HERE_AppID <> '') and (HERE_AppCode <> '') then begin
     // Registration required to access HERE maps:
@@ -911,20 +1049,13 @@ begin
     // restart the demo.
     HERE1 := 'http://%serv%.base.maps.api.here.com/maptile/2.1/maptile/newest/';
     HERE2 := '/%z%/%x%/%y%/256/png8?app_id=' + HERE_AppID + '&app_code=' + HERE_AppCode;
-    AddMapProvider('Here Maps', HERE1 + 'normal.day' + HERE2,
+    AddMapProvider('Here WeGo Map', ptEPSG3857, HERE1 + 'normal.day' + HERE2, 1, 19, 4, @GetSvrBase1);
-      1, 19, 4, @GetSvrBase1);
+    AddMapProvider('Here WeGo Grey Map', ptEPSG3857, HERE1 + 'normal.day.grey' + HERE2, 1, 19, 4, @GetSvrBase1);
-    AddMapProvider('Here Maps Grey', HERE1 + 'normal.day.grey' + HERE2,
+    AddMapProvider('Here WeGo Reduced Map', ptEPSG3857, HERE1 + 'reduced.day' + HERE2, 1, 19, 4, @GetSvrBase1);
-      1, 19, 4, @GetSvrBase1);
+    AddMapProvider('Here WeGo Transit Map', ptEPSG3857, HERE1 + 'normal.day.transit' + HERE2, 1, 19, 4, @GetSvrBase1);
-    AddMapProvider('Here Maps Reduced', HERE1 + 'reduced.day' + HERE2,
+    AddMapProvider('Here WeGo POI Map', ptEPSG3857, HERE1 + 'normal.day' + HERE2 + '&pois', 1, 19, 4, @GetSvrBase1);
-      1, 19, 4, @GetSvrBase1);
+    AddMapProvider('Here WeGo Pedestrian Map', ptEPSG3857, HERE1 + 'pedestrian.day' + HERE2, 1, 19, 4, @GetSvrBase1);
-    AddMapProvider('Here Maps Transit', HERE1 + 'normal.day.transit' + HERE2,
+    AddMapProvider('Here WeGo DreamWorks Map', ptEPSG3857, HERE1 + 'normal.day' + HERE2 + '&style=dreamworks', 1, 19, 4, @GetSvrBase1);
-      1, 19, 4, @GetSvrBase1);
-    AddMapProvider('Here POI Maps', HERE1 + 'normal.day' + HERE2 + '&pois',
-      1, 19, 4, @GetSvrBase1);
-    AddMapProvider('Here Pedestrian Maps', HERE1 + 'pedestrian.day' + HERE2,
-      1, 19, 4, @GetSvrBase1);
-    AddMapProvider('Here DreamWorks Maps', HERE1 + 'normal.day' + HERE2 + '&style=dreamworks',
-      1, 19, 4, @GetSvrBase1);
   end;
 
   if (OpenWeatherMap_ApiKey <> '') then begin
@@ -932,50 +1063,18 @@ begin
     //   https://home.openweathermap.org/users/sign_up
     // Store the API key found on the website in the ini file of the demo under
     // key [OpenWeatherMap] and API_Key and restart the demo
-    AddMapProvider('OpenWeatherMap Clouds',
+    AddMapProvider('OpenWeatherMap Clouds', ptEPSG3857, 'https://tile.openweathermap.org/map/clouds_new/%z%/%x%/%y%.png?appid=' + OpenWeatherMap_ApiKey, 1, 19, 1, nil);
-      'https://tile.openweathermap.org/map/clouds_new/%z%/%x%/%y%.png?appid=' + OpenWeatherMap_ApiKey,
+    AddMapProvider('OpenWeatherMap Precipitation', ptEPSG3857, 'https://tile.openweathermap.org/map/precipitation_new/%z%/%x%/%y%.png?appid=' + OpenWeatherMap_ApiKey, 1, 19, 1, nil);
-      1, 19, 1, nil);
+    AddMapProvider('OpenWeatherMap Pressure', ptEPSG3857, 'https://tile.openweathermap.org/map/pressure_new/%z%/%x%/%y%.png?appid=' + OpenWeatherMap_ApiKey, 1, 19, 1, nil);
-    AddMapProvider('OpenWeatherMap Precipitation',
+    AddMapProvider('OpenWeatherMap Temperature', ptEPSG3857, 'https://tile.openweathermap.org/map/temp_new/%z%/%x%/%y%.png?appid=' + OpenWeatherMap_ApiKey, 1, 19, 1, nil);
-      'https://tile.openweathermap.org/map/precipitation_new/%z%/%x%/%y%.png?appid=' + OpenWeatherMap_ApiKey,
+    AddMapProvider('OpenWeatherMap Wind', ptEPSG3857, 'https://tile.openweathermap.org/map/wind_new/%z%/%x%/%y%.png?appid=' + OpenWeatherMap_ApiKey, 1, 19, 1, nil);
-      1, 19, 1, nil);
-    AddMapProvider('OpenWeatherMap Pressure',
-      'https://tile.openweathermap.org/map/pressure_new/%z%/%x%/%y%.png?appid=' + OpenWeatherMap_ApiKey,
-      1, 19, 1, nil);
-    AddMapProvider('OpenWeatherMap Temperature',
-      'https://tile.openweathermap.org/map/temp_new/%z%/%x%/%y%.png?appid=' + OpenWeatherMap_ApiKey,
-      1, 19, 1, nil);
-    AddMapProvider('OpenWeatherMap Wind',
-      'https://tile.openweathermap.org/map/wind_new/%z%/%x%/%y%.png?appid=' + OpenWeatherMap_ApiKey,
-      1, 19, 1, nil);
   end;
 
-  { The Ovi Maps (former Nokia maps) are no longer available.
-
-  AddMapProvider('Ovi Normal',
-    'http://%serv%.maptile.maps.svc.ovi.com/maptiler/v2/maptile/newest/normal.day/%z%/%x%/%y%/256/png8',
-    0, 20, 5, @GetLetterSvr);
-  AddMapProvider('Ovi Satellite',
-    'http://%serv%.maptile.maps.svc.ovi.com/maptiler/v2/maptile/newest/satellite.day/%z%/%x%/%y%/256/png8',
-    0, 20, 5, @GetLetterSvr);
-  AddMapProvider('Ovi Hybrid',
-    'http://%serv%.maptile.maps.svc.ovi.com/maptiler/v2/maptile/newest/hybrid.day/%z%/%x%/%y%/256/png8',
-    0, 20, 5, @GetLetterSvr);
-  AddMapProvider('Ovi Physical',
-    'http://%serv%.maptile.maps.svc.ovi.com/maptiler/v2/maptile/newest/terrain.day/%z%/%x%/%y%/256/png8',
-    0, 20, 5, @GetLetterSvr);
-  }
-
-  {
-  AddMapProvider('Yahoo Normal','http://maps%serv%.yimg.com/hx/tl?b=1&v=4.3&.intl=en&x=%x%&y=%y%d&z=%d&r=1'        , 0,20,3,@GetYahooSvr, nil, @getYahooY, @GetYahooZ); //(Z+1]));
-  AddMapProvider('Yahoo Satellite','http://maps%serv%.yimg.com/ae/ximg?v=1.9&t=a&s=256&.intl=en&x=%d&y=%d&z=%d&r=1', 0,20,3,@GetYahooSvr, nil, @getYahooY, @GetYahooZ); //[Random(3)+1, X, YahooY(Y), Z+1]));
-  AddMapProvider('Yahoo Hybrid','http://maps%serv%.yimg.com/ae/ximg?v=1.9&t=a&s=256&.intl=en&x=%x%&y=%y%&z=%z%&r=1', 0,20,3,@GetYahooSvr, nil, @getYahooY, @GetYahooZ); //[Random(3)+1, X, YahooY(Y), Z+1]));
-  AddMapProvider('Yahoo Hybrid','http://maps%serv%.yimg.com/hx/tl?b=1&v=4.3&t=h&.intl=en&x=%x%&y=%y%&z=%z%&r=1'    , 0,20,3,@GetYahooSvr, nil, @getYahooY, @GetYahooZ); //[Random(3)+1, X, YahooY(Y), Z+1]));
-  }
 end;
 
 function TMapViewerEngine.ScreenToLonLat(aPt: TPoint): TRealPoint;
 begin
-  Result := MapWinToLonLat(MapWin, aPt);
+  Result := MapPixelsToDegrees(MapWin, aPt);
 end;
 
 procedure TMapViewerEngine.SetActive(AValue: boolean);
@@ -1043,6 +1142,7 @@ begin
   begin
     MapWin.MapProvider := TMapProvider(lstProvider.Objects[idx]);
     ConstraintZoom(MapWin);
+    CalculateWin(MapWin);
   end
   else
     MapWin.MapProvider := nil;
@@ -1102,7 +1202,7 @@ begin
     begin
        Cache.GetFromCache(EnvTile.Win.MapProvider, EnvTile.Tile, img);
        X := EnvTile.Win.X + EnvTile.Tile.X * TILE_SIZE; // begin of X
-       Y := EnvTile.Win.Y + EnvTile.Tile.Y * TILE_SIZE; // begin of X
+       Y := EnvTile.Win.Y + EnvTile.Tile.Y * TILE_SIZE; // begin of Y
        DrawTile(EnvTile.Tile, X, Y, img);
     end;
   finally
@@ -1154,8 +1254,8 @@ begin
   BottomRight.Y := tmpWin.Height;
   Repeat
     CalculateWin(tmpWin);
-    visArea.TopLeft := MapWinToLonLat(tmpWin, TopLeft);
+    visArea.TopLeft := MapPixelsToDegrees(tmpWin, TopLeft);
-    visArea.BottomRight := MapWinToLonLat(tmpWin, BottomRight);
+    visArea.BottomRight := MapPixelsToDegrees(tmpWin, BottomRight);
     if AreaInsideArea(aArea, visArea) then
       break;
     dec(tmpWin.Zoom);
@@ -1224,8 +1324,8 @@ end;
 
   Each part can exhibit a unit identifier, such as °, ', or ". BUT: they are
   ignored. This means that an input string 50°30" results in the output value 50.5
-  although the second part is marked as seconds, not minutes! 
+  although the second part is marked as seconds, not minutes!
-  
+
   Hemisphere suffixes ('N', 'S', 'E', 'W') are supported at the end of the input string.
 }
 function TryStrToGps(const AValue: String; out ADeg: Double): Boolean;
@@ -1368,7 +1468,7 @@ begin
     d_radians := PI * 2.0
   else
     d_radians := arccos(arg);
-  Result := EARTH_RADIUS * d_radians;
+  Result := EARTH_EQUATORIAL_RADIUS * d_radians;
 
   case AUnits of
     duMeters: ;

components/lazmapviewer/source/mvmapprovider.pas

@@ -28,6 +28,7 @@ type
 
   TGetSvrStr = function (id: integer): string;
   TGetValStr = function (const Tile: TTileId): String;
+  TProjectionType = (ptEPSG3857, ptEPSG3395);
 
   TMapProvider = class;
 
@@ -49,6 +50,7 @@ type
       idServer: Array of Integer;
       FName: String;
       FUrl: Array of string;
+      FProjectionType: Array of TProjectionType;
       FNbSvr: Array of integer;
       FGetSvrStr: Array of TGetSvrStr;
       FGetXStr: Array of TGetValStr;
@@ -59,13 +61,14 @@ type
       FTiles:array of TBaseTile;
       FTileHandling: TRTLCriticalSection;
       function GetLayerCount: integer;
+      function GetProjectionType: TProjectionType;
       procedure SetLayer(AValue: integer);
     public
       constructor Create(AName: String);
       destructor Destroy; override;
       function AppendTile(aTile: TBaseTile): integer;
       procedure RemoveTile(aTile: TBaseTile);
-      procedure AddURL(Url: String; NbSvr, aMinZoom, aMaxZoom: integer;
+      procedure AddURL(Url: String; ProjectionType: TProjectionType; NbSvr, aMinZoom, aMaxZoom: integer;
         GetSvrStr: TGetSvrStr; GetXStr: TGetValStr; GetYStr: TGetValStr;
         GetZStr: TGetValStr);
       procedure GetZoomInfos(out AZoomMin, AZoomMax: integer);
@@ -74,6 +77,7 @@ type
       property Name: String read FName;
       property LayerCount: integer read GetLayerCount;
       property Layer: integer read FLayer write SetLayer;
+      property ProjectionType: TProjectionType read GetProjectionType;
   end;
 
 
@@ -93,6 +97,9 @@ const
 
 implementation
 
+uses
+  TypInfo;
+
 function GetSvrLetter(id: integer): String;
 begin
   Result := Char(Ord('a') + id);
@@ -154,7 +161,12 @@ end;
 
 function TMapProvider.GetLayerCount: integer;
 begin
-  Result:=length(FUrl);
+  Result := Length(FUrl);
+end;
+
+function TMapProvider.GetProjectionType: TProjectionType;
+begin
+  Result := FProjectionType[layer];
 end;
 
 procedure TMapProvider.SetLayer(AValue: integer);
@@ -179,6 +191,7 @@ var
 begin
   Finalize(idServer);
   Finalize(FName);
+  Finalize(FProjectionType);
   Finalize(FUrl);
   Finalize(FNbSvr);
   Finalize(FGetSvrStr);
@@ -228,15 +241,16 @@ begin
     end;
 end;
 
-procedure TMapProvider.AddURL(Url: String; NbSvr, aMinZoom, aMaxZoom: integer;
+procedure TMapProvider.AddURL(Url: String; ProjectionType: TProjectionType;
-  GetSvrStr: TGetSvrStr; GetXStr: TGetValStr; GetYStr: TGetValStr;
+  NbSvr, aMinZoom, aMaxZoom: integer; GetSvrStr: TGetSvrStr;
-  GetZStr: TGetValStr);
+  GetXStr: TGetValStr; GetYStr: TGetValStr; GetZStr: TGetValStr);
 var
   nb: integer;
 begin
   nb := Length(FUrl)+1;
   SetLength(IdServer, nb);
   SetLength(FUrl, nb);
+  SetLength(FProjectionType, nb);
   SetLength(FNbSvr, nb);
   SetLength(FGetSvrStr, nb);
   SetLength(FGetXStr, nb);
@@ -246,6 +260,7 @@ begin
   SetLength(FMaxZoom, nb);
   nb := High(FUrl);
   FUrl[nb] := Url;
+  FProjectionType[nb] := ProjectionType;
   FNbSvr[nb] := NbSvr;
   FMinZoom[nb] := aMinZoom;
   FMaxZoom[nb] := aMaxZoom;
@@ -304,6 +319,7 @@ begin
   node := ADoc.CreateElement('map_provider');
   node.SetAttribute('name', FName);
   AParentNode.AppendChild(node);
+
   for i:=0 to LayerCount-1 do begin
     layerNode := ADoc.CreateElement('layer');
     node.AppendChild(layernode);
@@ -312,6 +328,11 @@ begin
     layerNode.SetAttribute('maxZoom', IntToStr(FMaxZoom[i]));
     layerNode.SetAttribute('serverCount', IntToStr(FNbSvr[i]));
 
+    s := GetEnumName(TypeInfo(TProjectionType), Ord(FProjectionType[i]));
+    if s.StartsWith('pt') then
+      s := s.Substring(2);
+    layerNode.SetAttribute('projection', s);
+
     if FGetSvrStr[i] = @GetSvrLetter then
       s := SVR_LETTER
     else if FGetSvrStr[i] = @GetSvrBase1 then