Implement auto resolution for remote grids for plotting modules

[PaulWessel]

Description of proposed changes

So far, the user has needed to specify which resolution they want for a remote grid. Unfortunately, knowing what is a good resolution when the goal is to make a map that will be 15 cm is not straightforward, so users tend to select too high resolutions given the size of their plots, resulting in extra download and processing time, plus bloated plots. This PR implements a feature for plotting where the resolution of a remote grid is optional. If so, we compute the resolution that is most compatible with an anticipated dots-per-unit of the plotted image. For non-plotting purposes (e.g., grdmath, grdproject, grdsample, ...) the user must specify the resolution. The key benefit of this PR is the huge simplification in making maps that uses the GMT remote grids for backgrounds in that no matter what the plot region, projection, and size is selected, the result will look good.

I have implemented a simple algorithm for this. Given -R -J I can compute two different representative lengths from the plot:

• A representative length D in spherical degrees (e.g., length of the diagonal or for some maps just the longitude extent)
• A representative length L in inches (e.g., length of the plot diagonal or for some maps just the width)

Then, given the desired dots-per-unit (new GMT default GMT_GRAPHICS_DPU [300i) we seek to find a grid whose number of nodes per degree, n, is the closest to satisfying this equation: L * dpi = n * D.

If the user has a reason to prefer a particular grid registration then they can still append _g or _p, but if no such version exist it will generate an error as no files can be found.

Here is a few examples of such a plots:

gmt grdimage @earth_relief -R0/60/0/60 -JM15c -B -pdf map
gmt grdimage @earth_relief -R270/20/305/25+r -JOc280/25.5/22/69/12c -B -pdf map
gmt grdimage @earth_age -Rg -JH0/15c -B -png t
gmt grdimage @earth_day -Rg -JH0/15c -B -png t

etc. Please try to break it. For now I am printing out an information line so we can see how it decides, e.g.

gmt grdimage @earth_age -Rg -JH0/15c -B -png t
gmt [NOTICE]: Given -Rg -JH0/15c, D_g = 360 degrees and D_m = 5.90551 inches and given dpu = 300i so we replace @earth_age by @earth_age_15m_p

I encourage you to play with the regions, plot size, and desired DPU to see how it affects the result, and let me know if you cind some combination that works poorly.

Note if I do

gmt grdimage @earth_night_g -Rg -JH0/15c -B -png t

it will give an error since there are no gridline version of these images...

Note: In causing this to fail I noticed that when a one-liner code fails in gmt_init_module and returns NULL, there is no place to terminate the modern session and the next modern one-liner will say

gmt [ERROR]: Cannot run a one-liner modern command within an existing modern mode session

The solution to that is to use the new terminate function I wrote for #5739 which is not completed, so this will go away later.

[PaulWessel]

Once this is tested and approved we should update all examples where we make plots with remote files to make it simpler.

[Esteban82]

Very nice feature!!

I already have an script with a zoom effect. I tested it. It looks good. The grids from 15m to 02m were used.

@PaulWessel with this script I got the message I got the warning message of the forum.

#!/usr/bin/env bash
cat << EOF > pre.sh
gmt begin
	gmt math -T3200/7900/50 -o1 T -I = times.txt
	gmt makecpt -Cgeo -H > topo.cpt
gmt end
EOF
cat << EOF > main.sh
gmt begin
	gmt grdimage @earth_relief -Rd -JG-58:26/-34:36/\${MOVIE_COL0}/0/0/0/60/60/14.9c -Yc -Xc -B0 -Ctopo.cpt
gmt end
EOF
gmt movie main.sh -Sbpre.sh -NZoom_Relief -Ttimes.txt -C15cx15cx100 -D24 -M0,png -Fmp4

Zoom_Relief.mp4

[PaulWessel]

Thanks, I fixed the issue that made the lower altitudes take a very long time (specially 500 and up). Now I can do 500/8000 but when I get even closer and need 03s and eventually 01s I have an issue of not finding (?) or downloading tiles. So more testing. I think making anim16 doing a zoom from "outer space" all the way down to 1 s resolution would be nice so your script is a good start.

[Esteban82]

I try to run the same script but with @earth_day and and it tries to download earth_day_30s_p (which doesn't exist, right?)

Error message:

gmt [NOTICE]: Given -Rd -JG-58:26/-34:36/1500/0/0/0/60/60/13.4ch, D_g = 21.7051 degrees and D_m = 7.46081 inches and given dpu = 300i so we replace @earth_day by @earth_day_30s_p
movie_master.sh: línea 13: 169150 Terminado (killed)      gmt grdimage @earth_day -Rd -JG-58:26/-34:36/${MOVIE_COL0}/0/0/0/60/60/13.4ch -Yc -Xc -B0

[Esteban82]

The first frames with earth_day look blurry.

Zoom_day2.mp4

#!/usr/bin/env bash
cat << EOF > pre.sh
gmt begin
	gmt math -T3000/10000/100 -o1 T -I = Altitudes.txt
gmt end
EOF
cat << EOF > main.sh
gmt begin
	gmt grdimage @earth_day -Rd -JG-58:26/-34:36/\${MOVIE_COL0}/0/0/0/60/60/13.4ch -Yc -Xc -B0
gmt end
EOF
gmt movie main.sh -Sbpre.sh -NZoom_day3 -TAltitudes.txt -C13.5cx13.5cx100 -Ml,png

[anbj]

This is a fantastic feature.
All those times a good plot has been destroyed by too many details (..).

Does this mean that, for plotting purposes, one is discouraged to set the resolution, and instead let gmt decide?

Also, I would like to emphasize how great the remote grids and auto download that gmt offers is and has been for me; it substantially reduces the friction of making (good) maps.
If there are other public datasets that would be a nice addition to the remote datasets offered by gmt, I would highly encourage it to be included.
In my opinion, this is something that reduces the bar a great deal, and attracts new users. If you eliminate the factor of having to harvest datasets (at least, some), that is a great advantage. Making a good map when you have all the data is hard enough!

[PaulWessel]

Once we are happy with how it works, I would say definitively. I doubt most GMT users could figure out a suitable resolution so they use too high a resolution to make sure it looks good. But yes, there may be times when you

• Require a specific resolution for other reasons
• Intend to zoom in on the figure so you want it to have more detail that can be exposed by zooming (PDF)
• Want a lower resolution to have a vaguer background that wont dominate over what you overlay

[maxrjones]

It would be nice to have a way to query the auto-resolution determination without actually downloading anything. In additional to potential benefits for users, this would enable testing of the feature without slowing things down with more images-based tests and a larger cache. This sub-feature could be added as a new option or modifier for grdcut, since this is currently the recommend way for getting a grid from the remotes.

[PaulWessel]

As you know, we can use grdcut with -R -J now to extract a rectangular encompassing grid given an oblique area. Of course, until now we would do this with a known resolution. Seems the only thing I need to allow for is to let grdcut, if both -R -J are given, to allow the missing resolution. The rest will happen as it does for grdimage. As for an option or modifier to just report the region and resolution: That is a good idea which would also be helpful even for the old oblique -R -J case. Given available options, perhaps -D for dry-run is OK? (I dont like -I for inquiry (used in histogram) since -I should really be saved for grid increments if we can help it.)

[maxrjones]

Given available options, perhaps -D for dry-run is OK?

Sounds good

[PaulWessel]

So adding a -D works something like this:

gmt grdcut @earth_relief -Rd  -JG13.0550/47.8095/2000/0/0/0/60/60/14.9c  -Gt.grd -D
gmt [NOTICE]: Given -Rd -JG13.0550/47.8095/2000/0/0/0/60/60/14.9c, D_g = 36.5382 degrees and D_m = 8.29598 inches and given dpu = 300i so we replace @earth_relief by @earth_relief_01m_p
gmt [NOTICE]: Extracted grid region will be -5.21667/31.3333/36.9833/57.3333 for increment 01m

The first message will be a -Vi message once this is done. What do we want the second line to be? A -R string like coast -E sends to stdout? A notice like this to stderr? Since presumably both -R and -I are interesting here so we write -Rw/e/s/n -Iinc? Thoughts? If anyone is hoping to script things to recover the two settings then it is hard if it is a stderr message. Or, -D could take a filename where to write the info, either instead of the message or in addition to it.

[maxrjones]

So adding a -D works something like this:

gmt grdcut @earth_relief -Rd  -JG13.0550/47.8095/2000/0/0/0/60/60/14.9c  -Gt.grd -D
gmt [NOTICE]: Given -Rd -JG13.0550/47.8095/2000/0/0/0/60/60/14.9c, D_g = 36.5382 degrees and D_m = 8.29598 inches and given dpu = 300i so we replace @earth_relief by @earth_relief_01m_p
gmt [NOTICE]: Extracted grid region will be -5.21667/31.3333/36.9833/57.3333 for increment 01m

The first message will be a -Vi message once this is done. What do we want the second line to be? A -R string like coast -E sends to stdout? A notice like this to stderr? Since presumably both -R and -I are interesting here so we write -Rw/e/s/n -Iinc? Thoughts? If anyone is hoping to script things to recover the two settings then it is hard if it is a stderr message. Or, -D could take a filename where to write the info, either instead of the message or in addition to it.

I think the default should be stdout. There could be a +ttxtfile modifier for instead writing to file and perhaps +i for increment only.

Here's how I would imagine the output:

gmt grdcut @earth_relief -Rd  -JG13.0550/47.8095/2000/0/0/0/60/60/14.9c  -Gt.grd -D
-R-5.21667/31.3333/36.9833/57.3333 -I01m

gmt grdcut @earth_relief -Rd  -JG13.0550/47.8095/2000/0/0/0/60/60/14.9c  -Gt.grd -D+i
-I01m

[joa-quim]

D_m = 8.29598 inches

SI, please.

We should not need to set a dump output grid (-Dt.grd)

-D+n for externals (and others) like grdinfo -C with numerical [w e s n inc]

[PaulWessel]

D_m = 8.29598 inches

SI, please.

Yes, per PROJ_LENGTH_UNIT setting. This is just temporary (and because internal values are in inches).

We should not need to set a dump output grid (-Dt.grd)

No grid output via -D was considered. I mean did we want to write the information (text) to a file via D, or to stdout. Did you means something else?

-D+n for externals (and others) like grdinfo -C with numerical [w e s n inc]

That I like. We can emulate grdinfo -C on this without the baggage.

[joa-quim]

No grid output via -D was considered. I mean did we want to write the information (text) to a file via D, or to stdout. Did you means something else?

Yes, my bad, I meant the -Gt.grd in your example. I take it to mean nothing an only please the parser for the mandatory -G, but this can be avoided.

[PaulWessel]

My thinking was that -D means -G is not allowed (or could be ignored).

The only tricky thing is that the new stuff happens in gmt_init_module before parse so I am fishing the -D inside there. But then grdcut itself basically exits after parse. However, for some of the other modes of grdcut we may also wish to know the new region without getting a grid, so then I would need to know if I wrote stuff during gmt_init_module or not. We'll see. I think teh changed KEYS will be

#define THIS_MODULE_KEYS "<G{,FD(=,>D),GG}"

and then a special check in GMT_Encode_Options to turn off implicit output (-G) if -D is given. This is to enable external usage like

info = gmt ('grdcut', '@earth_relief -R... -J... -D');

to receive the dataset and not a grid.

[PaulWessel]

Any other concerns with this PR?

[joa-quim]

Still need to batptize -D. information is too vague. It's information about autoresolution, which too big to put in a word.

[maxrjones]

Federico's comment about earth_day still needs to be addressed: #5753 (comment)

[PaulWessel]

Federico's comment about earth_day still needs to be addressed: #5753 (comment)

OK, I am doing a similar zoom-in from 800 km to 10 km and because of #5763 I am not there yet. I will try to work on that one - we can let this one wait a bit longer. It is possible the algorithm could be a bit more generous when lower-resolution grids/images are selected to avoid this.

[maxrjones]

Apart from the earth_day comment, the only other issue that remains is adding a test that includes each remote dataset type using grdcut -D.

[PaulWessel]

Test checking the result of grdcut -D has been added, @meghanrjones.

[maxrjones]

Test checking the result of grdcut -D has been added, @meghanrjones.

Great, thanks!

It doesn't seem that merging in the code from #5776 fixed the earth_day resolution issue.

[PaulWessel]

Now that I have correct -R I can look at that to see if the algorithm needs improvement.