A dihedral projection

[Milo]

Here's a new hemisphere-in-a-square projection. I consider it a special case of a polyhedral projection, based on a rectangular dihedron. The parameters are the four corners of the rectangle, that form two pairs of antipodes (so there are only four scalar degrees of freedom, as choosing two corners will also set the other two). The square form, which is usually the most interesting, is obtained when all four edges of the rectangle have equal length (setting this constraint leaves three degrees of freedom: one corner can be set arbitrarily, and then you can choose direction, but not distance, to the next).

Note that here "rectangle" refers to a spherical rectangle, specifically a degenerate rectangle that lies entirely on a great circle, with interior angles of 180°. It is rectangular in the sense that opposite edges must have the same lengths, but adjacent edges may have different lengths. Regardless, its projection will still take the form of a Euclidean square.

Here's a picture of the projection in square Guyou aspect (explained below) with the hemispheres separated at 20°W and 160°E:

yay.png (223.83 KiB) Viewed 14795 times

Contrast the actual Guyou projection, which is conformal. They look very similar, to the point that you'd have to view them side-by-side to tell the difference, but this is a compromise projection, with (slightly) worse angles but (slightly) better areas than the conformal Guyou projection.

To my knowledge, this is the first-ever non-conformal polyhedral projection to maintain smoothness everywhere except the vertices.

Okay, so if it's not a conformal projection, then what's so great about it?

Other than being smooth, of course.

Well, it meets the theoretical optimum resolution-efficiency for any rectangular dihedral projection (which, as I explained there, is valuable for ensuring that no information is lost, as well as making good compromise projections in general). In fact, it has the same resolution-efficiency as the plate carree projection, while having less extreme singularities. (Indeed, the plate carree can be obtained as a special case of this one, in the limit as two corners approach the same point.) That means I would seriously recommend considering this projection for any application where you'd consider using the plate carree projection!

It's definitely my favorite so far of any of the projections I've invented myself. (Admittedly, there aren't many of those that weren't jokes.)

Although it is provably optimal in its category, it is not necessarily the only optimal projection. There could be others tied with it, that might be considered better in some more nebulous aspect. (In fact, I already have developed another resolution-optimal projection some time ago, although it is patently worse since it isn't smooth.) However, this projection certainly looks pretty good to me.

Also, it doesn't require elliptic functions to calculate. That's one major advantage over the Guyou projection.

In fact, although it took quite some time to stumble across the idea and work out the details, once I was done the final projection turned out to be surprisingly simple.

The x coordinate of the projection of a point is given by one of:

• The semiminor axis of the ellipse containing the point whose focus points are the top-left and bottom-left corners.
• The semiminor axis of the ellipse containing the point whose focus points are the top-right and bottom-right corners.
• The reciprocal eccentricity of the hyperbola containing the point whose focus points are the top-left and top-right corners.
• The reciprocal eccentricity of the hyperbola containing the point whose focus points are the bottom-left and bottom-right corners.

These aren't actually equal, but they are all congruent transformations of each other (adding or subtracting a constant and/or flipping the sign), so they result in essentially the same projection.

The y coordinate of the projection of a point is given by one of:

• The semiminor axis of the ellipse containing the point whose focus points are the top-left and top-right corners.
• The semiminor axis of the ellipse containing the point whose focus points are the bottom-left and bottom-right corners.
• The reciprocal eccentricity of the hyperbola containing the point whose focus points are the top-left and bottom-left corners.
• The reciprocal eccentricity of the hyperbola containing the point whose focus points are the top-right and bottom-right corners.

Again, these are all congruent transformations of each other. Which is more convenient depends in part on whether you're treating the positive y direction as being up (as conventional in mathematics), or treating the positive y direction as being down (as conventional in computer graphics).

It follows trivially that horizontal and vertical lines in the projection will correspond to ellipses on the sphere.

The semiminor axis of an ellipse is given by:
semiminor = acos(cos(sum_of_distances_to_foci / 2) / cos(distance_between_foci / 2))
The reciprocal eccentricity of a hyperbola (which is also an ellipse, but shh, don't tell anyone) is given by:
1/eccentricity = asin(sin(difference_between_distances_to_foci / 2) / sin(distance_between_foci / 2))
...Or at least, that's what I think this value should be called, based on trying to find the closest planar analogue. I initially just called it a semiminor axis too, but that appears to not be correct. Regardless, whatever the above value is called is what we're interested in. Up to congruence (in the choice of applying asin versus acos), it's the equivalent of the semiminor axis of the ellipse when one focus is flipped to its antipode.

The definitions in terms of acos are more convenient if you want the Cartesian coordinate (0,0) to lie at a corner of the square (as you typically would when generating computer graphics), while the definitions in terms of asin are more convenient if you want the Cartesian coordinate (0,0) to lie at the center of the square.

Thus, using the computer graphics convention, we get:
x(P) = acos(cos((d(P,TL) + d(P,BL)) / 2) / cos(d(TL,BL) / 2)) = acos(sin((d(P,TR) - d(P,TL)) / 2) / sin(d(TL,TR) / 2))
y(P) = acos(cos((d(P,TL) + d(P,TR)) / 2) / cos(d(TL,TR) / 2)) = acos(sin((d(P,BL) - d(P,TL)) / 2) / sin(d(TL,BL) / 2))
...where d is the usual great-circle distance. It's convenient to use the first form for x and the second form for y, or vice versa, since this means only two distances need to be computed.

The inverse form is trickier. Working backwards from these equations, we can find:
d(P,TL) = acos(cos(y) * cos(d(TL,TR) / 2)) - asin(cos(x) * sin(d(TL,TR) / 2))
d(P,TR) = acos(cos(y) * cos(d(TL,TR) / 2)) + asin(cos(x) * sin(d(TL,TR) / 2))
One can then use the following formulae to convert the two-point equidistant projection to an azimuthal equidistant projection:
distance = acos((cos(d(P,TL)) + cos(d(P,TR))) / (2 * cos(d(TL,TR) / 2)))
azimuth = acos((cos(d(P,TL)) - cos(d(P,TR))) / (2 * sin(d(TL,TR) / 2)) / sin(distance))
In fact, it's possible to cut down on the trigonometric functions a little:
distance = acos(cos(y) * √(1 - (cos(x) * sin(d(TL,TR) / 2))²))
azimuth = acos(cos(x) * √(1 - (cos(y) * cos(d(TL,TR) / 2))²) / sin(distance))
This gives azimuthal coordinates centered on the midpoint of the TL-TR edge, not on the midpoint of the square. It can then be converted to standard longitude-latitude coordinates using the usual oblique detransformation:
longitude = atan2(sin(distance)*sin(azimuth-top_center_bearing), cos(top_center_latitude)*cos(distance) + sin(top_center_latitude)*sin(distance)*cos(azimuth-top_center_bearing)) + top_center_longitude
latitude = asin(sin(top_center_latitude)*cos(distance) - cos(top_center_latitude)*sin(distance)*cos(azimuth-top_center_bearing))
Phew! Possibly this can be crunched down further, but this is the version which I use and it works.

The Tissot distortion parameters are, unfortunately, difficult to calculate. I figured them out, see below! We can make the following observations:

• The Tissot semiminor axis is, at all points, greater than or equal to one. This is, in fact, the property that makes the projection to be resolution-efficiency-optimal (in combination with its total area). However, it is not constant one (the way the plate carree or azimuthal equidistant projections are), so it is not a constant-resolution projection. I'm not sure if trying to find a constant-resolution projection would actually improve it, since the decrease in semiminor axes across the projection would necessarily come with an increase in semimajor axes, probably resulting in an increase in angle distortion, but it's something to look out for, anyway.
• Both semiminor and semimajor axis, and consequently scale, are unbounded and approach infinity near the vertices. However, they remain much lower across the majority of the projection, and of course the projection as a whole still has a finite (and quite compact) area.
• Angle distortion, on the other hand, is sharply bounded: in the square case, the semimajor axis is never more than √2 times the semiminor axis. (In the rectangular case, the ratio is 1/sin(half length of shorter edge). Notice how this tends to infinity in the limit.) Therefore, although not conformal, angles are preserved quite well overall. (Since I mentioned that resolution-optimal projections aren't unique, I have to wonder if this projection has minimax angle distortion within its category, which would give a solid motivation for using it.) In the square case, the projection is conformal on the diagonals (but not constant scale, of course), while the worst angle distortion lies at the midpoints of the edges, away from the vertices (which have the worst area distortion). The rectangular case is more coplicated, as angle distortion will drop to normal at the longer edge, but increase to infinity at the shorter edge.

Oh yeah, I promised to explain the aspects? I feel that the usual classification of projections into "normal", "transverse", and "oblique" aspects isn't quite applicable to polyhedral projections. Instead, for regular polyhedra, there are three types of non-oblique aspects, depending one chooses to place a face center, an edge center, or a vertex at the pole. (In the case of a non-centrally-symmetric polyhedron, such as a tetrahedron or an odd dihedron, one also has to specify which pole. The square dihedron used here is centrally-symmetric, so one gets the same type of thing at both poles. For more complicated polyhedra like the Catalan solids, one may also need to specify which type of vertex is being centered.) In the special case of the square or rectangular dihedron, I also refer to these three aspects as the Peirce, Guyou, and Adams aspects, in honor of the first people to apply them in the conformal case (but they're just as relevant to non-conformal projections).

• In the face-polar, or Peirce, aspect, all four corners lie in the equator. In the square case, this allow one degree of freedom (a longitude). In the rectangular case, this allows two degrees of freedom (two longitudes). The square case is obtained when the longitudes are 90° apart. In the limit as the longitudes approach either equality or antipodality, one obtains a transverse plate carree projection.
• In the edge-polar, or Guyou, aspect, two corners lie at equal distances on opposite sides of the north pole, and the other two likewise for the south pole. In the square case, this allows one degree of freedom (a longitude). In the rectangular case, this allows two degrees of freedom (a longitude and a latitude). The square case is obtained when the latitude equals 45°. In the limit case as the latitude approaches 90°, one obtains a normal plate carree projection. In the limit case as the latitude approaches 0°, one obtains a transverse plate carree projection.
• In the vertex-polar, or Adams, aspect, two corners lie at the poles, and the other two lie on the equator. In the square case, this allows one degree of freedom (a longitude). There is no natural way to extend this to the rectangular case. (One could place a rectangles with two vertices at the poles and the other two oblique, but this lacks the symmetry of rectangular projections in Peirce or Guyou aspect. If one did try it, the limit case would be a normal plate carree projection.)
• All aspects other than those three are considered oblique. This allows three degrees of freedom in the square case and four in the rectangular case, so always two more degrees of freedom than any of the above "normal" aspects.

My implementation supports all of these, but I treat the Guyou aspect as the "default", motivated by its particularly natural rectangular forms that produce the normal plate carree projection in the limit case. That's just here, though. For other polyhedra I would recommend other aspects. (The Guyou aspect also simplifies the inverse form formulae I gave above, as it allows the skipping of a nasty deobliquing transformation, although in the cases of Peirce and Adams I'd still say it's just the computation that's oblique, rather than the projection itself.)

Although, as usual, I haven't implemented graticules into my projector, here is a rendition of the projection's graticules made in gnuplot:

graticule.png (56.81 KiB) Viewed 14795 times

The code used for this is:

Code: Select all

set parametric
set size ratio -1
unset xtics
unset ytics
unset xzeroaxis
unset yzeroaxis
unset key
set multiplot layout 1,2
x(longitude, latitude) = asin(sin((+ acos((sin(latitude) - cos(latitude)*sin(longitude)) / sqrt(2)) - acos((+ sin(latitude) + cos(latitude)*sin(longitude)) / sqrt(2))) / 2) * sqrt(2))
y(longitude, latitude) = asin(sin((+ acos((sin(latitude) - cos(latitude)*sin(longitude)) / sqrt(2)) - acos((- sin(latitude) - cos(latitude)*sin(longitude)) / sqrt(2))) / 2) * sqrt(2))
set title 'Guyou'
plot [-pi/2:+pi/2] [-pi/2:+pi/2] [-pi/2:+pi/2] for[i=-6:+6] x(t, pi*i/12), y(t, pi*i/12) lc 1, for[i=-5:+5] x(pi*i/12, t), y(pi*i/12, t) lc 1
x(longitude, latitude) = asin(-cos((acos(cos(latitude)*cos(longitude)) + acos(cos(latitude)*sin(longitude))) / 2) * sqrt(2))
y(longitude, latitude) = asin( sin((acos(cos(latitude)*cos(longitude)) - acos(cos(latitude)*sin(longitude))) / 2) * sqrt(2))
set title 'Peirce'
plot [-pi/2:+pi/2] [-pi/2:+pi/2] [-pi/2:+pi/2] for[i=1:6] x(t*2, pi*i/12), y(t*2, pi*i/12) lc 1, for[i=0:23] x(pi*i/12, t/2+pi/4), y(pi*i/12, t/2+pi/4) lc 1

Note how all these curves are smooth, and how the parallels in the Guyou aspect, and meridians in the Peirce aspect, are perpendicular to the edges of the square.

Concentric circles around the center of the square, such as parallels in the Peirce aspect, are mapped as curves satisfying cos(x) * cos(y) = cos(r). I discovered this formula quite some time ago, so I'm happy to finally find a practical use for it!

Finally, here's my implementation, with Linux binary and C source code:

dihedral_rect.zip

(13.49 KiB) Downloaded 889 times

What's left to do:

• Figure out how to extend the principle used here to other polyhedra besides the rectangular dihedron.
• Figure out analogous equal-area projections.
• Convince people to actually use this projection!

[Milo]

Here's some more aspects (all of the square form).

This is the Adams aspect, with the hemispheres separated at 20°W and 160°E (same as the Guyou above):

resolution-preserving_adams.png (243.24 KiB) Viewed 14820 times

Note the extreme distortion of Antarctica. Also note that the Adams aspect is often rotated so the poles are at the top and bottom, but there is no convenient way to do this that shows both hemispheres at once.

This is the Peirce aspect, with corners at 20°W, 70°E, 160°E, and 110°W on the equator:

resolution-preserving_peirce.png (213.36 KiB) Viewed 14820 times

Here is the same aspect, spliced and rotated to the quincuncial arrangement favored by Peirce (in an image editing program, since my projector does not currently support this):

resolution-preserving_peirce_quincuncial.png (213.98 KiB) Viewed 14820 times

[Milo]

I had the realization that some of the calculations I made during my failed attempt to derive a confocal equal-area projection could be repurposed to help figure out the Tissot indicatrices here. The trick was to ignore longitude and latitude and work entirely in elliptic coordinates.

Using this, I have obtained the following formulae, where ψ is half the distance between the top-left and top-right (or bottom-left and bottom-right) corners:
x scale = √((1 - cos(x)² sin(ψ)²) / (1 - cos(x)² sin(ψ)² - cos(y)² cos(ψ)²))
y scale = √((1 - cos(y)² cos(ψ)²) / (1 - cos(x)² sin(ψ)² - cos(y)² cos(ψ)²))
Whichever of these is larger is the Tissot semimajor axis and whichever is smaller is the Tissot semiminor axis. The product is areal scale (flation), and the ratio determines angle distortion.
(This uses the useful-in-computer-graphics convention that (x,y)=(0,0) is the upper-left corner. If you want (x,y)=(0,0) to be the center, which may be more convenient mathematically, just swap cos(x) and cos(y) with sin(x) and sin(y).)

Plotting these against my earlier numerically-computed solutions in Maxima shows a perfect match, up to 5E-12 or so floating point rounding errors, so I'm confident I got it right.

In the square case, ψ = π/4, and this simplifies to:
x scale = √((1 + sin(x)²) / (sin(x)² + sin(y)²))
y scale = √((1 + sin(y)²) / (sin(x)² + sin(y)²))

It is obvious from these formulae that x scale and y scale can both never be smaller than 1 (as the numerator is smaller than the denominator, being the same value minus another positive term), therefore proving resolution-optimality.

Furthermore, we find that the projection is conformal whenever cos(x)² sin(ψ)² = cos(y)² cos(ψ)², or in the square case, whenever cos(x)² = cos(y)², i.e. on the diagonals of the square!

Here is a gnuplot rendition of the Tissot indicatrices:

tissot.png (32.02 KiB) Viewed 14815 times

Code: Select all

set size ratio -1
unset key
unset xtics
unset ytics
set samples 13
set isosamples 13,13
x_scale(x,y) = sqrt((1 - (sin(x)*sin(psi))**2) / (1 - (sin(x)*sin(psi))**2 - (sin(y)*cos(psi))**2))
y_scale(x,y) = sqrt((1 - (sin(y)*cos(psi))**2) / (1 - (sin(x)*sin(psi))**2 - (sin(y)*cos(psi))**2))
psi=pi/4; plot [-pi/2:+pi/2] [-pi/2:+pi/2] sample [u=-pi/2:+pi/2] [v=-pi/2:+pi/2] '++' using ($1):($2):(abs($1)==pi/2&&abs($2)==pi/2?0: x_scale($1,$2)/15):(abs($1)==pi/2&&abs($2)==pi/2?0: y_scale($1,$2)/15) with ellipses fill solid  # I don't know why gnuplot throws up an error message if I name the variable anything other than u.

[Atarimaster]

Hello Milo,

since I’ve always liked the Guyou projection, it’s not surprising that I also like your compromise counterpart.
Good work!

[Milo]

Thanks!

Here's an example of a truly oblique projection. You may have noticed that, in the square form, the projection is interrupted over a total arc length of 270° (for the left hemisphere: 90° on the bottom, 90° on the left, and 90° on the top, leaving only the 90° on the right where it's attached to the other hemisphere). According to Wikipedia, the longest uninterrupted geodesic lying entirely in the ocean runs from 25°25'N 66°25'E to 59°38'N 163°24'E (viewable here), which has a length of 32042 km... in excess of 270° (which is approximately 30000 km). So by setting the interruption inside this geodesic, we can obtain a map that is nowhere interrupted through land! First, I used the 90° arc lying exactly in the middle of this geodesic, giving a top-right coordinate of 178.867°E 55.619°N and a bottom-right coordinate of 61.345°E 17.546°N (view the interruption here):

no-aleutian.png (230.32 KiB) Viewed 14804 times

Unfortunately, this cuts through the Aleutian islands! Since they're an island arc rather than a contiguous landmass, the computation of the longest geodesic is allowed to pass through them, but the islands are close enough together it does seem desirable to keep them grouped. So to fix this, I adjusted the location of the interruption by 7°, obtaining a top-right coordinate of 171.608°W 51.482°N and a bottom-right coordinate of 65.163°E 23.566°N (view the interruption here):

aleutian.png (233.93 KiB) Viewed 14804 times

This interruption is cutting it very close, but it works! Note that the scale exaggeration at the vertices makes it seem like the interruption is cutting it less close than it really is. Nonetheless, all landmasses look pretty good, with the most noticeable bending being around Miani Hor in southwest Asia, but even that isn't enough to make the map look blatantly wrong. Shame that Madagascar's location relative to Africa is obscured, though.

(Note: all coordinates calculated by me use the spherical approximation, so are probably not totally accurate on the ellipsoid.)

Finally, to show more than just the square form... here's an animation of a smooth transition from the normal to transverse plate carree projection and back!

animation-24f.png (838.21 KiB) Viewed 14804 times

Note that the hemispheres stay the same (the top-left and bottom-right squares are the hemisphere between the 200°W and 20°W meridians, while the top-right and bottom-left squares are the hemisphere between 20°W and 160°E meridians), only the projection within each hemisphere/square shifts.

...And I think that's enough images. Gotta give Daan's poor server a rest

[mapnerd2022]

I really like your work here in the mapthematics forums! Maybe I'll present a projection or two here as well!

[Milo]

I really like your work here in the mapthematics forums! Maybe I'll present a projection or two here as well!

Looking forward to seeing what you can come up with! I'll warn you, though, that people have been developing map projections for thousands of years, so coming up with a genuinely good one that no-one has discovered yet is pretty hard!

On another note, I discovered a minor bug in my implementation where attempting to project a Guyou aspect with negative latitude would give a projection rotated 180° from how it should be. (It's technically still the same projection, just wrong side up!) I uploaded a new version that fixes this.

[Atarimaster]

One question:

This is the Adams aspect

I like that one, too … but why is it called “Adams aspect”?
I can’t remember to have seen this aspect, attributed to Adams… or anywhere else, for that matter. Although I probably did see it somewhere at some point – after all, I’ve to the Gringorten projection in a very similar aspect on my website. I also can’t remember how I came up with this one, but I doubt that it was my own idea, obtained by trial-and-error…

[Milo]

why is it called “Adams aspect”?

The conformal projection in that aspect is named the "Adams hemisphere-in-a-square projection" on both Wikipedia and Mapthematics... and for that matter, on your own site. Getting a little forgetful?

All three sources credit it to a publication by Oscar Sherman Adams 1925. That makes it the last of the three aspects to be presented, as Émile Guyou published in 1887 and Charles Sanders Peirce published in 1879.

Incidentally, I realized that the "quincuncial" arrangement, though typically used with the Peirce aspect, is actually even more appropriate for the Adams aspect, as it naturally involves rotating the hemisphere to a diamond orientation.

The use of these names to represent the analogous aspects on non-conformal projections is my own invention and not (yet) seen anywhere else.

[Atarimaster]

The conformal projection in that aspect is named the "Adams hemisphere-in-a-square projection" on both Wikipedia and Mapthematics... and for that matter, on your own site. Getting a little forgetful?

Not that forgetful (since I’ve recently rendered another projection based on Adams hemisphere in a square, which will be soon added to my website); I merely failed to see that the projection actually IS in equatorial aspect. Drawing graticule lines in my head still is a thing that I’m not good at.

Thanks for clearing that up!