Introduce EPSG3857Utils for accurate Web Mercator handling

[fdiazcarsi]

Related with issue #327.

This pull request introduces a new utility class, EPSG3857Utils, and refactors existing map bounds and tile caching logic to ensure accurate handling of the EPSG:3857 (Web Mercator) projection.

The EPSG:3857 projection presents distance and area distortions, which vary with latitude. Traditional linear calculations for map scale and bounding boxes are inaccurate in this projection, especially at higher latitudes. This change addresses these inaccuracies by centralizing and applying geodetic correction methods.

Key changes and their impact:

• New Utility: EPSG3857Utils.java: A new utility class has been added to encapsulate specific calculations required for EPSG:3857. This includes methods for:

• Identifying if a CRS is EPSG:3857.

• Computing a scaling factor to compensate for Mercator distortion.

• Calculating the true orthodromic width of a bounding box.

• Computing a ReferencedEnvelope from a center, scale, and paint area, accounting for distortion.

• WMTSLayer.java (WMTSTileCacheInfo constructor): Modified to use EPSG3857Utils.computeScalingFactor() when the map projection is EPSG:3857. This ensures that the imageBufferScaling correctly accounts for Mercator distortion, leading to more accurate tile rendering.

• BBoxMapBounds.java (getScale method): Refactored to use EPSG3857Utils.computeOrthodromicWidthInInches() for EPSG:3857 projections. This replaces a potentially inaccurate linear calculation with a geodetically correct one, ensuring the computed map scale is precise in the center of image.

• CenterScaleMapBounds.java (toReferencedEnvelope method): Updated to utilize EPSG3857Utils.computeReferencedEnvelope() for EPSG:3857. This ensures that the bounding box is accurately derived from the center and scale, correctly accounting for the projection's distortions.

These changes improve the accuracy of map rendering and scale calculations for EPSG:3857.

See JIRA issue: EPSG-3857.

[sbrunner]

Shouldn't we provide a generic solution?

[fdiazcarsi]

I don't quite understand what you're referring to with a "generic solution". The problem appears to be exclusively with the EPSG:3857 projection, or at least that's the one with which we have detected it.

The reason for the specific check for EPSG:3857 is due to the nature of this projection and how distances are calculated in it:

1. Geodetic Projections (units in degrees, like EPSG:4326):
   A degree of longitude does not represent a constant distance in meters; it becomes shorter as you move away from the equator. Therefore, you cannot simply convert the "width in degrees" to meters. A GeodeticCalculator must be used to calculate the orthodromic distance (the actual shortest distance over the curved surface of the Earth) between the left and right edges of the map.

2. "Normal" Planar Projections (units in meters, feet, etc.):
   In most projections (like UTM), the units are linear and have minimal distortion in the area of interest. In this case, it is safe to take the width of the BBox in its native units (e.g., meters) and convert it directly to inches for the scale calculation.

3. The Special Case: EPSG:3857 (Web Mercator):

   • The Problem: Although the units of EPSG:3857 are nominally "meters," they are not real meters on the Earth's surface. They are "Mercator meters," a measurement that stretches and distorts enormously as latitude increases.
   • The Consequence: If the code treated EPSG:3857 as a normal planar projection, it would use the width in "Mercator meters" to calculate the scale. This would give a completely incorrect result for any location not near the equator. For example, a map in Norway would have a very different calculated scale compared to its real scale.
   • The Solution: Just like with degree-based projections, the code cannot rely on the BBox width. Instead, it must delegate the calculation to a specialized function that calculates the actual orthodromic distance for a BBox in this particular projection.

The distinction between Geodetic and Planar projections was already being made in some methods, such as the getScale method in the BBoxMapBounds class. In the places where this distinction was already being made, we have introduced the necessary checks and calculations to manage the idiosyncrasy of EPSG:3857.

[sebr72]

@fdiazcarsi

1. We agree with improvement you are introducing.
2. We believe that the Normal case is no longer always the default case of choice (because of the new change you are introducing)
3. When looking at your code, It does not appear that the transformation is specific to 3857 (it is mentioned in if (EPSG3857Utils.is3857(crs)){ ) but the transformation itself does not reference 3857 directly. So we suggest that your solution does not reference 3857 directly but introduce a new parameter in the map part of requestData.json to say we want to use Geodetic transformation
4. We do not understand the difference between your option 1 (Geodetic Projections) and 3.The Special Case: EPSG:3857. To us, the seem identical. Please explain,

Finally, please add tests ;-) and rebase your code to clean the history.

[fdiazcarsi]

Allow me to try to explain the differences between the 3 cases to get to the point of why the Pseudo-Mercator projection is a special case:

• Case 1 (Degrees, geographic):
  as for example: https://spatialreference.org/ref/epsg/4326/

  • Type: GEOGRAPHIC_2D_CRS
  • WGS84 Bounds: -180.0, -90.0, 180.0, 90.0
  • UoM: degree

  It is not a projection in the strict sense; it works with angles. This implies using the geodetic calculator (org.geotools.referencing.GeodeticCalculator) to calculate the orthodromic distance between two points. Based on the latitude and longitude values in degrees of those points, the shortest distance between them on the globe is calculated.

• Case 2 (linear units, projected):
  as for example, https://spatialreference.org/ref/epsg/25830/

  • Type: PROJECTED_CRS
  • WGS84 Bounds: -6.01, 35.26, 0.0, 80.49
  • UoM: metre

  It is projected, works in metric units, and is defined for a specific area with specific parameters to minimize distortions within its area of validity.
  To calculate the distance between two points within its area of validity, it is sufficient to use Euclidean geometry.

• Case 3: Pseudo-Mercator
  as for example, https://spatialreference.org/ref/epsg/3857/

  • Type: PROJECTED_CRS
  • WGS84 Bounds: -180.0, -85.06, 180.0, 85.06
  • UoM: metre

  It is projected, so it works in meters, but it is defined for the entire globe.
  Although it is conformal for angles, it is not for distances or areas, so the difference between the real world and the projected one increases significantly as we move away from the equator and approach the poles.
  Since it works in meters, it might seem that using Euclidean geometry would be sufficient to calculate the distance between two points. However, the distances between two points near the poles would yield a much larger result than the distance between two points near the equator, even if they were the same real-world distance apart.
  We cannot use the same method already used for Case 1 (e.g., computeGeodeticBBox in BBoxMapBounds) because it works with degrees (latitude/longitude), and this projection is in meters. Therefore, this case is very special and must be handled on its own path. Although at first glance the code in the new class's methods is very similar to what already exists, upon closer inspection, there are very subtle but important differences. We have always preferred to respect the code that already works, and we assume correctly, to avoid introducing potential errors.

I fully understand the concern about changing the default behavior and trying to find a generic solution that doesn't depend on a specific if for "EPSG:3857". Backward compatibility and a clean design are very important.

However, we would like to insist that the current handling of Pseudo-Mercator projections is not a simple inaccuracy, but a fundamental error in calculation. The result it produces for maps far from the equator (such as Greenland) is not a matter of preference, but is objectively incorrect and can lead to the creation of maps with completely misleading scales. We are aware that the calculations we propose do not fix the problem, but they do mitigate it significantly. For maps that are not very wide in latitude, the mitigation can even solve it, whereas for maps that are wide in latitude, differences in measurements at the edges will still be observed when compared to the scale bar (since the measurements would only be correct at the center of the map, and they lose precision as they move north or south).

The problem is tied to the Pseudo-Mercator projection with the "WGS 84" datum, but right now I only know of one non-deprecated CRS associated with it, "EPSG:3857". In fact, in the EPSG3857Utils.is3857 method, the check performed is whether it is "WGS 84 / Pseudo-Mercator" (perhaps we should introduce checks for "WGS 84 / Web Mercator", which is an alias we might still encounter). We could generalize it to that datum with that projection and change the class name to PseudoMercatorUtils and the method to isPseudoMercator. If in the future there are other CRSs associated with that projection, we could handle them in the same place. The calculations performed in the utility class should be common to any CRS with that datum and projection. I understand that the class and method names were probably not entirely well-chosen and could be misleading.

I also understand that introducing the change to improve measurement accuracy on the map when using that specific CRS may cause users who are currently generating maps to suddenly get different maps, and this could cause them some problems. I imagine that's why the suggestion was made to introduce something like useGeodeticCalculations in the "map". It's a safeguard so that what is already in use continues to work as it has until now (even if it's not correct). In the three places where we have now introduced special handling for CRSs based on "WGS 84 / Pseudo-Mercator", and there is an if( EPSG3857Utils.is3857(crs) )..., the check with useGeodeticCalculations could also be added: if( useGeodeticCalculations && PseudoMercatorUtils.isPseudoMercator(crs) )....

We are going to spend some time trying to introduce the useGeodeticCalculations parameter into the "map", but I'm not sure how to do it. I suppose we will have to consider that it can be in the YAML, in the JSON, and figure out in which Java class it should end up "living" and how it can be passed to the points where we need it. We would appreciate any hints on this last point.

[sbrunner]

It should be in the spec.json. The first thing to do is to add the parameter on this class https://github.com/mapfish/mapfish-print/blob/master/core/src/main/java/org/mapfish/print/attribute/map/GenericMapAttribute.java#L327 he will be automatically filed.
In the CreateMapProcessor.execute he will be available on param.map.geodetic.

Can you also clarify the @sebr72 point 4 (differences between cases 1 and 3 except for the additional reprojection)?

[fdiazcarsi]

Regarding the clarification of point 4, about the "differences between cases 1 and 3 except for the additional reprojection," we have been reviewing the code that we have introduced in:

• Class BBoxMapBounds, method getScale
• Class "WMTSLayer.WMTSTileCacheInfo", in the constructor.
• Class CenterScaleMapBounds, method toReferencedEnvelope (computeGeodeticBBox)

We will discuss each one.

BBoxMapBounds, getScale method

We proposed something like:

     ...
     if (projUnit == DistanceUnit.DEGREES) {
       GeodeticCalculator calculator = new GeodeticCalculator(getProjection());
       final double centerY = bboxAdjustedToScreen.centre().y;
       calculator.setStartingGeographicPoint(bboxAdjustedToScreen.getMinX(), centerY);
       calculator.setDestinationGeographicPoint(bboxAdjustedToScreen.getMaxX(), centerY);
       DistanceUnit.fromString(calculator.getEllipsoid().getAxisUnit().toString());
 
       geoWidthInInches = ellipsoidUnit.convertTo(geoWidthInEllipsoidUnits, DistanceUnit.IN);
    } else if (EPSG3857Utils.is3857(crs)){
        geoWidthInInches = EPSG3857Utils.computeOrthodromicWidthInInches(bboxAdjustedToScreen);
   ...

Reviewing it more carefully, the code in EPSG3857Utils.computeOrthodromicWidthInInches is similar to the one used for calculations in degrees, except that a transformation is performed before the calculations. In this case, we could unify the two cases by adding a new method, computeGeodeticWidthInInches, to the BBoxMapBounds class with something like:

private double computeGeodeticWidthInInches(final ReferencedEnvelope bbox) {
    try {
        CoordinateReferenceSystem crs = bbox.getCoordinateReferenceSystem();
        GeodeticCalculator calculator = new GeodeticCalculator(crs); // Use the original CRS for the ellipsoid

        Coordinate startDegrees;
        Coordinate endDegrees;

        if (DistanceUnit.fromProjection(crs) == DistanceUnit.DEGREES) {
            // Already in degrees, use directly
            startDegrees = new Coordinate(bbox.getMinX(), bbox.centre().y);
            endDegrees = new Coordinate(bbox.getMaxX(), bbox.centre().y);
        } else { // Assuming it is is3857() from the if statement above
            // Needs reprojection
            final MathTransform transform = CRS.findMathTransform(crs, GenericMapAttribute.parseProjection("EPSG:4326", true));
            startDegrees = JTS.transform(new Coordinate(bbox.getMinX(), bbox.centre().y), null, transform);
            endDegrees = JTS.transform(new Coordinate(bbox.getMaxX(), bbox.centre().y), null, transform);
        }

        // --- Common Logic ---
        calculator.setStartingGeographicPoint(startDegrees.x, startDegrees.y);
        calculator.setDestinationGeographicPoint(endDegrees.x, endDegrees.y);
        final double orthodromicWidth = calculator.getOrthodromicDistance();
        final DistanceUnit ellipsoidUnit = DistanceUnit.fromString(calculator.getEllipsoid().getAxisUnit().toString());

        return ellipsoidUnit.convertTo(orthodromicWidth, DistanceUnit.IN);

    } catch (Exception e) {
        throw new PrintException("Failed to compute geodetic width", e);
    }
}

And leaving in getScale something like:

public Scale getScale(final Rectangle paintArea, final double dpi) {
    final ReferencedEnvelope bboxAdjustedToScreen = toReferencedEnvelope(paintArea);
    DistanceUnit projUnit = DistanceUnit.fromProjection(getProjection());
    
    double geoWidthInInches;
    CoordinateReferenceSystem crs = getProjection();

    // If it is geodetic OR it is a special case requiring geodetic calculation
    if (projUnit == DistanceUnit.DEGREES || EPSG3857Utils.is3857(crs)) {
        geoWidthInInches = computeGeodeticWidthInInches(bboxAdjustedToScreen);
    } else {
        // Original logic for planar projections
        geoWidthInInches = projUnit.convertTo(bboxAdjustedToScreen.getWidth(), DistanceUnit.IN);
    }

    // ... rest of the getScale method ...
}

This way, it would be integrated into the BBoxMapBounds class, and we would remove the computeOrthodromicWidthInInches method from EPSG3857Utils.

Constructor of "WMTSLayer.WMTSTileCacheInfo"

In this constructor, an "extra" scale factor is simply being calculated to apply to the image when dealing with the "Pseudo-Mercator" projection. To calculate it, the EPSG3857Utils.computeScalingFactor method is called. We can move that method as a private method into the WMTSTileCacheInfo class so that its scope is limited to where it is used; however, the logic of the method would remain the same.

CenterScaleMapBounds, toReferencedEnvelope method

This is the most conflicting case.

In the CenterScaleMapBounds class, we have the computeGeodeticBBox method, which works with CRS in degrees (e.g., EPSG:4326), whereas in EPSG3857Utils we have introduced the computeReferencedEnvelope method, which works on a Pseudo-Mercator projection in meters. The necessary calculations for each are essentially different. The main differences would be:

1. Different Concepts of Scale and Input Parameters

The two methods are triggered by different scale definitions, which determines their starting units.

• Case 1 (Geographic, e.g., EPSG:4326): Starts from a physical scale linked to the output medium (DPI). The calculation begins with this.scale.getResolutionInInches() * paintArea.width, which results in a physical dimension in inches. The objective is to determine the geographic Bounding Box that corresponds to a specific physical size on the screen or on paper.

• Case 3 (Pseudo-Mercator): Starts from an abstract cartographic scale. The calculation begins with scale.getResolution() * paintArea.width, which results in a world dimension in meters. The objective is to determine the projected Bounding Box that corresponds to a given map scale (e.g., 1:1,000,000).

This initial difference in intent (physical scale vs. cartographic scale) means we cannot simply reuse the first method, as the input parameters are conceptually and numerically different.

2. Fundamentally Different Algorithm for the Y-Axis

This is the most critical algorithmic difference.

• Case 1 (computeGeodeticBBox): Calculates the north and south boundaries by asking the GeodeticCalculator for a destination coordinate:

  calc.setDirection(north, geoHeight / 2.0);
  double maxGeoY = calc.getDestinationPosition().getOrdinate(1); // Returns a latitude in degrees

  This works because it operates in an angular space (based on degrees).

• Case 3 (our computeReferencedEnvelope): Cannot do this. Due to the extreme distortion of the Y-axis in Pseudo-Mercator, we must use a different approach. We calculate the actual orthodromic distance and then apply it as a linear offset to the Y-coordinate (in Pseudo-Mercator meters) of the center:

  calc.setDirection(north, geoHeight / 2.0);
  double northHeight = calc.getOrthodromicDistance(); // Returns a distance in meters
  double maxGeoY = calc.getStartingPosition().getOrdinate(1) + halfHeight; // Arithmetic operation on Mercator meters

  This is a specific correction for the properties of the Pseudo-Mercator projection. Using the getDestinationPosition() logic here would produce an incorrect result.

3. Different Final BBox Construction (Output Processing)

Finally, the construction of the ReferencedEnvelope is different because the coordinate systems behave differently.

• Case 1 (Geographic): Must use rollLongitude() and rollLatitude() on the final coordinates. This is essential for handling wrapping around the globe (e.g., normalizing 190° longitude to -170°). It is a mandatory step for angular and spherical coordinate systems.

• Case 3 (Pseudo-Mercator): Operates in a Cartesian and planar space where coordinates are linear meters. There is no "wrap-around" at 180 meters. Applying a function like rollLongitude to these coordinates in meters would be meaningless and would corrupt the Bounding Box.

In summary:

The two methods differ in their initial inputs (physical vs. cartographic scale), their main algorithm for the Y-axis (geodetic destination vs. correction with a linear offset), and their final processing of the result (angular normalization vs. direct planar coordinates).

For these reasons, we conclude that creating a dedicated and specialized method for Pseudo-Mercator is the clearest and safest design choice. It correctly encapsulates the unique handling that this projection requires and avoids complicating the existing and proven computeGeodeticBBox method with conditional logic that would have to account for all these differences.

We could integrate the method we created, computeReferencedEnvelope, into the CenterScaleMapBounds class alongside computeGeodeticBBox to have the logic for this operation in the same class. It could be named computeGeodeticBBoxInPseudoMercator.

---

I don't know if this resolves the doubts about point 4, "differences between cases 1 and 3 except for the additional reprojection," or if they were related to something else.

On the other hand, if we make these changes to integrate the Pseudo-Mercator related calculations where they are needed, only the is3857 method would remain in the EPSG3857Utils utility class (which would be renamed to PseudoMercatorUtils and isPseudoMercator). I am not sure if it would be advisable to keep this class or move the isPseudoMercator method to another utility class that already exists in the project.

Regarding the introduction of the new parameter in param.map, we are looking into how to do it, and we would add the corresponding checks.

[sebr72]

@fdiazcarsi Thanks for all the explanations. We will get back to you shortly.

[fdiazcarsi]

We have made the changes we proposed to integrate our code with the project's base code. We added the useGeodeticCalculations attribute to the map parameters, added unit tests for the new methods, and included an example of how to use this attribute (use_geodetic_calculations_EPSG_3857).

We also performed the rebase and verified that the example works correctly.