[jira] [Commented] (SIS-547) Mercator projection should wraparound longitude values

["Martin Desruisseaux (Jira)" <ji...@apache.org>]

    [ https://issues.apache.org/jira/browse/SIS-547?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17575799#comment-17575799 ] 

Martin Desruisseaux commented on SIS-547:
-----------------------------------------

Below are some reply to email, copied here in the hope to clarify.

{quote}
(…snip…) for a local projection and for local geometries the consensus seems to be that eastings or longitudes should be continuous.
{quote}

Yes, keeping eastings and longitudes continuous is the difficulty. For achieving this goal, sometime a wraparound is necessary and sometime it is a disturbance. The difficulty is to find rules for deciding when to apply the wrapapound.

For some map projections, the longitude value is used directly in a sine or cosine function. In those lucky cases we don't need to care. Not only the wraparound is applied automatically by sine/cosine, but in addition the result is smooth (e.g. no sudden jump from non-zero positive eastings to negative eastings). This is the case of "Polar Stereographic", "Orthographic", "Lambert Azimuthal Equal Area", "Azimuthal Equidistant (Spherical)" and "Modified Azimuthal Equidistant" projections among others. For those map projections, the problem is moot.

But for some other map projections, the longitude values is not used in any trigonometric functions. It is simply multiplied by some factor. It is the case of "Mercator", "Cylindrical Equal Area", "Sinusoidal equal-area" and "Mollweide" projections among others. For those map projections, no wraparound occurs automagically; we have to do it explicitly. But the problem is that contrarily to the case in above paragraph (which goes smoothly through zero), applying a wraparound in this case causes a strong discontinuity in the mathematical function, where easting jumps suddenly from a large positive value to a large negative value. This is a situation that we try to avoid as much as possible, which is the cause of the difficulty mentioned in the first paragraph.

There is a third, more pernicious case where the longitude value is used in a sine or cosine function, but not directly. Instead the longitude is multiplied by some factor before to be used in the sine/cosine function. In that case a wraparound is applied by the sine/cosine, but that wraparound is not 360° of longitude. Instead the wraparound occurs with a period which depend on the multiplication factor _n_ in expressions like {{sin(n*λ)}}. It happens in "Albers Equal Area", "Oblique Mercator", "Oblique Stereographic" and "Lambert Conic Conformal" projections among  others. In those cases, a second 360° wraparound may need to be applied prior the wraparound performed by the sine/cosine functions.

If I remember correctly, the non-360° wraparound occurring with Lambert Conic Conformal was the cause of a transformation error reported a while ago in Alaska. The fix has been to apply a 360° wraparound prior the Lambert projection, but in a way as restrictive as possible. Instead of applying an unconditional 360° wraparound, the rules which have been set for the Lambert projection is as below:

* If the projection central meridian is 0°, do not apply any wraparound.
* If the projection central meridian is positive, apply a 360° wraparound only on longitude less than -180° + λ₀.
* If the projection central meridian is negative, apply a 360° wraparound only on longitude greater than 180° + λ₀.

The rational for the last two points is that the central meridian is subtracted from the longitude value before to be used in sine/cosine functions. So if the central meridian is -100° (negative), a longitude value of 100° may become 200°, which explain the "longitude greater than 180° + λ₀" condition.

Above set of rules apparently worked well for Lambert projection; I have not hear any issue since they have been established. But in the case of the Mercator projection, the fix attempt uses a different rule: it just applies an unconditional 360° wraparound for any longitude outside the -180°…+180° range, and regardless the projection central meridian. This policy caused a regression in codes that rely on a continuous behaviour for envelope or geometries projections.

After more thoughts, my current inclination is to generalize to all projections (except when unnecessary by the virtue of sine/cosine) the set of rules used for the Lambert projection. It should resolve the issue reported for "Mercator 41" (because its central meridian is 100°E, a wraparound will be applied on longitudes less than -80°) while avoiding the regression reported above (because its central meridian is 0, no wraparound is applied).

{quote}
 _[In the context of "Mercator 41", which has a central meridian at 100°E]_ The point  being that if you move increasingly eastward on the earth (e.g., from 179E to 179W) then the Easting should increase.
{quote}

It should have been fixed on 1.3-SNAPSHOT, but using the unconditional wraparound and for the forward projection only, not the inverse one. If the bug is still present on the forward projection, we will need to investigate why. But given that I would like to amend the current behaviour with above proposal (replace unconditional wraparound by more restrictive criteria), it may be better to re-test after that change. With above criteria, the issue described above should be okay, because 179°W is less than -80°, so it should be converted to 181°E before the "Mercator 41°" projection is applied. For the reverse projection, we need to define a set of rules as well.

{quote}
Input of 181E should be accepted by API and give the same XY output as 179W.  (for all f(LL)=XY=f(LL’) where L’ longitudes are modulo 360 of L).
{quote}

Well, this is the difficulty: input of 181°E should be accepted, but whether it should give the same XY than 179°W depends on what we want to compute. If we are projecting a geometry or if we are resampling the pixels of a raster, the sudden jump from large positive X to large negative X with Mercator projection is a problem. In remote sensing for example, it is very common to have rasters with longitude values crossing the 180° limit for some pixels. So sometime we want same XY, sometime not. It depends on the context, which is why the problem is hard.

The "less bad" approach I have found so far is to try to keep the mathematical functions as continuous as possible, even at the cost of XY values outside their normal range if the (φ,λ) input is outside the normal domain. This approach implies avoiding as much as possible the discontinuities (i.e. sudden jumps in XY values) caused by wraparounds. However wraparound is necessary when the projection central meridian is not zero, so the set of rules described above for the Lambert projection is the compromise I have found so far.

A rational for this approach is that, if an unconditional wraparound is desired, the user of Apache SIS API can apply the wraparound himself on the (φ,λ) inputs before to invoke a transform(…) method. But if a wraparound is not desired, for example for resampling a raster crossing the 180° meridian, then it would be very difficult for the user to remove the wraparound effect if it was applied unconditionally by the projection code.

{quote}
(…snip…) for geometry longitude in degrees needs to be continuous.   (..so longitude 178, 179, 180 then 181, 182, etc.(and same around -180).)
{quote}

Yes, this is exactly the problem. As said above, sometime we need 181° to be replaced by -179°, and sometime not. It depends on the context. Basically the golden rule that I'm trying to follow is: _Keep the (x,y) = P(φ,λ) mathematical function as continuous as possible._ It usually implies avoiding wraparound. The case where the projection central meridian is different than zero is an exception, because it is a case where applying the wraparound can actually remove the discontinuity caused by points crossing the ±180° line. But except for that case (which is the case that I need to fix in the same way for all map projections), I'm inclined to favour mathematical continuity, because it is easier for users to apply wraparound themselves than to remove discontinuities caused by undesired wraparounds.


> Mercator projection should wraparound longitude values
> ------------------------------------------------------
>
>                 Key: SIS-547
>                 URL: https://issues.apache.org/jira/browse/SIS-547
>             Project: Spatial Information Systems
>          Issue Type: Improvement
>          Components: Referencing
>    Affects Versions: 0.6, 0.7, 0.8, 1.0, 1.1, 1.2
>            Reporter: Martin Desruisseaux
>            Assignee: Martin Desruisseaux
>            Priority: Major
>             Fix For: 1.3
>
>         Attachments: image-2022-07-29-17-16-59-619.png, image-2022-07-29-17-21-47-669.png
>
>
> When the central meridian has a value different than zero, user may expect a range of 180° on the east of central meridian to produce positive easting values and conversely on for 180° on the west of central meridian. This is currently not the case:
> {code:java}
> ProjectedCRS crs = (ProjectedCRS) CRS.forCode("EPSG:3994");
> MathTransform prj = crs.getConversionFromBase().getMathTransform();
> double[] coordinates = {-41, 100, -41, 179, -41, 181, -41, -179};
> prj.transform(coordinates, 0, coordinates, 0, 4);
> System.out.println(java.util.Arrays.toString(coordinates));
> {code}
> Output (reformatted for readability):
> {code:none}
>         0   -3767132
>   6646680   -3767132
>   6814950   -3767132
> -23473717   -3767132
> {code}
> Other map projection implementations rely on trigonometric functions for applying an implicit wraparound. For example in a call to {{sin(λ)}} the λ argument value is implicitly reduced to a range of -π … +π around the λ₀ (central meridian). But it does not happen in the particular case of the Mercator projection, since the Easting value is just a multiplication factor without trigonometric functions involved. So we have to do the wraparound ourselves.



--
This message was sent by Atlassian Jira
(v8.20.10#820010)