From 46db12e26f476c3b01b83fef0817c40c7d2950e8 Mon Sep 17 00:00:00 2001 From: wp_xxyyzz Date: Wed, 30 Dec 2020 14:42:32 +0000 Subject: [PATCH] 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 | 428 +++++++++++------- .../lazmapviewer/source/mvmapprovider.pas | 31 +- 2 files changed, 290 insertions(+), 169 deletions(-) diff --git a/components/lazmapviewer/source/mvengine.pas b/components/lazmapviewer/source/mvengine.pas index 069608920..0a32582e6 100644 --- a/components/lazmapviewer/source/mvengine.pas +++ b/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; - MIN_LATITUDE = -85.05112878; - MAX_LATITUDE = 85.05112878; - MIN_LONGITUDE = -180; - MAX_LONGITUDE = 180; - SHIFT = 2 * pi * EARTH_RADIUS / 2.0; + EARTH_EQUATORIAL_RADIUS = 6378137; + EARTH_POLAR_RADIUS = 6356752.3142; + EARTH_CIRCUMFERENCE = 2 * pi * EARTH_EQUATORIAL_RADIUS; + EARTH_ECCENTRICITY = sqrt(1-( exp(2*ln(EARTH_POLAR_RADIUS)) / exp(2*ln(EARTH_EQUATORIAL_RADIUS)))); + // power() doesn't work in interface - -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 MapWinToLonLat(const aWin: TMapWindow; aPt : TPoint) : TRealPoint; - Procedure CalculateWin(var aWin: TMapWindow); - Procedure Redraw(const aWin: TmapWindow); + function DegreesToMapPixels(const AWin: TMapWindow; ALonLat: 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; 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 LonLatToWorldScreen(aPt: TRealPoint): TPoint; + function LonLatToScreen(ALonLat: 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; - MinZoom, MaxZoom, NbSvr: integer; GetSvrStr: TGetSvrStr; GetXStr: TGetValStr; - GetYStr: TGetValStr; GetZStr: TGetValStr): TMapProvider; +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 @@ -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; - res: Extended; - px, py: Int64; + PixelLocation: TPoint; // review: coth: Should it use Int64? begin - mx := aWin.Center.Lon * SHIFT / 180.0; - my := ln( tan((90 - aWin.Center.Lat) * pi / 360.0 )) / (pi / 180.0); - my := my * SHIFT / 180.0; + case AWin.MapProvider.ProjectionType of + ptEPSG3857: PixelLocation := DegreesToPixelsEPSG3857(AWin, AWin.Center); + ptEPSG3395: PixelLocation := DegreesToPixelsEPSG3395(AWin, AWin.Center); + else PixelLocation := DegreesToPixelsEPSG3857(AWin, AWin.Center); + end; - res := (2 * pi * EARTH_RADIUS) / (TILE_SIZE * (1 shl aWin.Zoom)); - px := Round((mx + shift) / res); - py := Round((my + shift) / res); - - aWin.X := aWin.Width div 2 - px; - aWin.Y := aWin.Height div 2 - py; + AWin.X := Int64(AWin.Width div 2) - PixelLocation.x; + AWin.Y := Int64(AWin.Height div 2) - PixelLocation.y; end; procedure TMapViewerEngine.CancelCurrentDrawing; @@ -572,88 +570,193 @@ begin (aTile.Y >= 0) and (aTile.Y <= tiles-1); end; -function TMapViewerEngine.LonLatToMapWin(const aWin: TMapWindow; - aPt: TRealPoint): TPoint; +function TMapViewerEngine.DegreesToMapPixels(const AWin: TMapWindow; + ALonLat: TRealPoint): TPoint; var - tiles: Int64; - circumference: Int64; - res: Extended; - tmpX,tmpY : Double; + pixelLocation: TPoint; begin - tiles := 1 shl aWin.Zoom; - circumference := tiles * TILE_SIZE; - tmpX := ((aPt.Lon+ 180.0)*circumference)/360.0; - - res := (2 * pi * EARTH_RADIUS) / circumference; - - tmpY := -aPt.Lat; - 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); + case AWin.MapProvider.ProjectionType of + ptEPSG3395: pixelLocation := DegreesToPixelsEPSG3395(AWin, ALonLat); + ptEPSG3857: pixelLocation := DegreesToPixelsEPSG3857(AWin, ALonLat); + else pixelLocation := DegreesToPixelsEPSG3857(AWin, ALonLat); + end; + Result.X := pixelLocation.x + AWin.X; + Result.Y := pixelLocation.y + AWin.Y; 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; - aPt: TPoint): TRealPoint; +function TMapViewerEngine.MapPixelsToDegrees(const AWin: TMapWindow; + APoint: TPoint): TRealPoint; var - tiles: Int64; - circumference: Int64; - lat: Extended; - res: Extended; + iMapWidth: Int64; mPoint : TPoint; begin - tiles := 1 shl aWin.Zoom; - circumference := tiles * TILE_SIZE; + // review: coth: respective projection check. move to subfunctions? figure out what did i mean here... + iMapWidth := Round(IntPower(2, AWin.Zoom)) * TILE_SIZE; - mPoint.X := aPt.X - aWin.X; - mPoint.Y := aPt.Y - aWin.Y; + mPoint.X := EnsureRange(APoint.X - AWin.X, 0, iMapWidth); + mPoint.Y := EnsureRange(APoint.Y - AWin.Y, 0, iMapWidth); - if mPoint.X < 0 then - mPoint.X := 0 - else - if mPoint.X > circumference then - mPoint.X := circumference; + case aWin.MapProvider.ProjectionType of + ptEPSG3857: Result := PixelsToDegreesEPSG3857(mPoint, AWin.Zoom); + ptEPSG3395: Result := PixelsToDegreesEPSG3395(mPoint, AWin.Zoom); + else Result := PixelsToDegreesEPSG3857(mPoint, AWin.Zoom); + end; +end; - if mPoint.Y < 0 then - mPoint.Y := 0 - else - if mPoint.Y > circumference then - mPoint.Y := circumference; +function TMapViewerEngine.PixelsToDegreesEPSG3857(APoint: TPoint; Zoom: Integer): TRealPoint; +const + MIN_LATITUDE = -85.05112878; + MAX_LATITUDE = 85.05112878; + 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; - lat := ((mPoint.Y * res - SHIFT) / SHIFT) * 180.0; + Result.Lat := Math.EnsureRange(Result.Lat, MIN_LATITUDE, MAX_LATITUDE); + Result.Lon := Math.EnsureRange(Result.Lon, MIN_LONGITUDE, MAX_LONGITUDE); +end; - lat := radtodeg (2 * arctan( exp( lat * pi / 180.0)) - pi / 2.0); - Result.Lat := -lat; +Function TMapViewerEngine.PixelsToDegreesEPSG3395(APoint: TPoint; Zoom: Integer): TRealPoint; - if Result.Lat > MAX_LATITUDE then - Result.Lat := MAX_LATITUDE - else - if Result.Lat < MIN_LATITUDE then - Result.Lat := MIN_LATITUDE; + 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; - if Result.Lon > MAX_LONGITUDE then - Result.Lon := MAX_LONGITUDE - else - if Result.Lon < MIN_LONGITUDE then - Result.Lon := MIN_LONGITUDE; +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; + 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', - 'http://mt%serv%.google.com/vt/lyrs=m@145&v=w2.104&x=%x%&y=%y%&z=%z%', - 0, 19, 4, nil); - AddMapProvider('Google Hybrid', - 'http://mt%serv%.google.com/vt/lyrs=h@145&v=w2.104&x=%x%&y=%y%&z=%z%', - 0, 19, 4, nil); - AddMapProvider('Google Physical', - 'http://mt%serv%.google.com/vt/lyrs=t@145&v=w2.104&x=%x%&y=%y%&z=%z%', - 0, 19, 4, nil); + // dev links + //https://gis-lab.info/forum/viewtopic.php?f=19&t=19763 + //https://a.tile.openstreetmap.org/16/51693/32520.png + //https://vec01.maps.yandex.net/tiles?l=map&x=51693+570&y=32520&z=16&scale=1&lang=ru_RU + //https://www.linux.org.ru/forum/development/9038716 + //https://wiki.openstreetmap.org/wiki/Tiles + //https://pubs.usgs.gov/pp/1395/report.pdf + //https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Tile_numbers_to_lon..2Flat. + //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, - 1, 19, 4, @GetSvrBase1); - AddMapProvider('Here Maps Grey', HERE1 + 'normal.day.grey' + HERE2, - 1, 19, 4, @GetSvrBase1); - AddMapProvider('Here Maps Reduced', HERE1 + 'reduced.day' + HERE2, - 1, 19, 4, @GetSvrBase1); - AddMapProvider('Here Maps Transit', HERE1 + 'normal.day.transit' + HERE2, - 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); + AddMapProvider('Here WeGo Map', ptEPSG3857, HERE1 + 'normal.day' + HERE2, 1, 19, 4, @GetSvrBase1); + AddMapProvider('Here WeGo Grey Map', ptEPSG3857, HERE1 + 'normal.day.grey' + HERE2, 1, 19, 4, @GetSvrBase1); + AddMapProvider('Here WeGo Reduced Map', ptEPSG3857, HERE1 + 'reduced.day' + HERE2, 1, 19, 4, @GetSvrBase1); + AddMapProvider('Here WeGo Transit Map', ptEPSG3857, HERE1 + 'normal.day.transit' + HERE2, 1, 19, 4, @GetSvrBase1); + AddMapProvider('Here WeGo POI Map', ptEPSG3857, HERE1 + 'normal.day' + HERE2 + '&pois', 1, 19, 4, @GetSvrBase1); + AddMapProvider('Here WeGo Pedestrian Map', ptEPSG3857, HERE1 + 'pedestrian.day' + HERE2, 1, 19, 4, @GetSvrBase1); + AddMapProvider('Here WeGo DreamWorks Map', ptEPSG3857, 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', - 'https://tile.openweathermap.org/map/clouds_new/%z%/%x%/%y%.png?appid=' + OpenWeatherMap_ApiKey, - 1, 19, 1, nil); - AddMapProvider('OpenWeatherMap Precipitation', - 'https://tile.openweathermap.org/map/precipitation_new/%z%/%x%/%y%.png?appid=' + OpenWeatherMap_ApiKey, - 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); + AddMapProvider('OpenWeatherMap Clouds', ptEPSG3857, 'https://tile.openweathermap.org/map/clouds_new/%z%/%x%/%y%.png?appid=' + OpenWeatherMap_ApiKey, 1, 19, 1, nil); + AddMapProvider('OpenWeatherMap Precipitation', ptEPSG3857, 'https://tile.openweathermap.org/map/precipitation_new/%z%/%x%/%y%.png?appid=' + OpenWeatherMap_ApiKey, 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 Temperature', ptEPSG3857, 'https://tile.openweathermap.org/map/temp_new/%z%/%x%/%y%.png?appid=' + OpenWeatherMap_ApiKey, 1, 19, 1, nil); + AddMapProvider('OpenWeatherMap Wind', ptEPSG3857, '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.BottomRight := MapWinToLonLat(tmpWin, BottomRight); + visArea.TopLeft := MapPixelsToDegrees(tmpWin, TopLeft); + 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: ; diff --git a/components/lazmapviewer/source/mvmapprovider.pas b/components/lazmapviewer/source/mvmapprovider.pas index d2abb3afa..bee6b6e66 100644 --- a/components/lazmapviewer/source/mvmapprovider.pas +++ b/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; - GetSvrStr: TGetSvrStr; GetXStr: TGetValStr; GetYStr: TGetValStr; - GetZStr: TGetValStr); +procedure TMapProvider.AddURL(Url: String; ProjectionType: TProjectionType; + NbSvr, aMinZoom, aMaxZoom: integer; GetSvrStr: TGetSvrStr; + 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