Tobler’s hyperelliptical projection

[mapnerd2022]

I (perhaps mis-) understood the “but formulas were far more complicated” as being the reason.

Yes, that's also what I understood, but this passage is referring only to Craster's hyperbolic projection (Putnins P6), and I wanted to know whether Snyder also gave a reason for Craster rejecting his elliptical projection (Putnins P2), which is much less similar to the sinusoidal.

He unfortunately didn't. He only writes that it is graphically similar to the Eumorphic. I really don't know if Putnins rejected the elliptical for the same reason, meaning, that his elliptical had more complicated formulas than the basic Sinusoidal.

[Atarimaster]

I go by the formulas in Snyder’s Flattening the Earth, in which the only substantive difference from Tobler’s paper is the sign correction that you note. The Wikipedia article matches Snyder. Of course, I don’t really know what Tobler’s intended formulas are, given that the paper could have typos, but I get results that appear to match Tobler’s…

There are things that I find rather irritating.
Here’s Geocart’s preferred parameterization, red lines are rendered by Geocart, the black projection in the background is from Tobler’s paper:

tobler-fig-9.jpg (139.63 KiB) Viewed 7512 times

You see that there are differences.
However, other map projection applications that I’ve tried (G.Projector, Justin Kunimune’s MapDesignerRaster and the d3 scripts) render this parameterization exactly like Geocart. So I guess that this is the correct result according to the published formulae – but why is Tobler’s original figure different?

I checked the other are examples – in Fig. 6 and 8, there are similar differences.
Again, the d3 scripts render this parameterizations exactly like Geocart, so I don’t have to show it.

tobler-fig-6-and-8.jpg (267.22 KiB) Viewed 7512 times

In Fig. 5 and 7, the differences are so obvious that I don’t even have to put the two version on top of each other.
The d3 result from Fig.5 is again identical to Geocart’s, Fig. 7 is … astonishingly different from both Tobler’s and Geocart’s versions, so I’m showing it here (blue lines):

tobler-fig-5-and-7.jpg (110.45 KiB) Viewed 7512 times

And as if all that wasn't confusing enough, Toblers writes in the header of TABLE 3:
“Distortion on a General Hyperelliptical Projection with α = 0, k = 2.5, and γ = 1.183136, Similar to That Shown in Figure 9”

Similar? Why similar? Shouldn’t it be “as Shown in Figure 9”?

Are there any explanations for this?
And does somebody have an idea how it is possible that, as far as I can tell, the d3 scripts render the hyperelliptical projections identical to Geocart but then suddenly in Figure 7 it’s totally different? (You can check that I’ve used the correct parameters using this example, although you have to round γ to 1.6)

[daan]

In Fig. 5 and 7, the differences are so obvious that I don’t even have to put the two version on top of each other.
The d3 result from Fig.5 is again identical to Geocart’s, Fig. 7 is … astonishingly different from both Tobler’s and Geocart’s versions, so I’m showing it here (blue lines):

tobler-fig-5-and-7.jpg


And as if all that wasn't confusing enough, Toblers writes in the header of TABLE 3:
“Distortion on a General Hyperelliptical Projection with α = 0, k = 2.5, and γ = 1.183136, Similar to That Shown in Figure 9”

Similar? Why similar? Shouldn’t it be “as Shown in Figure 9”?

Are there any explanations for this?

Numerics are hard. Tobler’s projection in the forward direction is a minefield of numerical difficulties. This means that projecting vectors in the hyperelliptic is hard. Tobler’s renderings were done in 1973, when getting time on a computer was expensive, programming the computer was tedious, executing the program was slow, and supporting libraries might be primitive or lacking.

In order to compute the forward projection, you have to deploy two difficult numerical operations, one nested within the other: quadrature (evaluating the integral) and root-finding. Getting 9-digit accuracy (Geocart’s goal for worst case) out of that takes an unreasonably long time, even on modern computers.

I’m guessing Tobler was working with limited precision and that the errors compounded due to the nested numerical operations to the point of yielding these obvious differences: he might haven barely eked out two digits of accuracy by the end of all that. This might also explain his “similar to” comment: He didn’t have enough accuracy to claim “is as”.

The amount work I put into implementing this projection (twice now) I would characterize as “unreasonable”. It’s by far the most difficult pseudocylindric I have ever implemented, both mathematically and numerically due to the integral being of a form that is hard to evaluate accurately when approaching the point-as-pole cases as well as the problem of finding a suitable seed value for the root-finding in the same circumstances.

I will look into the D3/Geocart discrepancy when I get a few minutes.

— daan

[daan]

I will look into the D3/Geocart discrepancy when I get a few minutes.

I haven’t tried to decipher D3’s code here; it’s not easy to follow. What I will note is that the implementation at observablehq cannot be right: There is no mechanism in the hyperelliptic projection to yield the vertical stretch that it gives with increasing 𝛾. Given that both Geocart implementations agree with each other and with Tobler’s rendering (approximately), I don’t think I’m off base here. I put an immense amount of analytical work into this wretched projection.

A second observation is that the observablehq implementation allows parameterizations that are not possible using the real projection definition. For example, any 𝛾 lower than 1.27323954473516 (= 4/π), yielding the Mollweide, can’t work with k = 2.

One thing I will note is that Geocart as shipping has a horizontal compression and vertical stretch to accord with Tobler’s recommendation for his preferred projection. Geocart applies that regardless of the parameterization, which really isn’t appropriate. You can eliminate or otherwise adjust it as you see fit in the Stretch and Rotation dialog, of course, but doing it that way is cumbersome because the compression you would be compensating for is built into the projection and therefore not visible in the Stretch and Rotation dialog.

In my new implementation, I have added two parameters to better control the experience. The one concerning “no pole line” is to let you avoid having to figure out by trial and error the right 𝛾 to eliminate the pole line given the k that you chose. Without it, you are stuck playing chicken with “incompatible parameters” trying to get there unless you just work it out analytically given the complicated formula I gave several postings ago.

Screen Shot 2023-05-18 at 10.13.38 AM.png (21.11 KiB) Viewed 10307 times

Cheers,
— daan

[Atarimaster]

Thanks for your two very informative postings!

I’m guessing Tobler was working with limited precision and that the errors compounded due to the nested numerical operations to the point of yielding these obvious differences: he might haven barely eked out two digits of accuracy by the end of all that. This might also explain his “similar to” comment: He didn’t have enough accuracy to claim “is as”.

That sounds like a sensible explanation.

I haven’t tried to decipher D3’s code here; it’s not easy to follow. What I will note is that the implementation at observablehq cannot be right

I was afraid you might say something like this.
Damn, it is very handy to have some sliders that you can drag to try out various configurations (that’s what I have done quite a lot for the generalized Wagner).

In my new implementation, I have added two parameters to better control the experience. (…)

This sounds and looks promising!

Kind regards,
Tobias

[daan]

I’m guessing Tobler was working with limited precision and that the errors compounded due to the nested numerical operations to the point of yielding these obvious differences: he might haven barely eked out two digits of accuracy by the end of all that. This might also explain his “similar to” comment: He didn’t have enough accuracy to claim “is as”.
— daan

I think this hypothesis is further bolstered by the value he gives for his preferred parameterization: 𝛾 = 1.183136. I think it’s clear that he intended for there not to be a pole-line, but his parameterization gives a hardly discernible pole line. To eliminate it, the value would be 𝛾 = 1.183104547. Not having to integrate all the way to the projection’s parametric limit greatly reduces the computation time because it puts the integrand into a numerically stable region. In this case, he’s up to 5 digits accuracy in having computed the 𝛾 value, but that’s a single value that he might have felt he could afford to spend the computational resources on, whereas creating a map means repeating such a calculation thousands of times. And, depending on how much analytical work he did, he would not have even needed to compute the integral for 𝛾, as I show with the formula for computing 𝛾 given k and the constraint of no pole-line.

One other thing I note is that you can mess up the integral calculation, and therefore the y coordinate, without compromising the projection’s equivalence (equal-area-ness) because the computation for x, which is based on whatever y you got, automatically gives you an equal-area projection regardless. It took me awhile to realize I was not getting good results during development, since the projections I was getting were still computationally equal-area up to high accuracy. What happens is that you lose the highest latitudes; they get shoved “over” the pole into never-never land. The results in mild cases look a lot like the difference between Geocart’s rendering and Tobler’s, where the obvious problem is the change in the perimeter’s curvature in mid latitudes.

— daan

[Milo]

Here’s Geocart’s preferred parameterization, red lines are rendered by Geocart, the black projection in the background is from Tobler’s paper:
tobler-fig-9.jpg
You see that there are differences.

Care to check how my version stacks up?

tobler_fig5.png (16.04 KiB) Viewed 10290 times

tobler_fig8.png (28.33 KiB) Viewed 10290 times

tobler_fig9.png (30.97 KiB) Viewed 10290 times

Tobler wasn't clear about the affine coefficients he used for most of his figures, so I stuck to a 2:1 aspect ratio for everything. I can provide other sizes or aspect ratios if necessary.

(Total time taken for computing all three projections put together, including preprocessing and overzealous PNG compression, is about 0.66 seconds.)

In my new implementation, I have added two parameters to better control the experience. The one concerning “no pole line” is to let you avoid having to figure out by trial and error the right 𝛾 to eliminate the pole line given the k that you chose. Without it, you are stuck playing chicken with “incompatible parameters” trying to get there unless you just work it out analytically given the complicated formula I gave several postings ago.

I really think that even if you wanted cropped projections in the vein of the non-minimal-ɣ projections (which, again, Tobler didn't, since all three of his examples with α = 0 also had minimal ɣ, or as close to minimal as he could get on the computing hardware of the time), there are better ways to parametrize that choice than Tobler's ɣ parameter. But then, I guess it wouldn't be the "Tobler projection" if you didn't support the same parameters that Tobler did?

At least when α = 0, the following parameters work well:
ɣ = 1 / ∫₀^δ (1 − z^k)^1/k dz
β^k + δ^k = 1

When α ≠ 0, it gets a little more complicated, and I favor the following tweak:
ɣ = α_Milo/δ + (1 − α_Milo) / ∫₀^δ (1 − z^k)^1/k dz
α_Tobler = α_Milo/(ɣδ)
β^k + δ^k = 1
(Maybe I should have scrounged up another Greek letter for my variant parameter instead of stealing Tobler's, but calling things α and β is so convenient...)
Assuming that I did my math right, this will result in a pole line and equator having a length ratio of 1−(1−α)(1−β) = α + β − αβ, with a point pole when both are 0.

My own code doesn't implement Tobler's parameters, but it would be pretty simple to adjust it to support these new parameters, theoretically allowing it to render all of the same projections as Tobler's official formulae, even though actually converting a Tobler ɣ parameter into something I can use would be a serious hassle. Not that I really care to, since I don't really like pole lines, and my own "different exponents for x and y" innovation gives the same advantages as a pole line without actually having one.

One other thing I note is that you can mess up the integral calculation, and therefore the y coordinate, without compromising the projection’s equivalence (equal-area-ness) because the computation for x, which is based on whatever y you got, automatically gives you an equal-area projection regardless.

That can't possibly be right. If it were, then you would be able to rearrange the parallel spacing in the Mollweide projection arbitrarily, while retaining an elliptical overall contour, and still have an equal-area projection.

Of course, if your computation errors are small, then your deviation from having equal area is also going to be small.

What happens is that you lose the highest latitudes; they get shoved “over” the pole into never-never land.

That's what happens when you compute the wrong value for ɣ (which, as discussed above, needs to be computed with an integral if you want a point pole). I don't think it's that simple if you mess up the integrals in calculating the actual projection. (For example, you could overestimate some steps in the numerical integration and underestimate others, so that your total integral over the whole projection ends up being just about correct, but some of your partial integrals are wildly off.)

Also, a moderation note: Since the last five pages have been almost entirely about the Tobler projection rather than experimental projections, perhaps a thread split is in order?

[daan]

Also, a moderation note: Since the last five pages have been almost entirely about the Tobler projection rather than experimental projections, perhaps a thread split is in order?

Done, along with two more splits.

— daan

[daan]

In my new implementation, I have added two parameters to better control the experience. The one concerning “no pole line” is to let you avoid having to figure out by trial and error the right 𝛾 to eliminate the pole line given the k that you chose. Without it, you are stuck playing chicken with “incompatible parameters” trying to get there unless you just work it out analytically given the complicated formula I gave several postings ago.

I really think that even if you wanted cropped projections in the vein of the non-minimal-ɣ projections (which, again, Tobler didn't, since all three of his examples with α = 0 also had minimal ɣ, or as close to minimal as he could get on the computing hardware of the time), there are better ways to parametrize that choice than Tobler's ɣ parameter. But then, I guess it wouldn't be the "Tobler projection" if you didn't support the same parameters that Tobler did?

Right. My sole concern was to implement Tobler’s published hyperelliptic projection as accurately as feasible. I don’t care to develop new pseudocylindric projections; that’s a broadly plowed field, and I don’t really care for pseudocylindric projections in general. In fact, I don’t think I’ve ever intentionally¹ developed one. I view them as tools on the way to more interesting projections. I know that’s not a popular opinion, but, well, I’m not a popular guy.


— daan

¹We’ve all made mistakes programming…

[daan]

One other thing I note is that you can mess up the integral calculation, and therefore the y coordinate, without compromising the projection’s equivalence (equal-area-ness) because the computation for x, which is based on whatever y you got, automatically gives you an equal-area projection regardless.

That can't possibly be right. If it were, then you would be able to rearrange the parallel spacing in the Mollweide projection arbitrarily, while retaining an elliptical overall contour, and still have an equal-area projection.

Of course, if your computation errors are small, then your deviation from having equal area is also going to be small.

What happens is that you lose the highest latitudes; they get shoved “over” the pole into never-never land.

That's what happens when you compute the wrong value for ɣ (which, as discussed above, needs to be computed with an integral if you want a point pole). I don't think it's that simple if you mess up the integrals in calculating the actual projection. (For example, you could overestimate some steps in the numerical integration and underestimate others, so that your total integral over the whole projection ends up being just about correct, but some of your partial integrals are wildly off.)

I don’t think I miswrote. Let’s say that you mess up the function that you integrate in this fashion: Instead of ∫₀^𝑦(𝛾^𝑘 − 𝑧^𝑘)^1/𝑘 d𝑧, you apply some arbitrary function to the non-constant part of the integrand, like so: ∫₀^𝑦(𝛾^𝑘 − 𝜃(𝑧))^1/𝑘 d𝑧. In this case, d𝑦/d𝜑 = 𝛾 cos 𝜑 / (𝛾^𝑘 − 𝜃(𝑦(𝜑)))^1/𝑘. Now, since you accidentally bungled the integrand, you will think that 𝜃(𝑦(𝜑)) is 𝑦(𝜑), unaware that it’s mangled, and you will use it in your computation for 𝑥, like so: 𝑥 = 𝜆 (𝛾^𝑘 − 𝜃(𝑦(𝜑)))^1/𝑘 / 𝛾. Of course d𝑦/d𝜆 is still 0, so d𝑥/d𝜑 doesn’t matter, but d𝑥/d𝜆 has the weird term in the denominator while d𝑦/d𝜑 has it in the numerator, so they cancel out under multiplication. You still get d𝑦/d𝜑 ∙ d𝑥/d𝜆 − d𝑥/d𝜑 ∙ d𝑦/d𝜆 = cos 𝜑. Therefore the projection is equivalent regardless.

Notice that I just infinitely expanded the universe of pseudocylindrical projections you can create with the same basic structure. I don’t care about pseudocylindrical projections, though.


— daan