Unfortunately, at the moment it falls apart when using a central meridian that is not 0° longitude. See the following example with lon0 = 20 So not really ready for a PR.
red - @timhume special version
blue - GMT default
Things get mirrored and other strange things happen … not really sure what’s going on. Maybe Tim can shed some light. I tried to understand what’s happening but failed.
I never tried the Van der Grinten with a non-zero central longitude.
Regarding the Winkel Tripel, I used a formula from Wikipedia. I saw the way it was done in the GMT code was different. I trust the GMT code more. I’ll double check with Snyder’s map projection book.
Unfortunately I’ll have to wait a bit to look into this more. It’s just turned into Monday here in Australia, and you all know what that means Hopefully I’ll get to look at it a bit more after work.
That work has to go to a showcase
Good morning Tim, thank you for your ongoing effort. May you have an easy start into the new week. In the Old World the weekend comes to an end as well now …
My interest in the non-zero central longitude is mainly due to artistic reasons – the van der Grinten map looks a bit more balanced with the continents a bit more centered.
See how the Pacific Ocean on the left half looks pretty empty while the Asia-Pacific region on the right half looks cramped when looking at the left plot? On the right plot the central meridian is at E030° longitude balancing the continents a bit nicer for a more pleasing map. At least in my eyes …
A longitudinal domain > 360° would also be great to have … I like some overlap for easier orientation when I have to change to the other side of the map. For a zero central longitude imagine a map starting “in the West” at E160° and spanning all the way “to the East” to W160° like this Mercator plot:
The map is extended for 20° longitude on both sides (highlighted in grey and blue). I like how it gives better context in the area around 180° longitude eg New Zealand and Fiji. Unfortunately this currently only works with some cylindrical projections.
Example code:
van der Grinten comparison:
gmt begin vdgcomp png
gmt coast -Rd -JV10c -Bg -Dc -W0.8p,red -G255/206/206 -A0/0/1
gmt coast -Rd -JV30/10c -Bg -Dc -W0.8p,red -G255/206/206 -A0/0/1 -X11c
gmt end show
400° longitude Mercator map:
gmt set PS_LINE_JOIN round
gmt set MAP_FRAME_TYPE plain
gmt set MAP_FRAME_PEN 3p,black
gmt begin 400lon png
gmt coast -JM2.5c -R160/180/-70/82 -Dc -A0/0/1 -Bg20 -BWSN -W2p,black -G206/206/206
gmt coast -JM45c -R-180/180/-70/82 -Dc -A0/0/1 -Bg20 -BSN -X2.5c -W2p,red -G255/206/206
gmt coast -JM2.5c -R-180/-160/-70/82 -Dc -A0/0/1 -Bg20 -BSEN -X45c -W2p,blue -G206/206/255
gmt end show
The Mercator chain idea was originally posted by @Andreas above.
OK, here’s the Van der Grinten code fixed to deal with the Equator issue, and also to have an arbitrary central longitude. First the awk script:
BEGIN{
pi = atan2(0, -1);
d2r = pi/180;
}
function abs(x) {
return (x < 0) ? -1*x : x;
}
function asind(x) {
return atan2(x, sqrt(1-x**2))/d2r;
}
function sind(x) {
return sin(x*d2r);
}
function cosd(x) {
return cos(x*d2r);
}
function tand(x) {
return sind(x)/cosd(x);
}
function lonlat_to_xy(lond, latd) {
theta = asind(abs(2*latd/180));
if (latd == 0){
x = lond/180;
if (lond < 0) x*=-1;
y = 0;
} else if (lond == 0){
x = 0;
y = tand(theta/2);
} else {
A = 0.5*abs(180/lond - lond/180);
if (abs(lond) > 180) A*=-1
A2 = A**2;
s = sind(theta);
c = cosd(theta);
G = c/(s + c - 1);
P = G*(2/s - 1);
P2 = P**2;
Q = A2 + G;
P2A2 = P2 + A2;
GP2 = G - P2;
iP2A2 = 1/P2A2;
x = (A*GP2 + sqrt(A2*GP2**2 - P2A2*(G**2-P2)))*iP2A2;
y = (P*Q - A*sqrt((A2+1)*P2A2 - Q**2))*iP2A2
}
if (lond < 0) x *= -1;
if (latd < 0) y *= -1;
print(x, y);
}
/^>/{
print($0);
next;
}
{
lond = $1+0;
latd = $2+0;
lonlat_to_xy(lond-lon0, latd);
}
And here’s the shell script:
#!/usr/bin/env bash
gmt coast -M -Dhigh -A0/0/1 -Rd -Wthin > lld.txt
cat lld.txt > ll.txt
awk ' /^>/{ print($0); next; } { print($1-360, $2); }' < lld.txt >> ll.txt
awk ' /^>/{ print($0); next; } { print($1+360, $2); }' < lld.txt >> ll.txt
gawk -v lon0=20 -f vdg.awk < ll.txt > xy.txt
gmt begin vdgrect png
gmt plot xy.txt -JX20 -R-1/1/-0.8/0.8 -Wthinnest,black
gmt end
The central longitude is specified on the awk command line (-v lon0=20).
And here is the one-true Van der Grinten map projection:
Though to be fair, I think your central longitude of 20E is better
I’m going to take a break from the Winkel Tripel for the time being; just getting the Van der Grinten sorted out was hard enough. If I get time, I’m kind of curious to implement the “extended” Van der Grinten in the Julia interface. This would allow for the formulae to be applied to bigger datasets, such as elevation data.
Edit: Well, the awk could work on elevation data converted to (x, y, z) in text format - but that’s pushing things a little too far.
Sorry, can chime on this but only after ~18 h GMT
But why don’t you guys use only mapproject?
Because the issue is with filling in the “corners” of the map. If you look at the photograph in the first post, you’ll see in the top left and right corners Alaska is repeated twice, in order to fill in the corners of the map (otherwise the map would be circular in shape, not rectangular).
If you pass a coordinate for some place in Alaska to mapproject (or grdproject for that matter), then it will provide (x, y) coordinates for just one of the places on the rectangular map, not both.
At least that’s my understanding of mapproject.
Initially, I wasn’t even sure if it was theoretically possible to fill in the “corners” of the map - except the original photo proved it was. But in general, I’m not sure if what can be done with the Van der Grinten can be applied to other projections (it certainly can for cylindrical projections though).
Yes, I understood that (the corners issue) but my point is if you are programming the projection in awk you might as well use mapproject to do the same, no?
Sorry, still struggling with Linux to let me build the GMT.jl docs without fscrewing (which works perfectly fine on Windows but than are the GH-actions who …)
Oh, you mean extend mapproject to be able to fill in the corners?
Yeah, that would be one route. However, being able to extend the map boundaries outside the range (-180, 180] (or [0, 360) ) without simply wrapping around (e.g. -190 becomes 170) would only make sense for some map projections - I think. For example, I don’t think extending the range outside of (-180, 180] makes sense for an orthographic map projection (maybe it does, but I cannot visualise it right now).
So essentially what the grubby little awk script does is replicate mapproject functionality for one specific map projection, and adds the “corners” feature as well.
You’ll need GMT.jl master to do this
cl = coast(dump=:true, res=:crude, region=:global);
cl_right = cl .+ 360;
cl_left = cl .- 360;
cl_vdg = lonlat2xy(cl, t_srs="+proj=vandg");
cl_right_vdg = lonlat2xy(cl_right, t_srs="+proj=vandg +over");
cl_left_vdg = lonlat2xy(cl_left, t_srs="+proj=vandg +over");
tmp = cat(cl_left_vdg, cl_vdg);
cl3_vdg = cat(tmp, cl_right_vdg);
plot(cl3_vdg, show=1)
The miraculous option is that +over
(see) that prevents the 360 wrapping.
Note, that option can be used from mapproject
as well with -J+proj=vandg+over
but my test with the ll.txt
file from the script above showed some lines running from the vertical sides to lon = 0.
@pwessel is there any easy place were we could implement an equivalent to +over
?
The advantage of that is we could in principle let pscoast
accept -R outside the 360 range and have continents/oceans painted, which is not possible with coastline dumps.
Clip rectangles at will
I doubt it. GMT has lots of checks on longitude range <= 360.
As @KristofKoch found, you can do the > 360 for some periodic cylindrical projections but not those special ones he is interested in.
Gentlemen, thank you for your input and your incredible work. I haven’t forgotten about you, Unfortunately that pesky bill paying work had me occupied quite a bit for the last days.
Too bad that +over
seems not to be feasible. I’ll experiment a bit with the coastline dumps and also the OSM coastlines to see how I can get them to be closed polygons.
Gridded data will be another issue … converting large grids into x y z text files … maybe I have a way around that. Needs more thought.
You can do +over with mapproject and gdalwarp
gdalwarp
is even better than what I had in mind. Thank you, @Joaquim
Can you guys elaborate a bit more on this, please? For somebody who is a slow thinker and far less experienced?
Thanks
@mkononets - Basically I try to plot data outside the usual 360° longitude domain. All geographic projections in GMT are constrained to the 360° longitude, 180° latitude domain.
You can daisy chain some cylindrical projections together to escape that limitation, but for other projections like van der Grinten this doesn’t work. See above posts for a discussion of this.
@timhume dove deep into the math behind the van der Grinten projection and, contrary to my efforts, got it to work beautifully. Unfortunately this is limited to xyz data in a text file. If I want to project gridded data I would have to either to convert the grid to a very inefficient text file or find some other way to deal with it. That grid to text abomination would be best done in Julia and @Joaquim was curious and kind enough to implement Tim’s math into the Julia GMT wrapper.
He also had the idea with gdalwrap
as it comes with the +over
option to escape that 360° longitude domain. If that really solves my problems remains to be seen but it is a great starting place for further experimentation.
I hope this helps a bit to understand what kind of shenanigans we are into in this thread.
Thanks, this was the missing part on my side.
It does but there are some minor issues in GDAL. It will error if output is netCDF but works for GeoTIFF, though it screams a lot.
@Joaquim it might be just me struggling with Julia but there is no way to get a center meridian >0?