Spilhaus projection

[AsgerPetersen]

Also see A story about the "Spilhaus projection" – a map projection that went viral in fall 2018 and will be supported in the next release of ArcGIS

The article finds that the spilhaus projection is a Adams World in a Square II projection with these params:

    Geographic Coordinate System: GCS WGS 1984
    Projection: Adams Square II
    False Easting: 0 m
    False Northing: 0 m
    Scale Factor: 1
    Azimuth: 40.17823482°
    Longitude of Center: 66.94970198°E
    Latitude of Center: 49.56371678°S
    XY Plane Rotation: 45°

It would be fun to have this projection in proj. This document suggests that the "Adams World in a Square II" projection was implemented in the libproject distribution of proj. So maybe it is even possible to find working code somewhere.

[kbevers]

Code is here: https://github.com/jeffbaumes/jeffbaumes-vtk/blob/master/Utilities/vtklibproj4/proj_guyou.c

Have fun translating it to modern day PROJ - Does seem to be five for one kind of deal, so there's that :-)

[AsgerPetersen]

I researched it a bit and wanted to document my findings in case someone would need it later. It seems to be a hot projection at the moment. At least ESRI is making a lot of noise about it at the moment because the new ArcGIS Pro supports it. See for instance this tweet by John Nelson. It even got a dedicated article in forbes a week ago.

[rouault]

Porting Adams Square II from libproj to PROJ wouldn't be hard. But there's more work to do since according to the ESRI story, Spilhaus adds "The Azimuth parameter is a direction from the North towards the top "pole" at the center on the conformal sphere"
There would be also other additional work because libproj Adams Square II is for spherical only:
"Adams only presented the forward equations for spherical Earth models. Most of today’s geospatial data is defined based on ellipsoidal models, such as WGS 1984 or GRS 1980. That data needs to be converted back from the projected surface to geographic coordinates and we needed to derive the forward and inverse equations for ellipsoidal Earth models. We achieved that by converting geodetic coordinates to a conformal sphere, conformally shrinking the model to a hemisphere, and resolving a complex elliptic integral of the first kind. "
So some cartographic expertise required here...

@beuan I see you mentionned as on of the author of the ESRI story. Any plan to contribute it to PROJ ?

[kbevers]

So some cartographic expertise required here...

I am fairly certain that expertise is reachable on our mailing list. You may not get lines of code but I am sure that someone can provide some insigt into as to how to solve the problem.

[stale]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.

[rouault]

Note: if support for Spilhaus is implemented, an action item would be to make sure that ESRI:54099 "WGS_1984_Spilhaus_Ocean_Map_in_Square" which was "blacklisted" in the mapping of ESRI WKT Adams_Square_II in #2157 is enabled

[stale]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.

[philippemiron]

I'm not familiar enough with PROJ, but I found this implementation. Is this something that could work if translated to proj syntax?

Adams World in a Square I + II

Porting Adams Square II from libproj to PROJ wouldn't be hard. But there's more work to do since according to the ESRI story, Spilhaus adds "The Azimuth parameter is a direction from the North towards the top "pole" at the center on the conformal sphere"
There would be also other additional work because libproj Adams Square II is for spherical only:
"Adams only presented the forward equations for spherical Earth models. Most of today’s geospatial data is defined based on ellipsoidal models, such as WGS 1984 or GRS 1980. That data needs to be converted back from the projected surface to geographic coordinates and we needed to derive the forward and inverse equations for ellipsoidal Earth models. We achieved that by converting geodetic coordinates to a conformal sphere, conformally shrinking the model to a hemisphere, and resolving a complex elliptic integral of the first kind. "
So some cartographic expertise required here...

@beuan I see you mentionned as on of the author of the ESRI story. Any plan to contribute it to PROJ ?

[stale]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.

[AJueling]

@rouault, this question may reveal my ignorance, but would it not be good enough to have the projection work only from spherical data for the majority of use cases which would be oceanographic and global in scope? I think this projection would be of great interest to oceanographers for predominantly illustrative (and perhaps less scientific) purposes.

[rouault]

would it not be good enough to have the projection work only from spherical data for the majority of use cases which would be oceanographic and global in scope?

A spherical version would be a good start, but someone would have to find the maths needed to implement the Azimuth and Latitude of Center parameters needed for Spilhaus, and not present in the PROJ implementation of Adams World in a Square II

[kerkovits]

would it not be good enough to have the projection work only from spherical data for the majority of use cases which would be oceanographic and global in scope?

A spherical version would be a good start, but someone would have to find the maths needed to implement the Azimuth and Latitude of Center parameters needed for Spilhaus, and not present in the PROJ implementation of Adams World in a Square II

This should be extremely easy. The map projection is conformal. If you want to map ellipsoidal data just use conformal latitude (formula 3-1 of https://pubs.usgs.gov/pp/1395/report.pdf) instead of geographic latitude and just use spherical formulae. This will be 100% conformal and works with every conformal map projection.

The "center of projection" and "azimuth" is just a graticule rotation before applying the map projection. (formula 5-7 and 5-8b of https://pubs.usgs.gov/pp/1395/report.pdf) I could reproduce the projection with alpha=30°, lambda_0=115°. Unfortunately the exact value of beta is unknown (it should be calculated using spherical trigonometry) but I could achieve an almost perfect result with beta=-28.8° (determined by trial and error).

So this should be possible.

[philippemiron]

would it not be good enough to have the projection work only from spherical data for the majority of use cases which would be oceanographic and global in scope?

A spherical version would be a good start, but someone would have to find the maths needed to implement the Azimuth and Latitude of Center parameters needed for Spilhaus, and not present in the PROJ implementation of Adams World in a Square II

This should be extremely easy. The map projection is conformal. If you want to map ellipsoidal data just use conformal latitude (formula 3-1 of https://pubs.usgs.gov/pp/1395/report.pdf) instead of geographic latitude and just use spherical formulae. This will be 100% conformal and works with every conformal map projection.

The "center of projection" and "azimuth" is just a graticule rotation before applying the map projection. (formula 5-7 and 5-8b of https://pubs.usgs.gov/pp/1395/report.pdf) I could reproduce the projection with alpha=30°, lambda_0=115°. Unfortunately the exact value of beta is unknown (it should be calculated using spherical trigonometry) but I could achieve an almost perfect result with beta=-28.8° (determined by trial and error).

So this should be possible.

Do you have a a working example?

[kerkovits]

Do you have a a working example?

Yes, I do but it is outside the PROJ ecosystem, and is available only in Hungarian. This program was created for personal use, so do not expect easy-to-understand program codes. I can only share the mathematical formula, I am bad at programming.

http://mercator.elte.hu/~kerkovits/projections/megnez.php?vetulet=adams&s1=115&s2=30&s3=-28.8

The image is rotated by 45° but it is not a big issue.

[kerkovits]

I found that the parameter beta in the formula I linked before should be atan2(-sin(azimuth),-tan(latitude_of_center)) This gives beta=-28.8013495219182°, which is fairly close to the value determined by trial and error.

[philippemiron]

Can we reopened this issue? And is there someone that could help me with this?

[kbevers]

@philippemiron if you are willing to take this on I can help you along as time allow it

[kerkovits]

I did some calculations and it seems that ellipsoidal formulae are a bit more complicated than I originally assumed but still possible to implement:

Step 1: calculate conformal latitude of of the center. Use formula (3-1) of https://pubs.usgs.gov/pp/1395/report.pdf substitute the latitude of center for phi, and the first eccentricity to e. I will call the result chi as the conformal_latitude_of_center. This step should be omitted for spherical calculations.

Step 2: calculate the following quantities:
alpha=asin(cos(conformal_latitude_of_center)*cos(azimuth)) (should be a bit less than +30° for ellipsoidal formulae and exactly +30° for a sphere)
lambda_0=longitude_of_center+atan2(tan(azimuth),-sin(conformal_latitude_of_center)) (should be exactly +115°)
beta=atan2(-sin(azimuth),-tan(conformal_latitude_of_center)) (should be approximately -29°)

Step 3: Calculate the conformal latitude of the mapped point by formula (3-1). Skip this step for spherical formula.

Step 4: Calculate phi_prime and lambda_prime from formulae (5-7) and (5-8b). Use constants alpha, beta, and lambda_0 from step 2 and use the conformal latitude from step 3 for phi. Lambda is the longitude of the mapped point.

Step 5: Apply the spherical Adams projection for phi_prime and lambda_prime.

Step 6: Rotate the result by 45° clockwise.

Step 7: Rejoice.

[noirchen]

@rtlemos This helps a lot.

[mdaeron]

May I ask about the current status of this feature request please? I'm not sure whether the fact that there are now implementations in js, Python and R is directly relevant to the inclusion of this projection in PROJ. In my professional community (paleoceanographers), even a spherical version would be super useful. I wish I could help but that seems beyond my expertise.

[kbevers]

May I ask about the current status of this feature request please?

This is were we keep track of the current status. If I am not mistaken, none of the regular PROJ contributors have the ressources nor the intereset to get into this topic. While I would like to see this projection implemented it is not a priority to me and I suspect it is the same for the rest of the group.

I wish I could help

You can! You can support the work by contracting someone to implement it on your behalf.

[AsgerPetersen]

To help potential sponsors better understand what contracting a developer might involve, could we come up with a very rough estimate of the effort required? Would it be on the scale of tens, hundreds, or thousands of hours? I realize this can vary significantly depending on the developer, but even an approximate indication could make it easier for someone to move forward with this.

[kbevers]

If I were to implement this tomorrow, my estimate would be around two weeks. I'm a bit rusty so it's likely someone like @rouault or @jjimenezshaw could do it faster.

[rouault]

not that I'm super eager to tackle the task, but if someone can find a budget, that would probably be closer to 3 days of my time.

[Heed725]

So like what you're saying it's possible for you guys to make a spilhaus projection in Proj4/Qgis within 72 hrs

[rouault]

within 72 hrs

I don't usually work 24 hours a day, so this is more 24 hours of work, cumulative, and not necessarily reflective of the delay :-)

[jjimenezshaw]

I have been working on this projection in this branch:
https://github.com/jjimenezshaw/PROJ/blob/spilhaus/src/projections/spilhaus.cpp

I followed the instructions above, and looked at the python implementation.
The output in QGIS looks like this:

However, using the checkpoints provided by @mdsumner

inadvertently, not seeing this comment I did this here: https://github.com/mdsumner/spilheist/tree/main/inst/extdata

I get the coordinates x and y multiplied by a factor of about 0.70847. All of checkpoints. Strange.
For instance, for the line 41.34837298,-117.28057438,12421875,-14140625, what I get is

echo -117.28057438 41.34837298 0 | ./cct +proj=spilhaus
8800603.1036  -10018297.0118        0.0000           inf

You can check the full file with this

cat grid214.csv | tail -n +2 | awk -F "," ' {print $2,$1, 0, $3, $4}' | ./build/bin/cct +proj=spilhaus | awk '{print $1/$4, $2/$5}' | sort -u

That factor is not sqrt(1/2), but not far.

Can anybody find the bug in my code? (or in ArcGIS ;)

See that this makes the images apparently correct because the scale is not noticed.

[mdsumner]

I will keep looking but I expect those checkpoints are off. (I never used them in anger I got bogged down in getting that to work).

If I use the implementation at rtlemos/spilhaus I get those same numbers as you (8800603.1036 -10018297.0118 ).

I'll do some more checks.

[jjimenezshaw]

This opens another question. Where is supposed to be Adams square II "true scale"? (a projected meter is actually a meter in the ellipsoid).
If I try this next to null island, it is not true scale:

echo 0.0 0.0 0 | ./cct +proj=adams_ws2
       0.0000         0.0000        0.0000           inf

echo 0.01 0.0 0 | ./cct +proj=adams_ws2
     556.5975         0.0000        0.0000           inf

echo 0 0 0.0 0.01 | ./geod +ellps=GRS80 -I
90d	-90d	1113.195

Maybe we are considering (involuntarily) a "scale factor" different than Esri.

[kerkovits]

The conformal latitude also introduces some distortion. I have no idea, if the implementation of Esri corrects for it. You can try to multiply both coordinates by (cos(Phi_0)/sqrt(1-e^2sin^2(Phi_0)))/cos(phi_0), where e is the first eccentricity, Phi_0 is the latitude of the centre, and phi_0 is the conformal latitude of the centre. This correction is very minor (will not solve the problem), but may give a hint for the scaling you need. I do not know whether this is necessary.

It might be a better idea to first eliminate inconsistencies of the scale on a spherical map projection, so that you can determine whether you need such a correction, too.

[jjimenezshaw]

@kerkovits that was interesting. I applied that correction, and now the difference with esri is 0.70710675... yes, that is sqrt(0.5)
See that we are doing a rotation of 45 deg between adams_ws and spilhaus, but I think the rotation is properly done.

Now the scale factor respect to a distance measured in the ellipsoid with the command geod between two points near the center of the projection is exactly 0.5. I would guess that for esri it is 0.707106.

Looking at a Tissot's indicatrix https://www.researchgate.net/figure/Tissots-indicatrix-of-the-Spilhaus-square-projection-showing-the-local-areal_fig6_371831598 you see that the minimum is at the center of the projection -ignore the crossing lines. It is an artifact of the world map crossing the antimeridian-

I did the same exercise in QGIS plotting the scale factor (being conformal, it is the same on easting and northing) for the Adams_ws2, and I got this map. The black rectangle is the line of factor = 1

The open questions are

• Do we want to apply the conformal correction? (This is easy: yes)
• Are we doing it correctly, and we simply have to decide an arbitrary scale factor?
• Which scale factor do we want to apply in Spilhaus?
• Do we want to align Spilhaus with esri? (and dis-align with other implementations)
• Do we want to change also the scale factor of Adams WS2?

[kerkovits]

@jjimenezshaw Your measurements of local scale agrees with the literature, so the implementation is probably correct.

[jjimenezshaw]

@jjimenezshaw Your measurements of local scale agrees with the literature, so the implementation is probably correct.

Thanks @kerkovits I have been looking for such literature and I couldn't find anything (my search skill are not the best ;)
Could you send me some links?

[kbevers]

@jjimenezshaw have you seen the write-up by Esri? It's an interesting read and as far as I can tell all the relevant academic references are listed at the bottom. The original papers by Spilhaus aren't exactly full of equations, though.

[rouault]

@jjimenezshaw Regarding the sqrt(2) factor difference between your implementation and ESRI's one, we probably need some scaling factor parmeter (would that be k_0 ?? or a custom one), so we can use it when mapping ESRI:54099 to PROJ string to get full interoperability with them.

[jjimenezshaw]

@rouault the protection now uses k_0 (independently of the factor we decide to use internally)

Does it make sense to also read the central latitude, central longitude and azimuth? And the rotation angle? (With the current default values)

https://spatialreference.org/ref/esri/54099/wkt.html

[rouault]

Does it make sense to also read the central latitude, central longitude and azimuth? And the rotation angle? (With the current default values)

I think that makes sense to have the parameters user-settable rather than hardcoded. No strong opinion if their default should be the boring ones (0) or the ones of ESRI;54099

[jjimenezshaw]

I think that makes sense to have the parameters user-settable rather than hardcoded.

I tried using the parameters lon_0, lat_0 and azi for the central longitude, latitude, and azimuth.
But there is some interference at least with lon_0.
If I run with +proj=spilhaus +lon_0=66.94970198 it should be the same output, but it is not.

One option is to use lon_c and lat_c. Any better idea?

[jjimenezshaw]

There is already a PR #4401 with the implementation and documentation (https://osgeo-proj--4401.org.readthedocs.build/en/4401/operations/projections/spilhaus.html)
If you have any comment, better doing it now before merging it into master.

[rouault]

Export ESRI:54099 WGS_1984_Spilhaus_Ocean_Map_in_Square as PROJ string in #4402

Regarding the sqrt(2) factor between PROJ and ESRI implementation, given that ESRI spilhaus is a particular case of adams_ws2, I would assume that the sqrt(2) factor would also exist for regular adams_ws2 . Could someone provide test point (ie longitude/latitude -> easting/northing) for ESRI:54098 WGS_1984_Adams_Square_II ?

[jjimenezshaw]

Regarding the sqrt(2) factor between PROJ and ESRI implementation, given that ESRI spilhaus is a particular case of adams_ws2, I would assume that the sqrt(2) factor would also exist for regular adams_ws2 . Could someone provide test point (ie longitude/latitude -> easting/northing) for ESRI:54098 WGS_1984_Adams_Square_II ?

Just to track it. https://github.com/kgjenkins/test-coordinates