From fec5780da2c3305180bc58a4d7ac306b8067a995 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Martin=20T=C5=AFma?= Date: Tue, 8 Jan 2019 21:42:28 +0100 Subject: [PATCH] Added support for polar stereographic projection closes #181 --- gpxsee.pro | 7 +- pkg/csv/pcs.csv | 1 + src/map/coordinatesystem.cpp | 1 + src/map/geotiff.cpp | 103 +++++++++--------- src/map/pcs.cpp | 1 + src/map/polarstereographic.cpp | 192 +++++++++++++++++++++++++++++++++ src/map/polarstereographic.h | 34 ++++++ src/map/projection.cpp | 7 ++ src/map/rectd.cpp | 53 +++++++++ src/map/rectd.h | 10 ++ src/map/wmsmap.cpp | 16 ++- src/map/wmsmap.h | 3 +- src/map/wmtsmap.cpp | 20 ++-- src/map/wmtsmap.h | 3 +- 14 files changed, 377 insertions(+), 74 deletions(-) create mode 100644 src/map/polarstereographic.cpp create mode 100644 src/map/polarstereographic.h create mode 100644 src/map/rectd.cpp diff --git a/gpxsee.pro b/gpxsee.pro index 7f7717ff..fde63907 100644 --- a/gpxsee.pro +++ b/gpxsee.pro @@ -148,7 +148,8 @@ HEADERS += src/common/config.h \ src/data/oziparsers.h \ src/data/locparser.h \ src/data/slfparser.h \ - src/data/dem.h + src/data/dem.h \ + src/map/polarstereographic.h SOURCES += src/main.cpp \ src/common/coordinates.cpp \ src/common/rectc.cpp \ @@ -257,7 +258,9 @@ SOURCES += src/main.cpp \ src/data/oziparsers.cpp \ src/data/locparser.cpp \ src/data/slfparser.cpp \ - src/data/dem.cpp + src/data/dem.cpp \ + src/map/polarstereographic.cpp \ + src/map/rectd.cpp RESOURCES += gpxsee.qrc TRANSLATIONS = lang/gpxsee_en.ts \ lang/gpxsee_cs.ts \ diff --git a/pkg/csv/pcs.csv b/pkg/csv/pcs.csv index 20479373..35d1bbe7 100644 --- a/pkg/csv/pcs.csv +++ b/pkg/csv/pcs.csv @@ -859,6 +859,7 @@ RT90 2.5 gon V,3021,4124,19929,9001,9807,4530,8801,0,9110,8802,15.48298,9110,880 RT90 0 gon,3022,4124,17336,9001,9807,4530,8801,0,9110,8802,18.03298,9110,8805,1,9201,8806,1500000,9001,8807,0,9001,,,,,, RT90 2.5 gon O,3023,4124,17337,9001,9807,4530,8801,0,9110,8802,20.18298,9110,8805,1,9201,8806,1500000,9001,8807,0,9001,,,,,, RT90 5 gon O,3024,4124,17338,9001,9807,4530,8801,0,9110,8802,22.33298,9110,8805,1,9201,8806,1500000,9001,8807,0,9001,,,,,, +WGS 84 / Antarctic Polar Stereographic,3031,4326,19992,9001,9829,4490,8806,0,9001,8807,0,9001,8832,-71,9102,8833,0,9102,,,,,,,,, WGS 84 / Australian Antarctic Lambert,3033,4326,19994,9001,9802,4400,8821,-50,9110,8822,70,9110,8823,-68.3,9110,8824,-74.3,9110,8826,6000000,9001,8827,6000000,9001,,, ETRS89 / LCC Europe,3034,4258,19985,9001,9802,4500,8821,52,9102,8822,10,9102,8823,35,9102,8824,65,9102,8826,4000000,9001,8827,2800000,9001,,, ETRS89 / LAEA Europe,3035,4258,19986,9001,9820,4532,8801,52,9102,8802,10,9102,8806,4321000,9001,8807,3210000,9001,,,,,,,,, diff --git a/src/map/coordinatesystem.cpp b/src/map/coordinatesystem.cpp index d2832644..0a7a3098 100644 --- a/src/map/coordinatesystem.cpp +++ b/src/map/coordinatesystem.cpp @@ -18,6 +18,7 @@ CoordinateSystem::CoordinateSystem(int code) case 4467: case 4469: case 4470: + case 4490: case 4495: case 4496: case 4497: diff --git a/src/map/geotiff.cpp b/src/map/geotiff.cpp index f29d85ae..ed4300f0 100644 --- a/src/map/geotiff.cpp +++ b/src/map/geotiff.cpp @@ -3,58 +3,60 @@ #include "geotiff.h" -#define TIFF_DOUBLE 12 -#define TIFF_SHORT 3 -#define TIFF_LONG 4 +#define TIFF_DOUBLE 12 +#define TIFF_SHORT 3 +#define TIFF_LONG 4 -#define ModelPixelScaleTag 33550 -#define ModelTiepointTag 33922 -#define ModelTransformationTag 34264 -#define GeoKeyDirectoryTag 34735 -#define GeoDoubleParamsTag 34736 +#define ModelPixelScaleTag 33550 +#define ModelTiepointTag 33922 +#define ModelTransformationTag 34264 +#define GeoKeyDirectoryTag 34735 +#define GeoDoubleParamsTag 34736 -#define GTModelTypeGeoKey 1024 -#define GTRasterTypeGeoKey 1025 -#define GeographicTypeGeoKey 2048 -#define GeogGeodeticDatumGeoKey 2050 -#define GeogPrimeMeridianGeoKey 2051 -#define GeogAngularUnitsGeoKey 2054 -#define GeogEllipsoidGeoKey 2056 -#define GeogAzimuthUnitsGeoKey 2060 -#define ProjectedCSTypeGeoKey 3072 -#define ProjectionGeoKey 3074 -#define ProjCoordTransGeoKey 3075 -#define ProjLinearUnitsGeoKey 3076 -#define ProjStdParallel1GeoKey 3078 -#define ProjStdParallel2GeoKey 3079 -#define ProjNatOriginLongGeoKey 3080 -#define ProjNatOriginLatGeoKey 3081 -#define ProjFalseEastingGeoKey 3082 -#define ProjFalseNorthingGeoKey 3083 -#define ProjFalseOriginLongGeoKey 3084 -#define ProjFalseOriginLatGeoKey 3085 -#define ProjFalseOriginEastingGeoKey 3086 -#define ProjFalseOriginNorthingGeoKey 3087 -#define ProjCenterLongGeoKey 3088 -#define ProjCenterLatGeoKey 3089 -#define ProjCenterEastingGeoKey 3090 -#define ProjCenterNorthingGeoKey 3091 -#define ProjScaleAtNatOriginGeoKey 3092 -#define ProjScaleAtCenterGeoKey 3093 -#define ProjAzimuthAngleGeoKey 3094 -#define ProjRectifiedGridAngleGeoKey 3096 +#define GTModelTypeGeoKey 1024 +#define GTRasterTypeGeoKey 1025 +#define GeographicTypeGeoKey 2048 +#define GeogGeodeticDatumGeoKey 2050 +#define GeogPrimeMeridianGeoKey 2051 +#define GeogAngularUnitsGeoKey 2054 +#define GeogEllipsoidGeoKey 2056 +#define GeogAzimuthUnitsGeoKey 2060 +#define ProjectedCSTypeGeoKey 3072 +#define ProjectionGeoKey 3074 +#define ProjCoordTransGeoKey 3075 +#define ProjLinearUnitsGeoKey 3076 +#define ProjStdParallel1GeoKey 3078 +#define ProjStdParallel2GeoKey 3079 +#define ProjNatOriginLongGeoKey 3080 +#define ProjNatOriginLatGeoKey 3081 +#define ProjFalseEastingGeoKey 3082 +#define ProjFalseNorthingGeoKey 3083 +#define ProjFalseOriginLongGeoKey 3084 +#define ProjFalseOriginLatGeoKey 3085 +#define ProjFalseOriginEastingGeoKey 3086 +#define ProjFalseOriginNorthingGeoKey 3087 +#define ProjCenterLongGeoKey 3088 +#define ProjCenterLatGeoKey 3089 +#define ProjCenterEastingGeoKey 3090 +#define ProjCenterNorthingGeoKey 3091 +#define ProjScaleAtNatOriginGeoKey 3092 +#define ProjScaleAtCenterGeoKey 3093 +#define ProjAzimuthAngleGeoKey 3094 +#define ProjStraightVertPoleLongGeoKey 3095 +#define ProjRectifiedGridAngleGeoKey 3096 -#define ModelTypeProjected 1 -#define ModelTypeGeographic 2 -#define ModelTypeGeocentric 3 +#define ModelTypeProjected 1 +#define ModelTypeGeographic 2 +#define ModelTypeGeocentric 3 -#define CT_TransverseMercator 1 -#define CT_ObliqueMercator 3 -#define CT_Mercator 7 -#define CT_LambertConfConic_2SP 8 -#define CT_LambertConfConic_1SP 9 -#define CT_LambertAzimEqualArea 10 -#define CT_AlbersEqualArea 11 +#define CT_TransverseMercator 1 +#define CT_ObliqueMercator 3 +#define CT_Mercator 7 +#define CT_LambertConfConic_2SP 8 +#define CT_LambertConfConic_1SP 9 +#define CT_LambertAzimEqualArea 10 +#define CT_AlbersEqualArea 11 +#define CT_PolarStereographic 15 #define IS_SET(map, key) \ @@ -261,6 +263,7 @@ bool GeoTIFF::readKeys(TIFFFile &file, Ctx &ctx, QMap &kv) const case ProjCenterNorthingGeoKey: case ProjFalseOriginEastingGeoKey: case ProjFalseOriginNorthingGeoKey: + case ProjStraightVertPoleLongGeoKey: if (!readGeoValue(file, ctx.values, entry.ValueOffset, value.DOUBLE)) return false; @@ -342,6 +345,8 @@ Projection::Method GeoTIFF::method(QMap &kv) return Projection::Method(9820); case CT_AlbersEqualArea: return Projection::Method(9822); + case CT_PolarStereographic: + return Projection::Method(9829); default: _errorString = QString("%1: unknown coordinate transformation method") .arg(kv.value(ProjCoordTransGeoKey).SHORT); @@ -407,6 +412,8 @@ bool GeoTIFF::projectedModel(QMap &kv) lon0 = au.toDegrees(kv.value(ProjCenterLongGeoKey).DOUBLE); else if (kv.contains(ProjFalseOriginLongGeoKey)) lon0 = au.toDegrees(kv.value(ProjFalseOriginLongGeoKey).DOUBLE); + else if (kv.contains(ProjStraightVertPoleLongGeoKey)) + lon0 = au.toDegrees(kv.value(ProjStraightVertPoleLongGeoKey).DOUBLE); else lon0 = NAN; if (kv.contains(ProjScaleAtNatOriginGeoKey)) diff --git a/src/map/pcs.cpp b/src/map/pcs.cpp index a43b50dd..5f4711e7 100644 --- a/src/map/pcs.cpp +++ b/src/map/pcs.cpp @@ -32,6 +32,7 @@ static bool parameter(int key, double val, int units, Projection::Setup &setup) case 8801: case 8811: case 8821: + case 8832: {AngularUnits au(units); if (au.isNull()) return false; diff --git a/src/map/polarstereographic.cpp b/src/map/polarstereographic.cpp new file mode 100644 index 00000000..6d790d84 --- /dev/null +++ b/src/map/polarstereographic.cpp @@ -0,0 +1,192 @@ +/* + * Based on libgeotrans with the following Source Code Disclaimer: + +1. The GEOTRANS source code ("the software") is provided free of charge by +the National Imagery and Mapping Agency (NIMA) of the United States +Department of Defense. Although NIMA makes no copyright claim under Title 17 +U.S.C., NIMA claims copyrights in the source code under other legal regimes. +NIMA hereby grants to each user of the software a license to use and +distribute the software, and develop derivative works. + +2. Warranty Disclaimer: The software was developed to meet only the internal +requirements of the U.S. National Imagery and Mapping Agency. The software +is provided "as is," and no warranty, express or implied, including but not +limited to the implied warranties of merchantability and fitness for +particular purpose or arising by statute or otherwise in law or from a +course of dealing or usage in trade, is made by NIMA as to the accuracy and +functioning of the software. + +3. NIMA and its personnel are not required to provide technical support or +general assistance with respect to the software. + +4. Neither NIMA nor its personnel will be liable for any claims, losses, or +damages arising from or connected with the use of the software. The user +agrees to hold harmless the United States National Imagery and Mapping +Agency. The user's sole and exclusive remedy is to stop using the software. + +5. NIMA requests that products developed using the software credit the +source of the software with the following statement, "The product was +developed using GEOTRANS, a product of the National Imagery and Mapping +Agency and U.S. Army Engineering Research and Development Center." + +6. For any products developed using the software, NIMA requires a disclaimer +that use of the software does not indicate endorsement or approval of the +product by the Secretary of Defense or the National Imagery and Mapping +Agency. Pursuant to the United States Code, 10 U.S.C. Sec. 2797, the name of +the National Imagery and Mapping Agency, the initials "NIMA", the seal of +the National Imagery and Mapping Agency, or any colorable imitation thereof +shall not be used to imply approval, endorsement, or authorization of a +product without prior written permission from United States Secretary of +Defense. + +*/ + +#include "ellipsoid.h" +#include "polarstereographic.h" + + +#define POLAR_POW(EsSin) pow((1.0 - EsSin) / (1.0 + EsSin), _es_OVER_2) + +PolarStereographic::PolarStereographic(const Ellipsoid *ellipsoid, + double latitudeOrigin, double longitudeOrigin, + double falseEasting, double falseNorthing) +{ + _two_a = 2.0 * ellipsoid->radius(); + + if (longitudeOrigin > M_PI) + longitudeOrigin -= 2*M_PI; + + if (latitudeOrigin < 0) { + _southernHemisphere = 1; + _originLatitude = -deg2rad(latitudeOrigin); + _originLongitude = -deg2rad(longitudeOrigin); + } else { + _southernHemisphere = 0; + _originLatitude = deg2rad(latitudeOrigin); + _originLongitude = deg2rad(longitudeOrigin); + } + _falseEasting = falseEasting; + _falseNorthing = falseNorthing; + + double es2 = ellipsoid->es(); + _es = sqrt(es2); + _es_OVER_2 = _es / 2.0; + + if (fabs(fabs(_originLatitude) - M_PI_2) > 1.0e-10) { + double slat = sin(_originLatitude); + double essin = _es * slat; + double pow_es = POLAR_POW(essin); + double clat = cos(_originLatitude); + _mc = clat / sqrt(1.0 - essin * essin); + _a_mc = ellipsoid->radius() * _mc; + _tc = tan(M_PI_4 - _originLatitude / 2.0) / pow_es; + } else { + double one_PLUS_es = 1.0 + _es; + double one_MINUS_es = 1.0 - _es; + _e4 = sqrt(pow(one_PLUS_es, one_PLUS_es) * pow(one_MINUS_es, + one_MINUS_es)); + } +} + +PointD PolarStereographic::ll2xy(const Coordinates &c) const +{ + double Easting, Northing; + double Longitude = deg2rad(c.lon()); + double Latitude = deg2rad(c.lat()); + + + if (fabs(fabs(Latitude) - M_PI_2) < 1.0e-10) { + Easting = 0.0; + Northing = 0.0; + } else { + if (_southernHemisphere != 0) { + Longitude *= -1.0; + Latitude *= -1.0; + } + + double dlam = Longitude - _originLongitude; + if (dlam > M_PI) + dlam -= 2*M_PI; + if (dlam < -M_PI) + dlam += 2*M_PI; + + double slat = sin(Latitude); + double essin = _es * slat; + double pow_es = POLAR_POW(essin); + double t = tan(M_PI_4 - Latitude / 2.0) / pow_es; + + double rho; + if (fabs(fabs(_originLatitude) - M_PI_2) > 1.0e-10) + rho = _a_mc * t / _tc; + else + rho = _two_a * t / _e4; + + Easting = rho * sin(dlam) + _falseEasting; + + if (_southernHemisphere != 0) { + Easting *= -1.0; + Northing = rho * cos(dlam) + _falseNorthing; + } else + Northing = -rho * cos(dlam) + _falseNorthing; + } + + return PointD(Easting, Northing); +} + +Coordinates PolarStereographic::xy2ll(const PointD &p) const +{ + double Latitude, Longitude; + double dy = p.y() - _falseNorthing; + double dx = p.x() - _falseEasting; + + if ((dy == 0.0) && (dx == 0.0)) { + Latitude = M_PI_2; + Longitude = _originLongitude; + } else { + if (_southernHemisphere != 0) { + dy *= -1.0; + dx *= -1.0; + } + + double rho = sqrt(dx * dx + dy * dy); + double t; + if (fabs(fabs(_originLatitude) - M_PI_2) > 1.0e-10) + t = rho * _tc / (_a_mc); + else + t = rho * _e4 / (_two_a); + double PHI = M_PI_2 - 2.0 * atan(t); + + double tempPHI = 0.0; + while (fabs(PHI - tempPHI) > 1.0e-10) { + tempPHI = PHI; + double sin_PHI = sin(PHI); + double essin = _es * sin_PHI; + double pow_es = POLAR_POW(essin); + PHI = M_PI_2 - 2.0 * atan(t * pow_es); + } + Latitude = PHI; + Longitude = _originLongitude + atan2(dx, -dy); + + if (Longitude > M_PI) + Longitude -= 2*M_PI; + else if (Longitude < -M_PI) + Longitude += 2*M_PI; + + if (Latitude > M_PI_2) + Latitude = M_PI_2; + else if (Latitude < -M_PI_2) + Latitude = -M_PI_2; + + if (Longitude > M_PI) + Longitude = M_PI; + else if (Longitude < -M_PI) + Longitude = -M_PI; + } + + if (_southernHemisphere != 0) { + Latitude *= -1.0; + Longitude *= -1.0; + } + + return Coordinates(rad2deg(Longitude), rad2deg(Latitude)); +} diff --git a/src/map/polarstereographic.h b/src/map/polarstereographic.h new file mode 100644 index 00000000..7ac606ff --- /dev/null +++ b/src/map/polarstereographic.h @@ -0,0 +1,34 @@ +#ifndef POLARSTEREOGRAPHIC_H +#define POLARSTEREOGRAPHIC_H + +#include "ct.h" + +class Ellipsoid; + +class PolarStereographic : public CT +{ +public: + PolarStereographic(const Ellipsoid *ellipsoid, double latitudeOrigin, + double longitudeOrigin, double falseEasting, double falseNorthing); + virtual CT *clone() const {return new PolarStereographic(*this);} + + virtual PointD ll2xy(const Coordinates &c) const; + virtual Coordinates xy2ll(const PointD &p) const; + +private: + double _originLatitude; + double _originLongitude; + double _falseEasting; + double _falseNorthing; + + double _a_mc; + double _es; + double _es_OVER_2; + double _two_a; + double _mc; + double _tc; + double _e4; + int _southernHemisphere; +}; + +#endif // POLARSTEREOGRAPHIC_H diff --git a/src/map/projection.cpp b/src/map/projection.cpp index b948c135..25f28d40 100644 --- a/src/map/projection.cpp +++ b/src/map/projection.cpp @@ -6,6 +6,7 @@ #include "albersequal.h" #include "lambertazimuthal.h" #include "krovak.h" +#include "polarstereographic.h" #include "latlon.h" #include "gcs.h" #include "pcs.h" @@ -25,6 +26,7 @@ Projection::Method::Method(int id) case 9819: case 9820: case 9822: + case 9829: _id = id; break; default: @@ -87,6 +89,11 @@ Projection::Projection(const PCS *pcs) : _gcs(pcs->gcs()), _units(pcs->units()), setup.longitudeOrigin(), setup.falseEasting(), setup.falseNorthing()); break; + case 9829: + _ct = new PolarStereographic(ellipsoid, setup.latitudeOrigin(), + setup.longitudeOrigin(), setup.falseEasting(), + setup.falseNorthing()); + break; default: _ct = 0; } diff --git a/src/map/rectd.cpp b/src/map/rectd.cpp new file mode 100644 index 00000000..13ddf792 --- /dev/null +++ b/src/map/rectd.cpp @@ -0,0 +1,53 @@ +#include "common/rectc.h" +#include "projection.h" +#include "rectd.h" + + +#define SAMPLE_POINTS 100 + +static void growRect(const Projection &proj, const Coordinates &c, RectD &rect) +{ + if (c.isNull()) + return; + + PointD p(proj.ll2xy(c)); + + if (rect.isNull()) + rect = RectD(p, p); + else { + if (p.x() < rect.left()) + rect.setLeft(p.x()); + if (p.x() > rect.right()) + rect.setRight(p.x()); + if (p.y() < rect.bottom()) + rect.setBottom(p.y()); + if (p.y() > rect.top()) + rect.setTop(p.y()); + } +} + +RectD::RectD(const RectC &rect, const Projection &proj) +{ + RectD prect; + double dx = (rect.right() - rect.left()) / SAMPLE_POINTS; + double dy = (rect.top() - rect.bottom()) / SAMPLE_POINTS; + + growRect(proj, rect.topLeft(), prect); + + if (dx > 0) { + for (int i = 0; i <= SAMPLE_POINTS; i++) { + double x = rect.left() + i * dx; + growRect(proj, Coordinates(x, rect.bottom()), prect); + growRect(proj, Coordinates(x, rect.top()), prect); + } + } + if (dy > 0) { + for (int i = 0; i <= SAMPLE_POINTS; i++ ) { + double y = rect.bottom() + i * dy; + growRect(proj, Coordinates(rect.left(), y), prect); + growRect(proj, Coordinates(rect.right(), y), prect); + } + } + + *this = prect; +} diff --git a/src/map/rectd.h b/src/map/rectd.h index 3b23547c..497f445c 100644 --- a/src/map/rectd.h +++ b/src/map/rectd.h @@ -3,12 +3,16 @@ #include "pointd.h" +class RectC; +class Projection; + class RectD { public: RectD() {} RectD(const PointD &topLeft, const PointD &bottomRight) : _tl(topLeft), _br(bottomRight) {} + RectD(const RectC &rect, const Projection &proj); PointD topLeft() const {return _tl;} PointD bottomRight() const {return _br;} @@ -18,6 +22,11 @@ public: double top() const {return _tl.y();} double bottom() const {return _br.y();} + void setLeft(double val) {_tl.rx() = val;} + void setRight(double val) {_br.rx() = val;} + void setTop(double val) {_tl.ry() = val;} + void setBottom(double val) {_br.ry() = val;} + double width() const {return (right() - left());} double height() const {return (top() - bottom());} @@ -26,6 +35,7 @@ public: && p.y() >= bottom());} bool isNull() const {return _tl.isNull() && _br.isNull();} + bool isValid() const {return !(_tl.isNull() || _br.isNull());} private: PointD _tl, _br; diff --git a/src/map/wmsmap.cpp b/src/map/wmsmap.cpp index 2a5ebc52..6d9dea85 100644 --- a/src/map/wmsmap.cpp +++ b/src/map/wmsmap.cpp @@ -65,12 +65,8 @@ void WMSMap::updateTransform() double pixelSpan = sd2res(_zooms.at(_zoom)); if (_projection.isGeographic()) pixelSpan /= deg2rad(WGS84_RADIUS); - double sx = _bbox.width() / pixelSpan; - double sy = _bbox.height() / pixelSpan; - - ReferencePoint tl(PointD(0, 0), _bbox.topLeft()); - ReferencePoint br(PointD(sx, sy), _bbox.bottomRight()); - _transform = Transform(tl, br); + _transform = Transform(ReferencePoint(PointD(0, 0), + _projection.ll2xy(_bbox.topLeft())), PointD(pixelSpan, pixelSpan)); } bool WMSMap::loadWMS() @@ -84,8 +80,8 @@ bool WMSMap::loadWMS() } _projection = wms.projection(); - _bbox = RectD(_projection.ll2xy(wms.boundingBox().topLeft()), - _projection.ll2xy(wms.boundingBox().bottomRight())); + _bbox = wms.boundingBox(); + _bounds = RectD(_bbox, _projection); _tileLoader->setUrl(tileUrl(wms.version())); if (wms.version() >= "1.3.0") { @@ -124,8 +120,8 @@ void WMSMap::clearCache() QRectF WMSMap::bounds() { - return QRectF(_transform.proj2img(_bbox.topLeft()) / _mapRatio, - _transform.proj2img(_bbox.bottomRight()) / _mapRatio); + return QRectF(_transform.proj2img(_bounds.topLeft()) / _mapRatio, + _transform.proj2img(_bounds.bottomRight()) / _mapRatio); } int WMSMap::zoomFit(const QSize &size, const RectC &rect) diff --git a/src/map/wmsmap.h b/src/map/wmsmap.h index 974a606f..7306ff95 100644 --- a/src/map/wmsmap.h +++ b/src/map/wmsmap.h @@ -55,7 +55,8 @@ private: Transform _transform; CoordinateSystem _cs; QVector _zooms; - RectD _bbox; + RectC _bbox; + RectD _bounds; int _zoom; qreal _mapRatio; diff --git a/src/map/wmtsmap.cpp b/src/map/wmtsmap.cpp index 98ee9df4..05ff6392 100644 --- a/src/map/wmtsmap.cpp +++ b/src/map/wmtsmap.cpp @@ -22,10 +22,10 @@ bool WMTSMap::loadWMTS() return false; } - _bounds = wmts.bounds(); _zooms = wmts.zooms(); _projection = wmts.projection(); _tileLoader->setUrl(wmts.tileUrl()); + _bounds = RectD(wmts.bounds(), _projection); if (_setup.coordinateSystem().axisOrder() == CoordinateSystem::Unknown) _cs = _projection.coordinateSystem(); @@ -77,14 +77,8 @@ void WMTSMap::updateTransform() double pixelSpan = sd2res(z.scaleDenominator()); if (_projection.isGeographic()) pixelSpan /= deg2rad(WGS84_RADIUS); - PointD tileSpan(z.tile().width() * pixelSpan, z.tile().height() * pixelSpan); - PointD bottomRight(topLeft.x() + tileSpan.x() * z.matrix().width(), - topLeft.y() - tileSpan.y() * z.matrix().height()); - - ReferencePoint tl(PointD(0, 0), topLeft); - ReferencePoint br(PointD(z.tile().width() * z.matrix().width(), - z.tile().height() * z.matrix().height()), bottomRight); - _transform = Transform(tl, br); + _transform = Transform(ReferencePoint(PointD(0, 0), topLeft), + PointD(pixelSpan, pixelSpan)); } QRectF WMTSMap::bounds() @@ -99,9 +93,11 @@ QRectF WMTSMap::bounds() * z.tile().height()), QSize(z.tile().width() * z.limits().width(), z.tile().height() * z.limits().height())); - bounds = _bounds.isValid() ? QRectF(ll2xy(_bounds.topLeft()), - ll2xy(_bounds.bottomRight())) : QRectF(); - return _bounds.isValid() ? tileBounds.intersected(bounds) : tileBounds; + if (_bounds.isValid()) + bounds = QRectF(_transform.proj2img(_bounds.topLeft()) + / coordinatesRatio(), _transform.proj2img(_bounds.bottomRight()) + / coordinatesRatio()); + return bounds.isValid() ? tileBounds.intersected(bounds) : tileBounds; } int WMTSMap::zoomFit(const QSize &size, const RectC &rect) diff --git a/src/map/wmtsmap.h b/src/map/wmtsmap.h index e6ed6e87..e86728fc 100644 --- a/src/map/wmtsmap.h +++ b/src/map/wmtsmap.h @@ -4,6 +4,7 @@ #include "transform.h" #include "projection.h" #include "map.h" +#include "rectd.h" #include "wmts.h" class TileLoader; @@ -50,7 +51,7 @@ private: QString _name; WMTS::Setup _setup; TileLoader *_tileLoader; - RectC _bounds; + RectD _bounds; QList _zooms; Projection _projection; Transform _transform;