{
  (c) 2014 ti_dic@hotmail.com

  Parts of this component are based on :
    Map Viewer Copyright (C) 2011 Maciej Kaczkowski / keit.co

  License: modified LGPL with linking exception (like RTL, FCL and LCL)

  See the file COPYING.modifiedLGPL.txt, included in the Lazarus distribution,
  for details about the license.

  See also: https://wiki.lazarus.freepascal.org/FPC_modified_LGPL
}


unit mvEngine;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, IntfGraphics, Controls, Math, GraphType, FPImage,
  mvTypes, mvJobQueue, mvMapProvider, mvDownloadEngine, mvCache, mvDragObj;

const
  EARTH_EQUATORIAL_RADIUS = 6378137;
  EARTH_POLAR_RADIUS = 6356752.3142;
  EARTH_CIRCUMFERENCE = 2 * pi * EARTH_EQUATORIAL_RADIUS;
  EARTH_ECCENTRICITY = sqrt(1 - sqr(EARTH_POLAR_RADIUS / EARTH_EQUATORIAL_RADIUS));

type
  TDrawTileEvent = procedure (const TileId: TTileId; X,Y: integer;
    TileImg: TPictureCacheItem) of object;

  TDrawStretchedTileEvent = procedure (const TileId: TTileId; X,Y: Integer;
    TileImg: TPictureCacheItem; const R: TRect) of object;

  TEraseBackgroundEvent = procedure (const R: TRect) of Object;

  TTileDownloadedEvent = procedure (const TileId: TTileId) of object;

  TTileIdArray = Array of TTileId;

  TDistanceUnits = (duMeters, duKilometers, duMiles);

  { TMapWindow }

  TMapWindow = Record
    MapProvider: TMapProvider;
    X: Int64;
    Y: Int64;
    Center: TRealPoint;
    Zoom: integer;
    ZoomCenter: TRealPoint;
    ZoomOffset: TPoint;
    Height: integer;
    Width: integer;
  end;


  { TMapViewerEngine }

  TMapViewerEngine = Class(TComponent)
    private
      FDragObj : TDragObj;
      Cache : TPictureCache;
      FActive: boolean;
      FBkColor: TFPColor;
      FCyclic: Boolean;
      FDownloadEngine: TMvCustomDownloadEngine;
      FDrawPreviewTiles: Boolean;
      FDrawTitleInGuiThread: boolean;
      FInDrag: Boolean;
      FOnCenterMove: TNotifyEvent;
      FOnChange: TNotifyEvent;
      FOnDrawTile: TDrawTileEvent;
      FOnDrawStretchedTile: TDrawStretchedTileEvent;
      FOnEraseBackground: TEraseBackgroundEvent;
      FOnTileDownloaded: TTileDownloadedEvent;
      FOnZoomChange: TNotifyEvent;
      FZoomMax: Integer;
      FZoomMin: Integer;
      lstProvider : TStringList;
      Queue : TJobQueue;
      MapWin : TMapWindow;
      FZoomToCursor: Boolean;
      function GetCacheItemClass: TPictureCacheItemClass;
      function GetCacheMaxAge: Integer;
      function GetCacheOnDisk: Boolean;
      function GetCachePath: String;
      function GetCenter: TRealPoint;
      function GetHeight: integer;
      class function GetProjectionType(const aWin: TMapWindow): TProjectionType;
      function GetMapProjectionType: TProjectionType;
      function GetMapProvider: String;
      function GetUseThreads: Boolean;
      function GetWidth: integer;
      function GetZoom: integer;
      function IsValidTile(const aWin: TMapWindow; const aTile: TTIleId): boolean;
      procedure MoveMapCenter(Sender: TDragObj);
      procedure SetActive(AValue: boolean);
      procedure SetBkColor(AValue: TFPColor);
      procedure SetCacheItemClass(AValue: TPictureCacheItemClass);
      procedure SetCacheMaxAge(AValue: Integer);
      procedure SetCacheOnDisk(AValue: Boolean);
      procedure SetCachePath(AValue: String);
      procedure SetCenter(ACenter: TRealPoint);
      procedure SetCyclic(AValue: Boolean);
      procedure SetDownloadEngine(AValue: TMvCustomDownloadEngine);
      procedure SetHeight(AValue: integer);
      procedure SetMapProvider(AValue: String);
      procedure SetUseThreads(AValue: Boolean);
      procedure SetWidth(AValue: integer);
      procedure SetZoom(AValue: Integer); overload;
      procedure SetZoom(AValue: integer; AZoomToCursor: Boolean); overload;
      function DegreesToMapPixels(const AWin: TMapWindow; APt: TRealPoint): TPoint;
      function MapPixelsToDegrees(const AWin: TMapWindow; APoint: TPoint): TRealPoint;
      function PixelsToDegreesEPSG3395(APoint: TPoint; Zoom: Integer): TRealPoint;
      function PixelsToDegreesEPSG3857(APoint: TPoint; Zoom: Integer): TRealPoint;
      procedure CalculateWin(var AWin: TMapWindow);
      function DegreesToPixelsEPSG3395(const AWin: TMapWindow; APt: TRealPoint): TPoint;
      function DegreesToPixelsEPSG3857(const AWin: TMapWindow; APt: TRealPoint): TPoint;
      procedure Redraw(const aWin: TMapWindow; const paintOnly: Boolean = False);
      function CalculateVisibleTiles(const aWin: TMapWindow; out Area: TArea): Boolean;
      function IsCurrentWin(const aWin: TMapWindow) : boolean;
    protected
      procedure AdjustZoomCenter(var AWin: TMapWindow);
      procedure ConstraintZoom(var aWin: TMapWindow);
      function GetTileName(const Id: TTileId): String;
      procedure evDownload(Data: TObject; Job: TJob);
      procedure TileDownloaded(Data: PtrInt);
      procedure EraseBackground(const R: TRect);
      procedure DrawTileFromCache(constref ATile: TTileId; constref AWin: TMapWindow);
      procedure DrawStretchedTile(const TileId: TTileID; X, Y: Integer; TileImg: TPictureCacheItem; const R: TRect);
      Procedure DrawTile(const TileId: TTileId; X,Y: integer; TileImg: TPictureCacheItem);
      Procedure DoDrag(Sender: TDragObj);
      Procedure DoEndDrag(Sender: TDragObj);
    public
      constructor Create(aOwner: TComponent); override;
      destructor Destroy; override;

      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;
      function CrossesDateline: Boolean;
      procedure GetMapProviders(AList: TStrings);
      function MapProviderByName(AName: String): TMapProvider;
      function LatLonToScreen(APt: TRealPoint): TPoint;
      function LonLatToScreen(APt: TRealPoint): TPoint; deprecated 'Use LatLonToScreen';
      function LatLonToWorldScreen(APt: TRealPoint): TPoint;
      function LonLatToWorldScreen(APt: TRealPoint): TPoint; deprecated 'Use LatLonToWorldScreen';
      function ReadProvidersFromXML(AFileName: String; out AMsg: String): Boolean;
      procedure Redraw;
      Procedure RegisterProviders;
      function ScreenToLatLon(aPt: TPoint): TRealPoint;
      function ScreenToLonLat(aPt: TPoint): TRealPoint; deprecated 'Use ScreenToLatLon';
      procedure SetSize(aWidth, aHeight: integer);
      function WorldScreenToLatLon(aPt: TPoint): TRealPoint;
      function WorldScreenToLonLat(aPt: TPoint): TRealPoint; deprecated 'Use WorldScreenToLatLon';
      procedure WriteProvidersToXML(AFileName: String);

      procedure DblClick(Sender: TObject);
      procedure MouseDown(Sender: TObject; Button: TMouseButton;
        {%H-}Shift: TShiftState; X, Y: Integer);
      procedure MouseMove(Sender: TObject; {%H-}Shift: TShiftState;
        X, Y: Integer);
      procedure MouseUp(Sender: TObject; Button: TMouseButton;
        {%H-}Shift: TShiftState; X, Y: Integer);
      procedure MouseWheel(Sender: TObject; {%H-}Shift: TShiftState;
        WheelDelta: Integer; {%H-}MousePos: TPoint; var Handled: Boolean);
      procedure ZoomOnArea(const aArea: TRealArea);
      procedure CopyMapWindowFrom(AEngine: TMapViewerEngine);

      property MapLeft: Int64 read MapWin.X;
      property MapTop: Int64 read MapWin.Y;
      property MapProjectionType: TProjectionType read GetMapProjectionType;

      property BkColor: TFPColor read FBkColor write SetBkColor;
      property Center: TRealPoint read GetCenter write SetCenter;
      property DrawPreviewTiles : Boolean read FDrawPreviewTiles write FDrawPreviewTiles;
      property InDrag: Boolean read FInDrag;
      property DragObj: TDragObj read FDragObj;
      property CacheMaxAge: Integer read GetCacheMaxAge write SetCacheMaxAge;
      property ZoomMin: Integer read FZoomMax write FZoomMin;
      property ZoomMax: Integer read FZoomMax write FZoomMax;

    published
      property Active: Boolean read FActive write SetActive default false;
      property CacheOnDisk: Boolean read GetCacheOnDisk write SetCacheOnDisk;
      property CachePath: String read GetCachePath write SetCachePath;
      property Cyclic: Boolean read FCyclic write SetCyclic default false;
      property DownloadEngine: TMvCustomDownloadEngine
        read FDownloadEngine write SetDownloadEngine;
      property DrawTitleInGuiThread: boolean
        read FDrawTitleInGuiThread write FDrawTitleInGuiThread;
      property Height: integer read GetHeight write SetHeight;
      property JobQueue: TJobQueue read Queue;
      property MapProvider: String read GetMapProvider write SetMapProvider;
      property UseThreads: Boolean read GetUseThreads write SetUseThreads;
      property Width: integer read GetWidth write SetWidth;
      property Zoom: integer read GetZoom write SetZoom;
      property ZoomToCursor: Boolean read FZoomToCursor write FZoomToCursor default True;
      property CacheItemClass: TPictureCacheItemClass read GetCacheItemClass write SetCacheItemClass;

      property OnCenterMove: TNotifyEvent read FOnCenterMove write FOnCenterMove;
      property OnChange: TNotifyEvent Read FOnChange write FOnchange; //called when visiable area change
      property OnDrawStretchedTile: TDrawStretchedTileEvent read FOnDrawStretchedTile write FOnDrawStretchedTile;
      property OnDrawTile: TDrawTileEvent read FOnDrawTile write FOnDrawTile;
      property OnEraseBackground: TEraseBackgroundEvent read FOnEraseBackground write FOnEraseBackground;
      property OnTileDownloaded: TTileDownloadedEvent read FOnTileDownloaded write FOnTileDownloaded;
      property OnZoomChange: TNotifyEvent read FOnZoomChange write FOnZoomChange;
  end;

function RealPoint(Lat, Lon: Double): TRealPoint;
function NormalizeLon(const Lon: Double): Double; inline;

function HaversineDist(Lat1, Lon1, Lat2, Lon2, Radius: Double): Double;
function CalcGeoDistance(Lat1, Lon1, Lat2, Lon2: double;
  AUnits: TDistanceUnits = duKilometers): double;

function CalcBearing(Lat1, Lon1, Lat2, Lon2: Double): Double;
procedure CalcLatLon(const Lat1, Lon1, ADist, ABearing: Double; out Lat2, Lon2: Double);
procedure CalcMidpoint(const Lat1, Lon1, Lat2, Lon2: Double; out Lat, Lon: Double);
procedure CalcIntermedPoint(const Lat1, Lon1, Lat2, Lon2, AFrac: Double; out Lat, Lon: 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;
function TryStrDMSToDeg(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, 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                                                      // 30..32
  );

function ZoomFactor(AZoomLevel: Integer): Int64;
begin
  if (AZoomLevel >= Low(ZOOM_FACTOR)) and (AZoomLevel < High(ZOOM_FACTOR)) then
    Result := ZOOM_FACTOR[AZoomLevel]
  else
    Result := round(IntPower(2, AZoomLevel));
end;

type

  { TEnvTile }

  TEnvTile = Class(TBaseTile)
  private
    Tile: TTileId;
    Win: TMapWindow;
  public
    constructor Create(const aTile: TTileId; const aWin: TMapWindow);reintroduce;
  end;


  { TMemObj }

  TMemObj = Class
  private
    FWin: TMapWindow;
  public
    constructor Create(const aWin: TMapWindow);
  end;

  constructor TMemObj.Create(const aWin: TMapWindow);
  begin
    FWin := aWin;
  end;


{ TEnvTile }

constructor TEnvTile.Create(const aTile: TTileId; const aWin: TMapWindow);
begin
  inherited Create(aWin.MapProvider);
  Tile := aTile;
  Win := aWin;
end;


{ TMapViewerEngine }

constructor TMapViewerEngine.Create(aOwner: TComponent);
begin
  DrawTitleInGuiThread := true;
  DrawPreviewTiles := true;
  FDragObj := TDragObj.Create;
  FDragObj.OnDrag := @DoDrag;
  FDragObj.OnEndDrag := @DoEndDrag;
  Cache := TPictureCache.Create(self);
  lstProvider := TStringList.Create;
  FBkColor := colWhite;
  RegisterProviders;
  Queue := TJobQueue.Create(8);

  inherited Create(aOwner);

  FZoomMin := 1;
  FZoomMax := 19;
  FZoomToCursor := true;
  ConstraintZoom(MapWin);
  CalculateWin(mapWin);
end;

destructor TMapViewerEngine.Destroy;
begin
  Queue.CancelAllJob(Nil);
  Queue.RemoveAsyncCalls(Self);
  FreeAndNil(Queue);
  ClearMapProviders;
  FreeAndNil(FDragObj);
  FreeAndNil(lstProvider);
  FreeAndNil(Cache);
  inherited Destroy;
end;

function TMapViewerEngine.AddMapProvider(OpeName: String; ProjectionType: TProjectionType;
  Url: String; MinZoom, MaxZoom, NbSvr: integer; GetSvrStr: TGetSvrStr;
  GetXStr: TGetValStr; GetYStr: TGetValStr; GetZStr: TGetValStr): TMapProvider;
var
  idx :integer;
Begin
  idx := lstProvider.IndexOf(OpeName);
  if idx = -1 then
  begin
    Result := TMapProvider.Create(OpeName);
    lstProvider.AddObject(OpeName, Result);
  end
  else
    Result := TMapProvider(lstProvider.Objects[idx]);
  Result.AddUrl(Url, ProjectionType, NbSvr, MinZoom, MaxZoom, GetSvrStr, GetXStr, GetYStr, GetZStr);
end;

procedure TMapViewerEngine.AdjustZoomCenter(var AWin: TMapWindow);
var
  ptMouseCursor: TPoint;
  rPtAdjustedCenter: TRealPoint;
begin
  ptMouseCursor := LatLonToScreen(AWin.ZoomCenter);
  rPtAdjustedCenter := ScreenToLatLon(ptMouseCursor.Add(AWin.ZoomOffset));
  AWin.Center := rPtAdjustedCenter;
  CalculateWin(AWin);
end;

// Returns whether the view area is not covered fully with a tiles
function TMapViewerEngine.CalculateVisibleTiles(const aWin: TMapWindow; out
  Area: TArea): Boolean;
var
  MaxX, MaxY, startX, startY: int64;
  WorldMax: Int64;
begin
  Area := Default(TArea);
  Result := (aWin.X <= 0) and (aWin.Y <= 0);
  WorldMax := 1 shl AWin.Zoom - 1;
  MaxX := Int64(aWin.Width) div TILE_SIZE + 1;
  MaxY := Int64(aWin.Height) div TILE_SIZE + 1;
  if (MaxX > WorldMax) or (MaxY > WorldMax) then
  begin
    Result := False;
    MaxX := Min(WorldMax, MaxX);
    MaxY := Min(WorldMax, MaxY);
  end;
  startX := -aWin.X div TILE_SIZE;
  startY := -aWin.Y div TILE_SIZE;
  if (startX < 0) or (startY < 0) then
  begin
    startX := Max(0, startX);
    startY := Max(0, startY);
    Result := False;
  end;
  Area.Left := startX;
  Area.Right := startX + MaxX;
  Area.Top := startY;
  Area.Bottom := startY + MaxY;
end;

procedure TMapViewerEngine.CalculateWin(var AWin: TMapWindow);
var
  PixelLocation: TPoint; // review: coth: Should it use Int64?
  PType: TProjectionType;
begin
  PType := GetProjectionType(AWin);
  case PType of
    ptEPSG3857: PixelLocation := DegreesToPixelsEPSG3857(AWin, AWin.Center);
    ptEPSG3395: PixelLocation := DegreesToPixelsEPSG3395(AWin, AWin.Center);
    else PixelLocation := DegreesToPixelsEPSG3857(AWin, AWin.Center);
  end;

  AWin.X := Int64(AWin.Width div 2) - PixelLocation.x;
  AWin.Y := Int64(AWin.Height div 2) - PixelLocation.y;
end;

procedure TMapViewerEngine.CancelCurrentDrawing;
var
  Jobs: TJobArray;
begin
  Jobs := Queue.CancelAllJob(self);
  Queue.WaitForTerminate(Jobs);
end;

procedure TMapViewerEngine.ClearMapProviders;
var
  i: Integer;
begin
  for i:=0 to lstProvider.Count-1 do
    TObject(lstProvider.Objects[i]).Free;
  lstProvider.Clear;
end;

procedure TMapViewerEngine.ConstraintZoom(var aWin: TMapWindow);
var
  zMin, zMax: integer;
begin
  if Assigned(aWin.MapProvider) then
  begin
    aWin.MapProvider.GetZoomInfos(zMin, zMax);
    if aWin.Zoom < zMin then
      aWin.Zoom := zMin;
    if aWin.Zoom > zMax then
      aWin.Zoom := zMax;
  end;
end;

{ Returns true when the visible window crosses the date line, i.e. the longitudes
  at the left of the window are > 0, and those at the right are < 0. }
function TMapViewerEngine.CrossesDateline: Boolean;
var
  visArea: TRealArea;
  mapWidth: Int64;
begin
  // Catch the case, that the screen is wider than the whole world
  mapWidth := ZoomFactor(MapWin.Zoom) * TILE_SIZE;
  Result := (MapWin.Width > mapWidth);
  if not Result then
  begin
    visArea.TopLeft := ScreenToLatLon(Point(0, 0));
    visArea.BottomRight := ScreenToLatLon(Point(Width, Height));
    Result := (visArea.TopLeft.Lon > 0) and (visArea.BottomRight.Lon < 0);
  end;
end;

procedure TMapViewerEngine.DblClick(Sender: TObject);
var
  pt: TPoint;
begin
  pt.X := FDragObj.MouseX;
  pt.Y := FDragObj.MouseY;
  SetCenter(ScreenToLatLon(pt));
end;

procedure TMapViewerEngine.DoDrag(Sender: TDragObj);
begin
  if Sender.DragSrc = self then
  begin
    FInDrag := True;
    MoveMapCenter(Sender);
  end;
end;

procedure TMapViewerEngine.DoEndDrag(Sender: TDragObj);
begin
  FInDrag := False;
end;

procedure TMapViewerEngine.DrawStretchedTile(const TileId: TTileID; X,
  Y: Integer; TileImg: TPictureCacheItem; const R: TRect);
begin
  if Assigned(FOnDrawStretchedTile) then
    FOnDrawStretchedTile(TileId, X, Y, TileImg, R);
end;

procedure TMapViewerEngine.DrawTile(const TileId: TTileId; X, Y: integer;
  TileImg: TPictureCacheItem);
begin
  if Assigned(FOnDrawTile) then
    FOnDrawTile(TileId, X, Y, TileImg);
end;

procedure TMapViewerEngine.evDownload(Data: TObject; Job: TJob);
var
  Id: TTileId;
  Url: String;
  Env: TEnvTile;
  MapO: TMapProvider;
  lStream: TMemoryStream;
begin
  Env := TEnvTile(Data);
  Id := Env.Tile;
  MapO := Env.Win.MapProvider;
  if Assigned(MapO) and Assigned(Cache) then
  begin
    if not Cache.InCache(MapO, Id) then
    begin
      if Assigned(FDownloadEngine) then
      begin
        Url := MapO.GetUrlForTile(Id);
        if Url <> '' then
        begin
          lStream := TMemoryStream.Create;
          try
            try
              FDownloadEngine.DownloadFile(Url, lStream);
              if Assigned(Cache) then Cache.Add(MapO, Id, lStream);
            except
            end;
          finally
            FreeAndNil(lStream);
          end;
        end;
      end;
    end;
  end;

  if Job.Cancelled then
    Exit;

  if DrawTitleInGuiThread then
    Queue.QueueAsyncCall(@TileDownloaded, PtrInt(Env))
  else
    TileDownloaded(PtrInt(Env));
end;

function TMapViewerEngine.GetCacheOnDisk: Boolean;
begin
  Result := Cache.UseDisk;
end;

function TMapViewerEngine.GetCacheItemClass: TPictureCacheItemClass;
begin
  Result := Cache.CacheItemClass;
end;

function TMapViewerEngine.GetCacheMaxAge: Integer;
begin
  Result := Cache.MaxAge;
end;

function TMapViewerEngine.GetCachePath: String;
begin
  Result := Cache.BasePath;
end;

function TMapViewerEngine.GetCenter: TRealPoint;
begin
  Result := MapWin.Center;
end;

function TMapViewerEngine.GetHeight: integer;
begin
  Result := MapWin.Height
end;

class function TMapViewerEngine.GetProjectionType(const aWin: TMapWindow
  ): TProjectionType;
begin
  if Assigned(AWin.MapProvider)
    then Result := AWin.MapProvider.ProjectionType
    else Result := ptEPSG3857; // Default
end;

function TMapViewerEngine.GetMapProjectionType: TProjectionType;
begin
  Result := GetProjectionType(MapWin);
end;

function TMapViewerEngine.GetMapProvider: String;
begin
  if Assigned(MapWin.MapProvider) then
    Result := MapWin.MapProvider.Name
  else
    Result := '';
end;

procedure TMapViewerEngine.GetMapProviders(AList: TStrings);
begin
  AList.Assign(lstProvider);
end;

function TMapViewerEngine.MapProviderByName(AName: String): TMapProvider;
var
  I: Integer;
begin
  I := lstProvider.IndexOf(AName);
  if I <> -1
    then Result := TMapProvider(lstProvider.Objects[I])
    else Result := Nil;
end;

function TMapViewerEngine.GetTileName(const Id: TTileId): String;
begin
  Result := IntToStr(Id.X) + '.' + IntToStr(Id.Y) + '.' + IntToStr(Id.Z);
end;

function TMapViewerEngine.GetUseThreads: Boolean;
begin
  Result := Queue.UseThreads;
end;

function TMapViewerEngine.GetWidth: integer;
begin
  Result := MapWin.Width;
end;

function TMapViewerEngine.GetZoom: integer;
begin
  Result := MapWin.Zoom;
end;

function TMapViewerEngine.IsCurrentWin(const aWin: TMapWindow): boolean;
begin
  Result := (aWin.Zoom = MapWin.Zoom) and
            (aWin.Center.Lat = MapWin.Center.Lat) and
            (aWin.Center.Lon = MapWin.Center.Lon) and
            (aWin.Width = MapWin.Width) and
            (aWin.Height = MapWin.Height);
end;

function TMapViewerEngine.IsValidTile(const aWin: TMapWindow;
  const aTile: TTIleId): boolean;
var
  tiles: int64;
begin
  tiles := 1 shl aWin.Zoom;
  Result := (aTile.X >= 0) and (aTile.X <= tiles-1) and
            (aTile.Y >= 0) and (aTile.Y <= tiles-1);
end;

function TMapViewerEngine.DegreesToMapPixels(const AWin: TMapWindow;
  APt: TRealPoint): TPoint;
var
  pixelLocation: TPoint;
  mapWidth: Int64;
  PType: TProjectionType;
begin
  PType := GetProjectionType(AWin);
  case PType of
    ptEPSG3395: pixelLocation := DegreesToPixelsEPSG3395(AWin, APt);
    ptEPSG3857: pixelLocation := DegreesToPixelsEPSG3857(AWin, APt);
    else pixelLocation := DegreesToPixelsEPSG3857(AWin, APt);
  end;
  Result.X := pixelLocation.x + AWin.X;
  if FCyclic and CrossesDateline then
  begin
    mapWidth := ZoomFactor(AWin.Zoom) * TILE_SIZE;
    while (Result.X < 0) do
      Result.X := Result.X + mapWidth;
    while (Result.X > AWin.Width) do
      Result.X := Result.X - mapWidth;
  end;
  Result.Y := pixelLocation.y + AWin.Y;
end;

// review: coth: Should it use Int64?
function TMapViewerEngine.DegreesToPixelsEPSG3857(const AWin: TMapWindow;
  APt: TRealPoint): TPoint;
const
  MIN_LATITUDE = -85.05112878;
  MAX_LATITUDE = 85.05112878;
  MIN_LONGITUDE = -180;
  MAX_LONGITUDE = 180;
  TWO_PI = 2.0 * pi;
var
  factor, 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(APt.Lat, MIN_LATITUDE, MAX_LATITUDE);
  pt.Lon := Math.EnsureRange(APt.Lon, MIN_LONGITUDE, MAX_LONGITUDE);

  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);
end;

// review: coth: Should it use Int64?
function TMapViewerEngine.DegreesToPixelsEPSG3395(const AWin: TMapWindow;
  APt: 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;
  zoomfac: Extended;  // 2**Z
begin
  // https://epsg.io/3395
  // https://pubs.usgs.gov/pp/1395/report.pdf, page 44
  pt.Lat := Math.EnsureRange(APt.Lat, MIN_LATITUDE, MAX_LATITUDE);
  pt.Lon := Math.EnsureRange(APt.Lon, MIN_LONGITUDE, MAX_LONGITUDE);

  Z := 23 - AWin.Zoom;
  zoomfac := ZoomFactor(Z);
  cfmpx := IntPower(2, 31);
  cfmpm := cfmpx / EARTH_CIRCUMFERENCE;
  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 / zoomfac;

  Result.x := Round(px);
  Result.y := Round(py);
end;

function TMapViewerEngine.LatLonToScreen(APt: TRealPoint): TPoint;
begin
  Result := DegreesToMapPixels(MapWin, APt);
end;

function TMapViewerEngine.LonLatToScreen(APt: TRealPoint): TPoint;
Begin
  Result := DegreesToMapPixels(MapWin, APt);
end;

function TMapViewerEngine.LatLonToWorldScreen(APt: TRealPoint): TPoint;
begin
  Result := LatLonToScreen(APt);
  Result.X := Result.X + MapWin.X;
  Result.Y := Result.Y + MapWin.Y;
end;

function TMapViewerEngine.LonLatToWorldScreen(APt: TRealPoint): TPoint;
begin
  Result := LatLonToScreen(APt);
end;

function TMapViewerEngine.MapPixelsToDegrees(const AWin: TMapWindow;
  APoint: TPoint): TRealPoint;
var
  mapWidth: Int64;
  mPoint : TPoint;
  PType: TProjectionType;
begin
  mapWidth := round(ZoomFactor(AWin.Zoom)) * TILE_SIZE;

  if FCyclic then
  begin
    mPoint.X := (APoint.X - AWin.X) mod mapWidth;
    while mPoint.X < 0 do
      mPoint.X := mPoint.X + mapWidth;
    while mPoint.X >= mapWidth do
      mPoint.X := mPoint.X - mapWidth;
  end else
    mPoint.X := EnsureRange(APoint.X - AWin.X, 0, mapWidth);
  mPoint.Y := EnsureRange(APoint.Y - AWin.Y, 0, mapWidth);

  PType := GetProjectionType(AWin);
  case PType of
    ptEPSG3857: Result := PixelsToDegreesEPSG3857(mPoint, AWin.Zoom);
    ptEPSG3395: Result := PixelsToDegreesEPSG3395(mPoint, AWin.Zoom);
    else        Result := PixelsToDegreesEPSG3857(mPoint, AWin.Zoom);
  end;
end;

function TMapViewerEngine.PixelsToDegreesEPSG3857(APoint: TPoint; Zoom: Integer): TRealPoint;
const
  MIN_LATITUDE = -85.05112878;
  MAX_LATITUDE = 85.05112878;
  MIN_LONGITUDE = -180;
  MAX_LONGITUDE = 180;
var
  zoomfac: Int64;
begin
  // https://epsg.io/3857
  // https://pubs.usgs.gov/pp/1395/report.pdf, page 41

  // 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) );
  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);
end;

function TMapViewerEngine.PixelsToDegreesEPSG3395(APoint: TPoint; Zoom: Integer
  ): TRealPoint;

  function PhiIteration(y, phi: Extended): Extended;
  var
    t: Extended;
    sin_phi: Extended;
    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;

const
  MIN_LATITUDE = -80;
  MAX_LATITUDE = 84;
  MIN_LONGITUDE = -180;
  MAX_LONGITUDE = 180;
  EPS = 1e-8;
var
  LonRad, LatRad: Extended;
  WorldSize: Int64;
  Cpm: Extended;
  Z: Integer;
  t, phi:  Extended;
  zoomFac: Int64;
  i: Integer;
begin
  // https://epsg.io/3395
  // https://pubs.usgs.gov/pp/1395/report.pdf, page 44

  Z := 23 - Zoom;
  zoomFac := ZoomFactor(Z);
  WorldSize := ZoomFactor(31);
  Cpm :=  WorldSize / EARTH_CIRCUMFERENCE;

  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));

  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;
  Shift: TShiftState; X, Y: Integer);
begin
  FDragObj.MouseDown(self,X,Y);
end;

procedure TMapViewerEngine.MouseMove(Sender: TObject; Shift: TShiftState;
  X, Y: Integer);
begin
  FDragObj.MouseMove(X,Y);
end;

procedure TMapViewerEngine.MouseUp(Sender: TObject; Button: TMouseButton;
  Shift: TShiftState; X, Y: Integer);
begin
  FDragObj.MouseUp(X,Y);
end;

procedure TMapViewerEngine.MouseWheel(Sender: TObject;
  Shift: TShiftState; WheelDelta: Integer; MousePos: TPoint;
  var Handled: Boolean);
var
  Val: Integer;
  nZoom: integer;
  bZoomToCursor: Boolean;
begin
  bZoomToCursor := False;
  Val := 0;
  if WheelDelta > 0 then
    Val := 1;
  if WheelDelta < 0 then
    Val := -1;
  nZoom := Zoom + Val;
  if (nZoom >= FZoomMin) and (nZoom <= FZoomMax) then
  begin
    if ZoomToCursor then
    begin
      MapWin.ZoomCenter := ScreenToLatLon(MousePos);
      MapWin.ZoomOffset := LatLonToScreen(Center).Subtract(MousePos);
      bZoomToCursor := True;
    end;
    SetZoom(nZoom, bZoomToCursor);
  end;
  Handled := true;
end;

procedure TMapViewerEngine.MoveMapCenter(Sender: TDragObj);
var
  old: TMemObj;
  nCenter: TRealPoint;
  aPt: TPoint;
Begin
  if Sender.LnkObj = nil then
    Sender.LnkObj := TMemObj.Create(MapWin);
  old := TMemObj(Sender.LnkObj);
  aPt.X := old.FWin.Width DIV 2-Sender.OfsX;
  aPt.Y := old.FWin.Height DIV 2-Sender.OfsY;
  nCenter := MapPixelsToDegrees(old.FWin,aPt);
  SetCenter(nCenter);
end;

function TMapViewerEngine.ReadProvidersFromXML(AFileName: String;
  out AMsg: String): Boolean;

  function GetSvrStr(AName: String): TGetSvrStr;
  var
    lcName: String;
  begin
    lcName := LowerCase(AName);
    if lcName = LowerCase(SVR_LETTER) then
      Result := @GetSvrLetter
    else if lcName = LowerCase(SVR_BASE1) then
      Result := @GetSvrBase1
    else
      Result := nil;
  end;

  function GetValStr(AName: String): TGetValStr;
  var
    lcName: String;
  begin
    lcName := Lowercase(AName);
    if lcName = LowerCase(STR_QUADKEY) then
      Result := @GetStrQuadKey
    else if lcName = LowerCase(STR_YAHOOY) then
      Result := @GetStrYahooY
    else if lcName = LowerCase(STR_YAHOOZ) then
      Result := @GetStrYahooZ
    else
      Result := nil;
  end;

  function GetAttrValue(ANode: TDOMNode; AttrName: String): String;
  var
    node: TDOMNode;
  begin
    Result := '';
    if ANode.HasAttributes then begin
      node := ANode.Attributes.GetNamedItem(AttrName);
      if Assigned(node) then Result := node.NodeValue;
    end;
  end;

var
  stream: TFileStream;
  doc: TXMLDocument = nil;
  node, layerNode: TDOMNode;
  providerName: String;
  projectionType: TProjectionType;
  url: String;
  minZoom: Integer;
  maxZoom: Integer;
  svrCount: Integer;
  s: String;
  svrProc: String;
  xProc: String;
  yProc: String;
  zProc: String;
  first: Boolean;
begin
  Result := false;
  AMsg := '';
  stream := TFileStream.Create(AFileName, fmOpenread or fmShareDenyWrite);
  try
    ReadXMLFile(doc, stream, [xrfAllowSpecialCharsInAttributeValue, xrfAllowLowerThanInAttributeValue]);
    node := doc.FindNode('map_providers');
    if node = nil then begin
      AMsg := 'No map providers in file.';
      exit;
    end;

    first := true;
    node := node.FirstChild;
    while node <> nil do begin
      providerName := GetAttrValue(node, 'name');
      layerNode := node.FirstChild;
      while layerNode <> nil do begin
        url := GetAttrValue(layerNode, 'url');
        if url = '' then
          continue;
        s := GetAttrValue(layerNode, 'minZom');
        if s = '' then minZoom := 0
          else minZoom := StrToInt(s);
        s := GetAttrValue(layerNode, 'maxZoom');
        if s = '' then maxzoom := 9
          else maxZoom := StrToInt(s);
        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');
        zProc := GetAttrValue(layerNode, 'zProc');
        layerNode := layerNode.NextSibling;
      end;
      if first then begin
        ClearMapProviders;
        first := false;
      end;
      AddMapProvider(providerName, projectionType,
        url, minZoom, maxZoom, svrCount,
        GetSvrStr(svrProc), GetValStr(xProc), GetValStr(yProc), GetValStr(zProc)
      );
      node := node.NextSibling;
    end;
    Result := true;
  finally
    stream.Free;
    doc.Free;
  end;
end;

procedure TMapViewerEngine.Redraw;
begin
  Redraw(MapWin, FDragObj.InDrag);
end;

procedure TMapViewerEngine.Redraw(const aWin: TMapWindow;
  const paintOnly: Boolean);
var
  TilesVis: TArea;
  x, y, px, py: Integer;
  iTile, numTiles, XShift: Integer;
  Tiles: TTileIdArray = nil;
  tile: TTileID;
  previewDrawn: Boolean;
  previewImg: TPictureCacheItem;
  R: TRect;

  procedure AddJob;
  var
    lTile: TEnvTile;
    lJob: TEventJob;
  begin
    lTile:=TEnvTile.Create(Tiles[iTile], aWin);
    lJob := TEventJob.Create(@evDownload, lTile, False,                                    // owns data
      GetTileName(Tiles[iTile]));
    if not Queue.AddUniqueJob(lJob, Self) then
    begin
      FreeAndNil(lJob);
      FreeAndNil(lTile);
    end;
  end;

  procedure EraseAround;
  var
    T, L, B, R: Integer;
  begin
    T := -AWin.Y div TILE_SIZE - Max(0, Sign(AWin.Y));
    B := T + AWin.Height div TILE_SIZE + 1;
    L := -AWin.X div TILE_SIZE - Max(0, Sign(AWin.X));
    R := L + AWin.Width div TILE_SIZE + 1;
    if T < TilesVis.top then  // Erase above top
      EraseBackground(Rect(0, 0, AWin.Width, AWin.Y + TilesVis.top * TILE_SIZE));
    if L < TilesVis.left then // Erase on the left
      EraseBackground(Rect(0, AWin.Y + TilesVis.top * TILE_SIZE,
        AWin.X + TilesVis.left * TILE_SIZE,
        AWin.Y + (TilesVis.bottom + 1) * TILE_SIZE));
    if R > TilesVis.right then // Erase on the right
      EraseBackground(Rect(AWin.X + (TilesVis.right + 1) * TILE_SIZE,
        AWin.Y + TilesVis.top * TILE_SIZE, AWin.Width,
        AWin.Y + (TilesVis.bottom + 1) * TILE_SIZE));
    if B > TilesVis.bottom then // Erase below
      EraseBackground(Rect(0, AWin.Y + (TilesVis.bottom + 1) * TILE_SIZE,
        AWin.Width, AWin.Height));
  end;

begin
  if not(Active) then
    Exit;
  if not Assigned(AWin.MapProvider) then
  begin
    EraseBackground(Rect(0, 0, AWin.Width, AWin.Height));
    Exit;
  end;

  if not CalculateVisibleTiles(AWin, TilesVis) then
    EraseAround;
  SetLength(Tiles, (TilesVis.Bottom - TilesVis.Top + 1) * (TilesVis.Right - TilesVis.Left + 1));
  iTile := Low(Tiles);
  numTiles := 1 shl AWin.Zoom;
  XShift := IfThen(aWin.X > 0, numTiles - aWin.X div TILE_SIZE - 1, 0);
  for y := TilesVis.Top to TilesVis.Bottom do
    for X := TilesVis.Left to TilesVis.Right do
    begin
      if FCyclic then
      begin // 0,1,2,3,4,5 --> 15,16,>0<,1,2,3
        Tiles[iTile].X := (X + XShift) mod numTiles;
        if Tiles[iTile].X < 0 then  //
          Tiles[iTile].X := Tiles[iTile].X + numTiles;
      end
      else
        Tiles[iTile].X := X;
      Tiles[iTile].Y := Y;
      Tiles[iTile].Z := AWin.Zoom;

      if Cache.InCache(AWin.MapProvider, Tiles[iTile]) then
        DrawTileFromCache(Tiles[iTile], AWin)
      else
      // Avoid tiling artefacts when a tile does not exist (lowest zoom) or
      // is not valid
      begin
        previewdrawn := False;
        py := AWin.Y + Y * TILE_SIZE;
        px := AWin.X + X * TILE_SIZE;
        if FDrawPreviewTiles then
        begin
          if IsValidTile(AWin, Tiles[iTile]) then  // Invalid tiles probably will not be found in the cache
          begin
            tile := Tiles[iTile];
            if Cache.GetPreviewFromCache(AWin.MapProvider, tile, R) then
            begin
              Cache.GetFromCache(AWin.MapProvider, tile, previewImg);
              DrawStretchedTile(Tiles[iTile], px, py, previewImg, R);
              previewDrawn := true;
            end;
          end;
        end;
        if not previewDrawn then
          DrawTile(Tiles[iTile], px, py, nil);  // Draw blank tile if preview cannot be generated

        if not paintOnly and IsValidTile(AWin, Tiles[iTile]) then
        begin
          AddJob;
          inc(iTile);
        end;
      end;
    end;
  SetLength(Tiles, iTile);
  if Length(Tiles) < 1 then
    Cache.CheckCacheSize(Nil)
end;

procedure TMapViewerEngine.RegisterProviders;
var
  HERE1, HERE2: String;
begin
  {$I mvengine_mapreg.inc}
  MapWin.MapProvider := Nil; // Undo the OSM Mapnik default selection
end;

function TMapViewerEngine.ScreenToLatLon(aPt: TPoint): TRealPoint;
begin
  Result := MapPixelsToDegrees(MapWin, aPt);
end;

function TMapViewerEngine.ScreenToLonLat(aPt: TPoint): TRealPoint;
begin
  Result := MapPixelsToDegrees(MapWin, aPt);
end;

procedure TMapViewerEngine.SetActive(AValue: boolean);
begin
  if FActive = AValue then Exit;
  FActive := AValue;

  // One problem at designtime is that the control is creating cache directories
  // at unexpected places when Cache.BasePath is relative
  if csDesigning in ComponentState then
    exit;

  if not FActive then
    Queue.CancelAllJob(self)
  else begin
    if Cache.UseDisk then ForceDirectories(Cache.BasePath);
    Redraw(MapWin);
  end;
end;

procedure TMapViewerEngine.SetBkColor(AValue: TFPColor);
begin
  if FBkColor = AValue then Exit;
  FBkColor := AValue;
  Redraw(MapWin);
end;

procedure TMapViewerEngine.SetCacheItemClass(AValue: TPictureCacheItemClass);
begin
  if Cache.CacheItemClass = AValue then
    Exit;
  Cache.CacheItemClass := AValue;
end;

procedure TMapViewerEngine.SetCacheMaxAge(AValue: Integer);
begin
  if Cache.MaxAge = AValue then
    Exit;
  Cache.MaxAge := AValue;
end;

procedure TMapViewerEngine.SetCacheOnDisk(AValue: Boolean);
begin
  if Cache.UseDisk = AValue then Exit;
  Cache.UseDisk := AValue;
end;

procedure TMapViewerEngine.SetCachePath(AValue: String);
begin
  Cache.BasePath := aValue;
end;

procedure TMapViewerEngine.SetCenter(ACenter: TRealPoint);
begin
  if (MapWin.Center.Lon <> aCenter.Lon) or (MapWin.Center.Lat <> aCenter.Lat) then
  begin
    Mapwin.Center := aCenter;
    CalculateWin(MapWin);
    if Assigned(OnCenterMove) then
      OnCenterMove(Self);
    if Assigned(OnChange) then
      OnChange(Self);
  end;
end;

procedure TMapViewerEngine.SetCyclic(AValue: Boolean);
begin
  if FCyclic = AValue then exit;
  FCyclic := AValue;
  if Assigned(MapWin.MapProvider) and CrossesDateLine then
    Redraw;
end;

procedure TMapViewerEngine.SetDownloadEngine(AValue: TMvCustomDownloadEngine);
begin
  if FDownloadEngine = AValue then Exit;
  FDownloadEngine := AValue;
  if Assigned(FDownloadEngine) then
    FDownloadEngine.FreeNotification(self);
end;

procedure TMapViewerEngine.SetHeight(AValue: integer);
begin
  if MapWin.Height = AValue then Exit;
  MapWin.Height := AValue;
  CalculateWin(MapWin);
  Redraw(MapWin);
end;

procedure TMapViewerEngine.SetMapProvider(AValue: String);
var
  Provider: TMapProvider;
begin
  Provider := MapProviderByName(AValue);
  if not ((aValue = '') or Assigned(Provider)) then
    raise Exception.Create('Unknow Provider: ' + aValue);
  if Assigned(MapWin.MapProvider) and (MapWin.MapProvider.Name = AValue) then Exit;
  MapWin.MapProvider := Provider;
  if Assigned(Provider) then
  begin
    ConstraintZoom(MapWin);
    CalculateWin(MapWin);
    Redraw(MapWin);
  end;
end;

procedure TMapViewerEngine.SetSize(aWidth, aHeight: integer);
begin
  if (MapWin.Width = aWidth) and (MapWin.Height = aHeight) then Exit;
  CancelCurrentDrawing;
  MapWin.Width := aWidth;
  MapWin.Height := aHeight;
  CalculateWin(MapWin);
  Redraw(MapWin);
  if Assigned(OnChange) then
    OnChange(Self);
end;

procedure TMapViewerEngine.SetUseThreads(AValue: Boolean);
begin
  if Queue.UseThreads = AValue then Exit;
  Queue.UseThreads := AValue;
  Cache.UseThreads := AValue;
end;

procedure TMapViewerEngine.SetWidth(AValue: integer);
begin
  if MapWin.Width = AValue then Exit;
  MapWin.Width := AValue;
  CalculateWin(MapWin);
  Redraw(MapWin);
end;

procedure TMapViewerEngine.SetZoom(AValue: Integer);
begin
  SetZoom(AValue, false);
end;

procedure TMapViewerEngine.SetZoom(AValue: integer; AZoomToCursor: Boolean);
begin
  if MapWin.Zoom = AValue then Exit;
  MapWin.Zoom := AValue;
  ConstraintZoom(MapWin);
  CalculateWin(MapWin);
  if AZoomToCursor then
    AdjustZoomCenter(MapWin);
  Redraw(MapWin, True);
  if Assigned(OnZoomChange) then
    OnZoomChange(Self);
  if Assigned(OnChange) then
    OnChange(Self);
end;

procedure TMapViewerEngine.TileDownloaded(Data: PtrInt);
var
  EnvTile: TEnvTile;
begin
  EnvTile := TEnvTile(Data);
  try
    if Assigned(FOnTileDownloaded) then
      FOnTileDownloaded(EnvTile.Tile);
  finally
    FreeAndNil(EnvTile);
  end;
end;

procedure TMapViewerEngine.EraseBackground(const R: TRect);
begin
  if Assigned(FOnEraseBackground) then
      FOnEraseBackground(R);
end;

procedure TMapViewerEngine.DrawTileFromCache(constref ATile: TTileId; constref
  AWin: TMapWindow);
var
  img: TPictureCacheItem;
  X, Y: integer;
  worldWidth : Integer;
  numTiles : Integer;
  baseX : Integer;
begin
  if IsCurrentWin(AWin)then
  begin
     Cache.GetFromCache(AWin.MapProvider, ATile, img);
     Y := AWin.Y + ATile.Y * TILE_SIZE; // begin of Y
     if Cyclic then
     begin
       baseX := AWin.X + ATile.X * TILE_SIZE; // begin of X
       numTiles := 1 shl AWin.Zoom;
       worldWidth := numTiles * TILE_SIZE;
       // From the center to the left (western) hemisphere
       X := baseX;
       while (X+TILE_SIZE >= 0) do
       begin
         DrawTile(ATile, X, Y, img);
         X := X - worldWidth;
       end;
       // From the center to the right (eastern) hemisphere
       X := baseX + worldWidth;
       while ((X-TILE_SIZE) <= AWin.Width) do
       begin
         DrawTile(ATile, X, Y, img);
         X := X + worldWidth;
       end;
     end else
     begin
       X := AWin.X + ATile.X * TILE_SIZE; // begin of X
       DrawTile(ATile, X, Y, img);
     end;
  end;
end;

function TMapViewerengine.WorldScreenToLatLon(aPt: TPoint): TRealPoint;
begin
  aPt.X := aPt.X - MapWin.X;
  aPt.Y := aPt.Y - MapWin.Y;
  Result := ScreenToLatLon(aPt);
end;

function TMapViewerEngine.WorldScreenToLonLat(aPt: TPoint): TRealPoint;
begin
  Result := WorldScreenToLatLon(aPt);
end;

procedure TMapViewerEngine.WriteProvidersToXML(AFileName: String);
var
  doc: TXMLDocument;
  root: TDOMNode;
  i: Integer;
  prov: TMapProvider;
begin
  doc := TXMLDocument.Create;
  try
    root := doc.CreateElement('map_providers');
    doc.AppendChild(root);
    for i := 0 to lstProvider.Count - 1 do begin
      prov := TMapProvider(lstProvider.Objects[i]);
      prov.ToXML(doc, root);
    end;
    WriteXMLFile(doc, AFileName);
  finally
    doc.Free;
  end;
end;

procedure TMapViewerEngine.ZoomOnArea(const aArea: TRealArea);
var
  tmpWin: TMapWindow;
  visArea: TRealArea;
  TopLeft, BottomRight: TPoint;
begin
  tmpWin := MapWin;
  tmpWin.Center.Lon := (aArea.TopLeft.Lon + aArea.BottomRight.Lon) / 2;
  tmpWin.Center.Lat := (aArea.TopLeft.Lat + aArea.BottomRight.Lat) / 2;
  tmpWin.Zoom := 18;
  TopLeft.X := 0;
  TopLeft.Y := 0;
  BottomRight.X := tmpWin.Width;
  BottomRight.Y := tmpWin.Height;
  Repeat
    CalculateWin(tmpWin);
    visArea.TopLeft := MapPixelsToDegrees(tmpWin, TopLeft);
    visArea.BottomRight := MapPixelsToDegrees(tmpWin, BottomRight);
    if AreaInsideArea(aArea, visArea) then
      break;
    dec(tmpWin.Zoom);
  until (tmpWin.Zoom = 2);
  MapWin := tmpWin;
  Redraw(MapWin);
end;

procedure TMapViewerEngine.CopyMapWindowFrom(AEngine: TMapViewerEngine);
var
  MP: TMapProvider;
begin
  MP := MapWin.MapProvider;
  MapWin := AEngine.MapWin;
  MapWin.MapProvider := MP; // Restore the old Map Provider
end;

//------------------------------------------------------------------------------

function RealPoint(Lat, Lon: Double): TRealPoint;
begin
  Result.Lon := Lon;
  Result.Lat := Lat;
end;

procedure SplitGps(AValue: Double; out ADegs, AMins: Double);
begin
  AValue := abs(AValue);
  AMins := frac(AValue) * 60;
  if abs(AMins - 60) < 1E-3 then
  begin
    AMins := 0;
    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);
begin
  SplitGps(AValue, ADegs, AMins);
  ASecs := frac(AMins) * 60;
  AMins := trunc(AMins);
  if abs(ASecs - 60) < 1E-3 then
  begin
    ASecs := 0;
    AMins := AMins + 1;
    if abs(AMins - 60) < 1e-3 then
    begin
      AMins := 0;
      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'' %.*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), AFormatSettings)
  else
    Result := Format('%.6f°',[abs(ALatitude)], AFormatSettings);
  if ALatitude > 0 then
    Result := Result + ' N'
  else
  if ALatitude < 0 then
    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), AFormatSettings)
  else
    Result := Format('%.6f°', [abs(ALongitude)], AFormatSettings);
  if ALongitude > 0 then
    Result := Result + ' E'
  else if ALongitude < 0 then
    Result := Result + ' W';
end;

{ Combines up to three parts of a GPS coordinate string (degrees, minutes, seconds)
  to a floating-point degree value. The parts are separated by non-numeric
  characters:

  three parts ---> d m s ---> d and m must be integer, s can be float
  two parts   ---> d m   ---> d must be integer, s can be float
  one part    ---> d     ---> d can be float

  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!

  Hemisphere suffixes ('N', 'S', 'E', 'W') are supported at the end of the input string.
}
function TryStrToGps(const AValue: String; out ADeg: Double): Boolean;
const
  NUMERIC_CHARS = ['0'..'9', '.', ',', '-', '+'];
var
  mins, secs: Double;
  i, j, len: Integer;
  n: Integer;
  s: String = '';
  res: Integer;
  sgn: Double;
begin
  Result := false;

  ADeg := NaN;
  mins := 0;
  secs := 0;

  if AValue = '' then
    exit;

  len := Length(AValue);
  i := len;
  while (i >= 1) and (AValue[i] = ' ') do dec(i);
  sgn := 1.0;
  if (AValue[i] in ['S', 's', 'W', 'w']) then sgn := -1;

  // skip leading non-numeric characters
  i := 1;
  while (i <= len) and not (AValue[i] in NUMERIC_CHARS) do
    inc(i);

  // extract first value: degrees
  SetLength(s, len);
  j := 1;
  n := 0;
  while (i <= len) and (AValue[i] in NUMERIC_CHARS) do begin
    if AValue[i] = ',' then s[j] := '.' else s[j] := AValue[i];
    inc(i);
    inc(j);
    inc(n);
  end;
  if n > 0 then begin
    SetLength(s, n);
    val(s, ADeg, res);
    if res <> 0 then
      exit;
  end;

  // skip non-numeric characters between degrees and minutes
  while (i <= len) and not (AValue[i] in NUMERIC_CHARS) do
    inc(i);

  // extract second value: minutes
  SetLength(s, len);
  j := 1;
  n := 0;
  while (i <= len) and (AValue[i] in NUMERIC_CHARS) do begin
    if AValue[i] = ',' then s[j] := '.' else s[j] := AValue[i];
    inc(i);
    inc(j);
    inc(n);
  end;
  if n > 0 then begin
    SetLength(s, n);
    val(s, mins, res);
    if (res <> 0) or (mins < 0) then
      exit;
  end;

  // skip non-numeric characters between minutes and seconds
  while (i <= len) and not (AValue[i] in NUMERIC_CHARS) do
    inc(i);

  // extract third value: seconds
  SetLength(s, len);
  j := 1;
  n := 0;
  while (i <= len) and (AValue[i] in NUMERIC_CHARS) do begin
    if AValue[i] = ',' then s[j] := '.' else s[j] := AValue[i];
    inc(i);
    inc(j);
    inc(n);
  end;
  if n > 0 then begin
    SetLength(s, n);
    val(s, secs, res);
    if (res <> 0) or (secs < 0) then
      exit;
  end;

  // If the string contains seconds then minutes and deegrees must be integers
  if (secs <> 0) and ((frac(ADeg) > 0) or (frac(mins) > 0)) then
    exit;
  // If the string does not contain seconds then degrees must be integer.
  if (secs = 0) and (mins <> 0) and (frac(ADeg) > 0) then
    exit;

  // If the string contains minutes, but no seconds, then the degrees must be integer.
  Result := (mins >= 0) and (mins < 60) and (secs >= 0) and (secs < 60);

  // A similar check should be made for the degrees range, but since this is
  // different for latitude and longitude the check is skipped here.
  if Result then
    ADeg := sgn * (abs(ADeg) + mins / 60 + secs / 3600);
end;

{ Convert Lat/Lon string to degrees.

  Recognized formats:

    Degrees, Minutes and Seconds:
      DDD°MM'SS.S"
      32°18'23.1"N 122°36'52.5"W

    Degrees and Decimal Minutes:
      DDD° MM.MMM'
      32°18.385'N 122°36.875'W

    Decimal Degrees:
      DDD.DDDDD°
      32.30642°N 122.61458°W
      +32.30642 -122.61458
}
function TryStrDMSToDeg(const AValue: String; out ADeg: Double): Boolean;
const
  NUMERIC_CHARS = ['0'..'9', '.', ',', '-', '+'];
  WS_CHARS = [' '];
var
  I, Len, N: Integer;
  S: String;
  D, Minutes, Seconds: Double;
  R: Word;
  Last, Neg: Boolean;

  function EOL: Boolean; inline;
  begin
    Result := Len < I;
  end;

  procedure SkipWS; inline;
  begin
    while not EOL and (AValue[I] in WS_CHARS) do
      Inc(I);
  end;

  function NextNum: String;
  begin
    Result := '';
    SkipWS;
    while not EOL and (AValue[I] in NUMERIC_CHARS) do
    begin
      if AValue[I] = ','
        then Result := Result + '.'
        else Result := Result + AValue[I];
      Inc(I);
    end;
  end;

begin
  Result := False;
  ADeg := NaN;
  Len := Length(AValue);
  I := 1;

  // Degrees
  S := NextNum;
  if S = '' then
    Exit;
  Val(S, ADeg, R);
  if R > 0 then
    Exit;

  // It must be the only part if negative or fractional
  Neg := ADeg < 0.0;
  Last := Neg or (Frac(ADeg) > 0.0);

  // Eat the degree symbol if present
  if not EOL and (Utf8CodePointLen(@AValue[I], 2, False) = 2)
    and (AValue[I] = #$c2) and (AValue[I + 1] = #$b0)
  then
    Inc(I, 2);

  Minutes := 0.0;
  Seconds := 0.0;

  N := 1;
  while not (EOL or Last) and (N < 3) do
  begin
    S := NextNum;
    if S = '' then
      Break
    else
    begin
      Val(S, D, R);
      // Invalid part or negative one
      if (R > 0) or (D < 0.0) then
        Exit;
      // No more parts when fractional
      Last := Frac(D) > 0.0;
      if not EOL then
        case AValue[I] of
          '''': // Munutes suffix
            begin
              Minutes := D;
              Inc(I); // Eat the '
            end;
          '"': // Seconds suffix
            begin
              Seconds := D;
              Last := True; // Last part
              Inc(I); // Eat the "
            end;
          otherwise
            if N = 1
              then Minutes := D
              else Seconds := D;
        end;
    end;
    Inc(N);
  end;

  // Merge parts
  ADeg := ADeg + Minutes / 60 + Seconds / 3600;

  // Check for N-S and E-W designators
  SkipWS;
  if not (EOL or Neg) and (AValue[I] in ['S', 's', 'W', 'w', 'N', 'n', 'E', 'e']) then
  begin
    if AValue[I] in ['S', 's', 'W', 'w']
      then ADeg := -1 * ADeg;
    Inc(I);
  end;
  SkipWS;

  // It must be entirely consumed
  Result := EOL;
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 NaN.
 Usage: CalcGeoDistance(51.53323, -2.90130, 51.29442, -2.27275, duKilometers);
}
function CalcGeoDistance(Lat1, Lon1, Lat2, Lon2: double;
  AUnits: TDistanceUnits = duKilometers): double;
begin
  // Validate
  if (Lat1 < -90.0) or (Lat1 > 90.0) then exit(NaN);
  if (Lat2 < -90.0) or (Lat2 > 90.0) then exit(NaN);

  Result := HaversineDist(Lat1, Lon1, Lat2, Lon2, EARTH_EQUATORIAL_RADIUS);
  case AUnits of
    duMeters: ;
    duKilometers: Result := Result * 1E-3;
    duMiles: Result := Result * 0.62137E-3;
  end;
end;

{ Calculate initial bearing (in °) from the start point Lat1,Lon1 to
  the end point Lat2,Lon2. No parameter checks, result normalized to 0°..360°.
}
function CalcBearing(Lat1, Lon1, Lat2, Lon2: Double): Double;
var
  latFrom, latTo, lonDiff: Double;
begin
  lonDiff := DegToRad(Lon2 - Lon1);
  latFrom := DegToRad(Lat1);
  latTo := DegToRad(Lat2);
  Result := ArcTan2(Sin(lonDiff) * Cos(latTo),
    Cos(latFrom) * Sin(latTo) - Sin(latFrom) * Cos(latTo) * Cos(lonDiff));
  Result := RadToDeg(Result);
  if Result < 0.0 then
    Result := Result + 360.0;
end;

{ Calculate end point Lat2,Lon2 by given start point Lat1,Lon1, distance
  ADist (in meters) and bearing ABearing (in °). No parameter checks,
  result Lon2 normalized to -180°..180°.
}
procedure CalcLatLon(const Lat1, Lon1, ADist, ABearing: Double; out Lat2,
  Lon2: Double);
var
  latFrom, lonFrom, brng, aD: Double;
begin
  latFrom := DegToRad(Lat1);
  lonFrom := DegToRad(Lon1);
  brng := DegToRad(ABearing);
  aD := ADist / EARTH_EQUATORIAL_RADIUS;
  Lat2 := ArcSin(Sin(latFrom) * Cos(aD) + Cos(latFrom) * Sin(aD) * Cos(brng));
  Lon2 := lonFrom + ArcTan2(Sin(brng) * Sin(aD) * Cos(latFrom),
    Cos(aD) - Sin(latFrom) * Sin(Lat2));
  Lat2 := RadToDeg(Lat2);
  Lon2 := NormalizeLon(RadToDeg(Lon2));
end;

{ Calculate midpoint Lat,Lon by given start point Lat1,Lon1 and end point
  Lat2,Lon2. No parameter checks, result Lon normalized to -180°..180°.
}
procedure CalcMidpoint(const Lat1, Lon1, Lat2, Lon2: Double; out Lat,
  Lon: Double);
var
  latFrom, lonDiff, latTo, lonTo, Bx, By: Double;
begin
  lonDiff := DegToRad(Lon2 - Lon1);
  latFrom := DegToRad(Lat1);
  latTo := DegToRad(Lat2);
  lonTo := DegToRad(Lon2);
  Bx := Cos(latTo) * Cos(lonDiff);
  By := Cos(latTo) * Sin(lonDiff);
  Lat := ArcTan2(Sin(latFrom) + Sin(latTo), Sqrt(Sqr(Cos(latFrom) + Bx) + Sqr(By)));
  Lon := lonTo + ArcTan2(By, Cos(latFrom) + By);
  Lat := RadToDeg(Lat);
  Lon := NormalizeLon(RadToDeg(Lon));
end;

{ Calculate intermediate point Lat,Lon by given start point Lat1,Lon1, end point
  Lat2,Lon2 and fraction AFrac (0.0-1.0). No parameter checks for
  Lat1,Lon1,Lat2 and Lon2. Result Lon normalized to -180°..180°.
}
procedure CalcIntermedPoint(const Lat1, Lon1, Lat2, Lon2, AFrac: Double; out
  Lat, Lon: Double);
var
  latFrom, lonFrom, latTo, lonTo: Double;
  A, B, aD, X, Y, Z: Double;
begin
  if (Lat1 = Lat2) and (Lon1 = Lon2) or (AFrac < 0.001) then
  begin
    Lat := Lat1;
    Lon := Lon1;
    Exit;
  end;
  if AFrac > 0.999 then
  begin
    Lat := Lat2;
    Lon := Lon2;
    Exit;
  end;
  aD := CalcGeoDistance(Lat1, Lon1, Lat2, Lon2) / EARTH_EQUATORIAL_RADIUS;
  latFrom := DegToRad(Lat1);
  lonFrom := DegToRad(Lon1);
  latTo := DegToRad(Lat2);
  lonTo := DegToRad(Lon2);
  A := Sin((1.0 - AFrac) * aD) / Sin(aD);
  B := Sin(AFrac * aD) / Sin(aD);
  X := A * Cos(latFrom) * Cos(lonFrom) + B * Cos(latTo) * Cos(lonTo);
  Y := A * Cos(latFrom) * Sin(lonFrom) + B * Cos(latTo) * Sin(lonTo);
  Z := A * Sin(latFrom) + B * Sin(latTo);
  Lat := ArcTan2(Z, Sqrt(Sqr(X) + Sqr(Y)));
  Lon := ArcTan2(Y, X);
  Lat := RadToDeg(Lat);
  Lon := NormalizeLon(RadToDeg(Lon));
end;

function NormalizeLon(const Lon: Double): Double;
begin
  if InRange(Lon, -180.0, 180.0)
    then Result := Lon
    else Result := FMod(Lon + 540.0, 360.0) - 180.0;
end;

{ Converts an angle given as degrees, minutes and seconds to a single
  floating point degrees value. }
function DMSToDeg(Deg, Min: Word; Sec: Double): Double;
begin
  Result := Deg + Min/60.0 + Sec/3600.0;
end;

end.