# "Best" projection for rectangular world map?

**URL:** <https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715>\
**Category:** Q&A\
**Created:** [February 27, 2023, 4:24pm UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715 "2023-02-27T16:24:31Z")\
**Posts on this page:** 20\
**Page:** 3

<div class="post-metadata">

**Author:** ![KristofKoch](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/kristofkoch/32/314_2.png) [@KristofKoch](https://forum.generic-mapping-tools.org/u/KristofKoch)\
**Post date:** [March 26, 2023, 7:38pm UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/41 "2023-03-26T19:38:26Z")

</div>

> [@Joaquim](#):
>
> I’m now waiting for a GMT-Julia PR.

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.

 ![vdglon20](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/2X/e/ef30677df9c7a44beacd779960e1d2d747de7d1c.jpeg)

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.

---

<div class="post-metadata">

**Author:** ![timhume](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/timhume/32/296_2.png) [@timhume](https://forum.generic-mapping-tools.org/u/timhume)\
**Post date:** [March 26, 2023, 8:03pm UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/42 "2023-03-26T20:03:40Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![PlanetGus](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/planetgus/32/303_2.png) [@PlanetGus](https://forum.generic-mapping-tools.org/u/PlanetGus)\
**Post date:** [March 26, 2023, 8:04pm UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/43 "2023-03-26T20:04:29Z")

</div>

That work has to go to a showcase

---

<div class="post-metadata">

**Author:** ![KristofKoch](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/kristofkoch/32/314_2.png) [@KristofKoch](https://forum.generic-mapping-tools.org/u/KristofKoch)\
**Post date:** [March 26, 2023, 9:57pm UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/44 "2023-03-26T21:57:25Z")

</div>

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 …

> [@timhume](#):
>
> I never tried the Van der Grinten with a non-zero central longitude.

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.

 ![vdgcomp](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/2X/f/f4fe043786a88346d1397b24194a2ecc91b64a77.png)

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:

 ![400lon](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/2X/d/dd7f1ad4880ee1bd9c9283ec3e2697ce76fe59e2.png)

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](https://flint.soest.hawaii.edu/t/best-projection-for-rectangular-world-map/3715/2).

---

<div class="post-metadata">

**Author:** ![timhume](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/timhume/32/296_2.png) [@timhume](https://forum.generic-mapping-tools.org/u/timhume)\
**Post date:** [March 27, 2023, 9:39am UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/45 "2023-03-27T09:39:58Z")

</div>

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:

```auto
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:

```auto
#!/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:

 ![vdgrect](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/2X/d/de95850a8530015af3221e9ff7202b7b37709326.png)

Though to be fair, I think your central longitude of 20E is better 🙂

---

<div class="post-metadata">

**Author:** ![timhume](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/timhume/32/296_2.png) [@timhume](https://forum.generic-mapping-tools.org/u/timhume)\
**Post date:** [March 27, 2023, 9:55am UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/46 "2023-03-27T09:55:13Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [March 27, 2023, 11:10am UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/47 "2023-03-27T11:10:50Z")

</div>

Sorry, can chime on this but only after ~18 h GMT

But why don’t you guys use only mapproject?

---

<div class="post-metadata">

**Author:** ![timhume](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/timhume/32/296_2.png) [@timhume](https://forum.generic-mapping-tools.org/u/timhume)\
**Post date:** [March 27, 2023, 8:22pm UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/48 "2023-03-27T20:22:21Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [March 27, 2023, 11:30pm UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/49 "2023-03-27T23:30:15Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![timhume](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/timhume/32/296_2.png) [@timhume](https://forum.generic-mapping-tools.org/u/timhume)\
**Post date:** [March 28, 2023, 2:39am UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/50 "2023-03-28T02:39:11Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [March 30, 2023, 12:00am UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/51 "2023-03-30T00:00:42Z")

</div>

You’ll need GMT.jl master to do this

```auto
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](https://proj.org/usage/projections.html#longitude-wrapping)) 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

 ![GMTjl_tmp](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/2X/0/04defdf5ef46c4bc9d847d0c2d50ef0bb564eb64.png)

---

<div class="post-metadata">

**Author:** ![pwessel](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/pwessel/32/30_2.png) [@pwessel](https://forum.generic-mapping-tools.org/u/pwessel)\
**Post date:** [March 31, 2023, 11:26am UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/52 "2023-03-31T11:26:58Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![KristofKoch](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/kristofkoch/32/314_2.png) [@KristofKoch](https://forum.generic-mapping-tools.org/u/KristofKoch)\
**Post date:** [April 8, 2023, 4:02am UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/53 "2023-04-08T04:02:29Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [April 8, 2023, 4:26am UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/54 "2023-04-08T04:26:15Z")

</div>

You can do +over with mapproject and gdalwarp

---

<div class="post-metadata">

**Author:** ![KristofKoch](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/kristofkoch/32/314_2.png) [@KristofKoch](https://forum.generic-mapping-tools.org/u/KristofKoch)\
**Post date:** [April 8, 2023, 12:44pm UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/55 "2023-04-08T12:44:45Z")

</div>

`gdalwarp` is even better than what I had in mind. Thank you, @Joaquim

---

<div class="post-metadata">

**Author:** ![mkononets](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/mkononets/32/4029_2.png) [@mkononets](https://forum.generic-mapping-tools.org/u/mkononets)\
**Post date:** [April 8, 2023, 1:20pm UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/56 "2023-04-08T13:20:38Z")

</div>

> [@Joaquim](#):
>
> You can do +over with mapproject and gdalwarp

Can you guys elaborate a bit more on this, please? For somebody who is a slow thinker and far less experienced?

Thanks

---

<div class="post-metadata">

**Author:** ![KristofKoch](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/kristofkoch/32/314_2.png) [@KristofKoch](https://forum.generic-mapping-tools.org/u/KristofKoch)\
**Post date:** [April 8, 2023, 2:09pm UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/57 "2023-04-08T14:09:10Z")

</div>

@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.

---

<div class="post-metadata">

**Author:** ![mkononets](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/mkononets/32/4029_2.png) [@mkononets](https://forum.generic-mapping-tools.org/u/mkononets)\
**Post date:** [April 8, 2023, 2:30pm UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/58 "2023-04-08T14:30:52Z")

</div>

> [@KristofKoch](#):
>
> He also had the idea with `gdalwrap` as it comes with the `+over` option to escape that 360° longitude domain.

Thanks, this was the missing part on my side.

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [April 8, 2023, 2:37pm UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/59 "2023-04-08T14:37:06Z")

</div>

> [@KristofKoch](#):
>
> If that really solves my problems remains to be seen

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.

---

<div class="post-metadata">

**Author:** ![KristofKoch](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/kristofkoch/32/314_2.png) [@KristofKoch](https://forum.generic-mapping-tools.org/u/KristofKoch)\
**Post date:** [April 18, 2023, 4:13pm UTC](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715/60 "2023-04-18T16:13:55Z")

</div>

@Joaquim it might be just me struggling with Julia but there is no way to get a center meridian \>0?

[Previous page](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715.md?page=2)

[Next page](https://forum.generic-mapping-tools.org/t/best-projection-for-rectangular-world-map/3715.md?page=4)