Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

EPSG 3395 Map shift #60

Open
MatveiStudentIITK opened this issue Feb 4, 2025 · 0 comments
Open

EPSG 3395 Map shift #60

MatveiStudentIITK opened this issue Feb 4, 2025 · 0 comments

Comments

@MatveiStudentIITK
Copy link

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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant