Global image rotation and placement

[PaulWessel]

Description of the problem

This issue also relates to the problems discussed in #4335.

When we read in a global grid (which may have a longitude range such as -180/180 or 0/360, typically) we are aware of the w/e/s/n of the desired plot area. If that area is a subset then we of course only read that subset. If, however, the w/e/s/n is the global domain but with a different central meridian, then we roll each column in the grid so that upon receipt from the file the grid now has the same w/e boundary as the desired plot [the fact that we do not read a file is one of the aspects that makes #4335 an issue].

However, it should be clear that the rolling of longitudes can only be done in whole steps of longitude increments. For instance, for a global grid with 1x1 degree resolution, we can roll a -180/180 range to 1/361, but surely we can never actually roll it to become 0.5/360.5 since half grid-nodes do not exist. Yet, because we have mostly been working with higher resolution data sets, the issues with fractional node rolling have largely gone unnoticed, Here, I document the effects by making some -JQ (default projection) global maps of @earth_relief_01d.grd (pixel grid), plus drawing 1x1 degree gridlines on top. Note: You may need to zoom in on these small samples to see the details.

Using -R0/360/-90/90:

Zooming in on the left (Greenwich) boundary near the UK, I show what the image looks like. This is correct since the new central meridian (or left boundary if you like) is a multiple of 1 degree:

Using -R1/361/-90/90:

Stepping a full degree is not a problem since that is the grid resolution:

Using -R0.5/360.5/-90/90:

We cannot roll by half a degree but we can of course plot. Grid is still global but something goes wrong:

As you can see, the first column from 0-1 longitude, which should be halfway clipped by the west boundary, is entirely missing, leaving a gap of width 0.5 degrees. Let's look at the east boundary:

All is well there: The "rolled" column that was 0-1 is now 360-361 and it is imaged, with half clipped by the east boundary.

In the fractional roll case, clearly the affected column is required to appear twice (with fractional clipping at either end), yet there is no repeated column in the pixel-registered grid. If I instead made that same map with the gridline-registered @earth_relief_01d_g.grd then it all works fine (at either boundary):

Conclusions:

When rolling global pixel-registered grids by a fractional node amount, the repeated column required by the periodic boundary for many map projections is naturally missing. Options to fix this problem are:

1. Warn users to use gridline-registered grids.
2. Detect case internally during reading and call grdsample -T to get equivalent gridline-registered version
3. Special handling by plotting the thin missing strip after the main image.

Each of these has challenges and problems:

1. Clearly, this is not helpful to users who as a data set they want to image and it happens to be pixel-registered. For our own remote data we could perhaps issue a message and switch from the default pixel to gridline version, if it exists, but for users grids, not so much.
2. When reading from file it is possible this could be fixed by adding a repeating column to a pixel grid (I am not sure if this runs afoul of other checks in GMT since we do not tolerate grids with > 360 extent and this would be 360 + xinc. If that worked then perhaps it would be OK. However, for read-only grids as in Matrix grid doesn't always work with different central meridians #4335 it would not help unless we make a copy.
3. This is very messy as it would apply to grdimage, grdview, grdcontour for sure, perhaps more. A per-module solution is not a good solution.

A final option is to basically ignore this except when the grid is course enough so that the plot width of one column exceeds some threshold, say one point or whatever. This means we limit the problems to courser grids. Also note I am pretty sure this problem is also only applicable to global plot domain with periodic map boundaries, so not every map case will fall into the bad category.
Also, it may be an opportunity for us to yet again curse the notion of pixel grids and reconsider our selection of the remote pixel grids as the default (we have waffled on this before). We are basically causing this problem for ourselves, but perhaps we should be happy to at least detect the problem. Yet, we cannot curse the notion of pixel images...

So there you have it. At this point, feedback from @GenericMappingTools/core would be helpful.

[WalterHFSmith]

Seems to me we shouldn’t allow someone to roll a fraction of a delta_x sample. Walter H F Smith Geophysicist, Laboratory for Satellite Altimetry NOAA NCWCP code E/RA-31 5830 University Research Court, room 3752 College Park MD 20740-3818 301-683-3377 (voice, my desk) 301-683-3300 (voice, reception) 301-683-3301 (fax, reception) [email protected] http://www.star.nesdis.noaa.gov/star/Smith_WHF.php

…

On Nov 10, 2020, at 3:07 PM, Paul Wessel ***@***.***> wrote: Description of the problem This issue also relates to the problems discussed in #4335. When we read in a global grid (which may have a longitude range such as -180/180 or 0/360, typically) we are aware of the w/e/s/n of the desired plot area. If that area is a subset then we of course only read that subset. If, however, the w/e/s/n is the global domain but with a different central meridian, then we roll each column in the grid so that upon receipt from the file the grid now has the same w/e boundary as the desired plot [the fact that we do not read a file is one of the aspects that makes #4335 an issue]. However, it should be clear that the rolling of longitudes can only be done in whole steps of longitude increments. For instance, for a global grid with 1x1 degree resolution, we can roll a -180/180 range to 1/361, but surely we can never actually roll it to become 0.5/360.5 since half grid-nodes do not exist. Yet, because we have mostly been working with higher resolution data sets, the issues with fractional node rolling have largely gone unnoticed, Here, I document the effects by making some -JQ (default projection) global maps of @earth_relief_01d.grd (pixel grid), plus drawing 1x1 degree gridlines on top. Note: You may need to zoom in on these small samples to see the details. Using -R0/360/-90/90: Zooming in on the left (Greenwich) boundary near the UK, I show what the image looks like. This is correct since the new central meridian (or left boundary if you like) is a multiple of 1 degree: Using -R1/361/-90/90: Stepping a full degree is not a problem since that is the grid resolution: Using -R0.5/360.5/-90/90: We cannot roll by half a degree but we can of course plot. Grid is still global but something goes wrong: As you can see, the first column from 0-1 longitude, which should be halfway clipped by the west boundary, is entirely missing, leaving a gap of width 0.5 degrees. Let's look at the east boundary: All is well there: The "rolled" column that was 0-1 is now 360-361 and it is imaged, with half clipped by the east boundary. In the fractional roll case, clearly the affected column is required to appear twice (with fractional clipping at either end), yet there is no repeated column in the pixel-registered grid. If I instead made that same map with the gridline-registered @earth_relief_01d_g.grd then it all works fine (at either boundary): Conclusions: When rolling global pixel-registered grids by a fractional node amount, the repeated column required by the periodic boundary for many map projections is naturally missing. Options to fix this problem are: • Warn users to use gridline-registered grids. • Detect case internally during reading and call grdsample -T to get equivalent gridline-registered version • Special handling by plotting the thin missing strip after the main image. Each of these has challenges and problems: • Clearly, this is not helpful to users who as a data set they want to image and it happens to be pixel-registered. For our own remote data we could perhaps issue a message and switch from the default pixel to gridline version, if it exists, but for users grids, not so much. • When reading from file it is possible this could be fixed by adding a repeating column to a pixel grid (I am not sure if this runs afoul of other checks in GMT since we do not tolerate grids with > 360 extent and this would be 360 + xinc. If that worked then perhaps it would be OK. However, for read-only grids as in #4335 it would not help unless we make a copy. • This is very messy as it would apply to grdimage, grdview, grdcontour for sure, perhaps more. A per-module solution is not a good solution. A final option is to basically ignore this except when the grid is course enough so that the plot width of one column exceeds some threshold, say one point or whatever. This means we limit the problems to courser grids. Also note I am pretty sure this problem is also only applicable to global plot domain with periodic map boundaries, so not every map case will fall into the bad category. Also, it may be an opportunity for us to yet again curse the notion of pixel grids and reconsider our selection of the remote pixel grids as the default (we have waffled on this before). We are basically causing this problem for ourselves, but perhaps we should be happy to at least detect the problem. Yet, we cannot curse the notion of pixel images... So there you have it. At this point, feedback from @GenericMappingTools/core would be helpful. — You are receiving this because you are on a team that was mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[PaulWessel]

That is indeed a possible solution: Limit a global plot w/e to a multiple of grid resolutions for global maps with periodic boundaries. Might be awkward for truly curse grids (5x5) that we cannot pick any -R. But way easier to implement than any other options.

[PaulWessel]

Actually, I should be clear on further limitations on this problems. It is only an issue when all these are true:

1. Pixel-registered grid (and possibly image)
2. Full 360-degree range plotting
3. Projection with repeated periodic boundaries (e.g., cylindrical)
4. Projection of grid not required as graticules are linear and orthogonal
5. Users did not request interpolation via -E.

I did not stress conditions 4 (or 5) above. Once you pick a projection where we must reproject the lon/lat lattice onto a new x/y lattice then we are again fine.

So it seems this problem is limited to -JQ only ( the "longlat" projection), in which case Walter's suggestion is appealing.