"Best" projection for rectangular world map?

Unfortunately still no joy:

julia> GMT.GMTdevdate
0001-01-01

(@v1.8) pkg> up
    Updating registry at `~/.julia/registries/General.toml`
  No Changes to `~/.julia/environments/v1.8/Project.toml`
  No Changes to `~/.julia/environments/v1.8/Manifest.toml`

julia> using GMT

julia> GMT.GMTdevdate
0001-01-01

Edit: Julia reports the following version:

julia> using Pkg

(@v1.8) pkg> st GMT
Status `~/.julia/environments/v1.8/Project.toml`
  [5752ebe1] GMT v0.44.8

This is new stuff in master. Do ] add GMT#master

Ah, thank you @Joaquim. Now I get

julia> GMT.GMTdevdate

2023-04-13

which isn’t the date you were showing in your post above.

Running your code:

G, cl = worldrectangular("@earth_relief_30m_g", pm=30, latlim=(80,-70), coast=:low);
imshow(G, shade=true, show=false, cmap=:geo)
plot!(cl, show=1)

Gives me a plot but not like yours:

Do I need to invoke the master branch somehow? I did the usual using GMT and also tried using GMT#master.

That’s fine. It is the date when GMT was built in BinaryBuilder. Mine is different because on Windows what is used is the installed GMT (the one I build) or our installer.

I’m having different errors from yours in WSL, which is also quite strange but it seem to indicate that I have non-committed stuff. Well I know I do, but didn’t think it was needed to reproduce this example. I’ll come back when I update with my local changes.

No, just do the usual using GMT

This is indeed very strange. I updated master but get the same as you in WSL. I wonder if I screw something with the instructions to build GMT in BinaryBuilder (Julia stuff) and those GMT binaries are not updated to what it needs for this. The shufflings with bits of the global grid are not even done and we don’t see a full filled rectangular grid.

The only way I see for now on unix is to see if these intended instructions to make use of system wide installed GMT do work.

# The file deps/deps.jl is created by compile from the deps/build.jl
# On Windows we use a system wide GMT if it is found from path or install it from a GMT installer. It is a MSVC binary.
# On Unix the default is to use the GMT_jll artifact. However this can be changed to use a system wide GMT installation.
# To swap to a system wide GMT installation, do (in REPL):
# 1- ENV["SYSTEMWIDE_GMT"] = 1;
# 2- import Pkg; Pkg.build("GMT")
# 3- restart Julia
#
# Note the above will work up until some other reason triggers a Julia recompile, where the JLL artifacts 
# will be used again. To make the ENV["SYSTEMWIDE_GMT"] = 1 solution permanent, declare a "SYSTEMWIDE_GMT"
# environment variable permanently in your .bashrc (or whatever).

But to keep teasing, this is what it was supposed to be possible by now. @pwessel how difficult would it be to let custom annotations have slanted ticks?

1 Like

I am surprised they are not slanted unless MAP_OBLIQUE_ANNOT has been set or changed (?)

They can’t slant because they are custom annotations and as far as they know, they are plotting on a cartesian axis.

Sure, if this is added as a Cartesian overlay then there are no slanting. I cannot see a simple method to overcome that.

@KristofKoch, the problem is more serious. The failure in unix is coming from the gdalwarp command not doing its job. Unfortunately, there are error messages issued by GDAL that I’ve not been able to capture in the Julia wrapper and maybe this is one of these cases. This makes even more needed to test if my recipe above makes it work.

It’s a GDAL bug. On Windows this works. It screams a lot but it works (it doesn’t for .nc)

gdalwarp -t_srs "+proj=vandg +over +pm=30" earth15.grd lixo.tiff

but not on linux. It screams equally but doesn’t work
… unless we repeat the command. Then no screaming and it works. That is, when overwriting it works but not when the file doesn’t exist.

What about inside annotations? Is it possible to have them oriented along meridians and parallels? Or does the same problem apply as with outside annotations?

Nope. Only what custom axes allow

To finish this (for now). I can’t make it work on Linux until GDAL fixes the bug (opened an issue). I thought that I could cheat it by running twice saving on disk and reloading but no way. Not even that worked.
On Win we can now do (function names are not definitive)

julia> G, cl = worldrectangular("@earth_relief_30m_g", pm=30, latlim=(80,-70), coast=:low);
julia> grid  = worldrectgrid(proj="+proj=vandg");
julia> grdimage(G, shade=true, plot=(data=cl,), cmap=:geo, B=:none)
julia> GMT.plotgrid!(G, grid);

(@Andreas your attempts at are perfectly doable with some changes in the functions that do this, and not only coastlines but grids too.)

2 Likes

For reference: the issue @Joaquim opened can be found over at GDAL under #7644.

As I unfortunately can’t experiment with your work – is it still possible to extend beyond the +/-180° longitude domain as in your earlier work above? To get a chart that spans 400° longitude at the Equator?

That example goes from [-540 540], so it should be a matter of calculating what are the projected coords at +/- 400 and trim the plot at those XX.

@KristofKoch, you saw Evan’s answer in the GDAL issue. Only way for now is to have on-the-edge GDAL+PROJ build and follow my recipe above to use a local GMT build.

I have now a minimal working code, though are still many issues that need future work. Namely, the default title offset is not enough to not overlap the custom annotations and those, easily clutter near the prime meridian.

G, cl = worldrectangular("@earth_relief_30m_p", proj="vandg", latlim=(-70,85), coast=true);
grid = worldrectgrid(G);
grdimage(G, shade=true, plot=(data=cl,), cmap=:geo, B=:none, title="Van der Grinten")
plotgrid!(G, grid, figname="vandg.png", show=true)

1 Like
G, cl = worldrectangular("@earth_relief_30m_p", proj="aitoff", coast=true, latlim=(-80,89));
grid = worldrectgrid(G, annot_x=[-180,-150,0,150,180]);
grdimage(G, shade=true, plot=(data=cl,), cmap=:geo, B=:none, title="Aitoff")
plotgrid!(G, grid, figname="aitoff.png", show=true)

1 Like

And there is something strange with Winkel Tripel. It takes much longer to compute and the file is > 3 times larger

G, cl = worldrectangular("@earth_relief_30m_p", proj="wintri", coast=true, latlim=(-80,89));
grid = worldrectgrid(G);
grdimage(G, shade=true, plot=(data=cl,), cmap=:geo, B=:none, title="Winkel Tripel")
plotgrid!(G, grid, figname="wintri.jpg", show=true)

2 Likes

That is very impressive @Joaquim - thank you! Now on to installing the bleeding edge GDAL+PROJ versions and then follow your recipe … that day job takes up way too much of my time …