[
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)
