-
Notifications
You must be signed in to change notification settings - Fork 71
Open
Description
Hi! I've a few questions about QGVProjection. I'm trying to create an EPSG3395 projection, but I'm having difficulty with maps. The projections in EPSG3758 and EPSG3395 have different projection points, but the same geographical projection point. I used formulas from https://wiki.openstreetmap.org/wiki/Mercator . Also, why is Y negative in the projection?
#include "QGVProjectionEPSG3395.h"
#include <QLineF>
#include <QtMath>
QGVProjectionEPSG3395::QGVProjectionEPSG3395()
: QGVProjection("EPSG3395",
"WGS84 World Mercator",
"QGVProjection used in web mapping applications like "
"Yandex/etc. Sometimes known as "
"EPSG:3395.")
{
mEarthRadius = 6378137.0;
mOriginShift = M_PI * mEarthRadius;
mEquatorialRadius = mEarthRadius;
mPolarRadius = 6356752.3142;
mEccentricity = qSqrt(1.0 - qPow(mPolarRadius / mEquatorialRadius, 2));
mCenterMass = mEccentricity / 2.0;
mGeoBoundary = QGV::GeoRect(85, -180, -85, +180);
mProjBoundary = geoToProj(mGeoBoundary);
}
QGV::GeoRect QGVProjectionEPSG3395::boundaryGeoRect() const
{
return mGeoBoundary;
}
QRectF QGVProjectionEPSG3395::boundaryProjRect() const
{
return mProjBoundary;
}
double QGVProjectionEPSG3395::degToRad(const double& deg) const
{
return deg * M_PI / 180.0;
}
QPointF QGVProjectionEPSG3395::geoToProj(const QGV::GeoPos& geoPos) const
{
const double lon = geoPos.longitude();
const double lat = (geoPos.latitude() > mGeoBoundary.topLeft().latitude()) ? mGeoBoundary.topLeft().latitude()
: geoPos.latitude();
const double x = mEquatorialRadius * degToRad(lon);
const double phi = degToRad(lat);
const double sinphi = qSin(phi);
double con = mEccentricity * sinphi;
con = qPow((1.0 - con) / (1.0 + con), mCenterMass);
const double ts = qTan(0.5 * (M_PI * 0.5 - phi)) / con;
const double y = mEquatorialRadius * qLn(ts);
return QPointF(x, y);
}
double QGVProjectionEPSG3395::radToDeg(const double& rad) const
{
return rad * 180.0 / M_PI;
}
QGV::GeoPos QGVProjectionEPSG3395::projToGeo(const QPointF& projPos) const
{
const double y = - projPos.y();
const double x = projPos.x();
const double lon = radToDeg(x) / mEquatorialRadius;
const double ts = qExp ( - y / mEquatorialRadius);
double phi = M_PI_2 - 2 * qAtan(ts);
double dphi = 1.0;
int i;
for (i = 0; fabs(dphi) > 0.000000001 && i < 15; i++) {
double con = mEccentricity * qSin (phi);
dphi = M_PI_2 - 2 * qAtan (ts * qPow((1.0 - con) / (1.0 + con), mCenterMass)) - phi;
phi += dphi;
}
const double lat = radToDeg(phi);
return QGV::GeoPos(lat, lon);
}
QRectF QGVProjectionEPSG3395::geoToProj(const QGV::GeoRect& geoRect) const
{
QRectF rect;
rect.setTopLeft(geoToProj(geoRect.topLeft()));
rect.setBottomRight(geoToProj(geoRect.bottomRight()));
return rect;
}
QGV::GeoRect QGVProjectionEPSG3395::projToGeo(const QRectF& projRect) const
{
return QGV::GeoRect(projToGeo(projRect.topLeft()), projToGeo(projRect.bottomRight()));
}
double QGVProjectionEPSG3395::geodesicMeters(const QPointF& projPos1, const QPointF& projPos2) const
{
const QGV::GeoPos geoPos1 = projToGeo(projPos1);
const QGV::GeoPos geoPos2 = projToGeo(projPos2);
const double latitudeArc = (geoPos1.latitude() - geoPos2.latitude()) * M_PI / 180.0;
const double longitudeArc = (geoPos1.longitude() - geoPos2.longitude()) * M_PI / 180.0;
const double latitudeH = qPow(sin(latitudeArc * 0.5), 2);
const double lontitudeH = qPow(sin(longitudeArc * 0.5), 2);
const double lonFactor = cos(geoPos1.latitude() * M_PI / 180.0) * cos(geoPos2.latitude() * M_PI / 180.0);
const double arcInRadians = 2.0 * asin(sqrt(latitudeH + lonFactor * lontitudeH));
return mEquatorialRadius * arcInRadians;
}
Metadata
Metadata
Assignees
Labels
No labels