Spilhaus projection #2529

MaceKuailv · 2025-05-10T20:25:55Z

Rationale

Implementing a cool-looking map projection for global ocean called the Spilhaus Adams world in a square II map projection.
This is in response to projection request #1376

Example

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

fig = plt.figure(figsize=(8, 8))

rots = [-45,45,-135,135]
for i in range(4):
    ax = plt.subplot(2,2,i+1,projection = ccrs.Spilhaus(rotation=rots[i]))
    ax.add_feature(cfeature.OCEAN,color = 'royalblue')
#     ax.axis('off')
    
plt.tight_layout()
plt.subplots_adjust(wspace=0.0, hspace=0.0)

Many thanks to @greglucas.

…s a http error when I first run the test

CLAassistant · 2025-05-10T20:26:03Z

All committers have signed the CLA.

MaceKuailv · 2025-05-10T20:45:24Z

Some notes:

I have tried both wkt and proj4 initialization, they creates the identical map projection, and they are both very slow (about 30 seconds).
The boundary listed here is very far off. That boundary is a tangential line to 180. I then use some optimization algorithm to find the actual boundary. The meta pole is not exactly 115E 30N, but very slightly off [115.00000024, 30.00000036]. The projected domain is also not a perfect square, and the aspect ratio is off by about 1/10,000.
The tests are all failing because the map projection is too new. All of the tests pass on my local machine with PROJ 9.6.0 and pyproj 3.7.1.

MaceKuailv · 2025-05-10T20:50:36Z

https://github.com/SciTools/cartopy/blob/1e1583334587ae16644ebb31f1b311b2ef36f23f/.github/workflows/ci-testing.yml#L43C9-L45C1
Creating the Spilhaus projection object requires newer version of pyproj. I think this minimum version requirement means all tests run on 3.3.1.
I am going to add a "skipif", but this means the new tests simply won't run.

rcomer · 2025-05-10T21:06:59Z

The minimum version specification only affects one of the test runners. The others should use the latest pyproj from PyPI, though it looks like that is older than Proj 9.6. So maybe we just need to wait for a pyproj release. Did you create your local environment with conda? Conda-forge does its own re-building independent of the package maintainers.

MaceKuailv · 2025-05-10T21:14:19Z

@rcomer Thanks for clarifying. Yes, I guess we will have to wait for a pip release from pyproj. Yes, I built a new environment using conda env create -f environment.yml.

OK, it seems like the only one I passed is the minimum requirement one...

greglucas

Nice! Glad you're having fun with this one :) I've got a few suggestions to mostly make this simpler on the Cartopy maintenance side and then users can extend it later with the extra names if they want.

lib/cartopy/crs.py

lib/cartopy/tests/crs/helpers.py

lib/cartopy/tests/crs/test_spilhaus.py

lib/cartopy/crs.py

Co-authored-by: Greg Lucas <greg.m.lucas@gmail.com>

lib/cartopy/crs.py

greglucas · 2025-05-11T20:41:44Z

lib/cartopy/crs.py

+        self.bounds = [
+            -11802684.083372328,
+            11802683.949222516,
+            -11801129.925928915,
+            11801129.925928915
+        ]


Should these be an actual square and therefore the same value in all positions?

Is there a pyproj way to get these exactly with pyproj_crs.area_of_use.bounds?

Well, they should be a actual square, but it is off by about 0.01%. I think this slight inaccuracy stem from upstream. The feature I love the most about this projection is that it is a square and can be used to tile the entire surface.

The Spilhaus projection is a global projection, but (-180,-90) and (180,90) are not at the corners. The transformed bounds defined here also does not produce the correct result some how...Theoretically the points furthest away from the center should be (115,30) but it is actually slightly off. My strategy was to use scipy minimize to find the actual pole, but I don't think it make sense to run the optimization every time we open the projection object.

It is a bit unsatisfying indeed. I don't know the math going into the map projection well enough to calculate the correct lon_0,lat_0, and azimuth, but may be I could design a optimization process to find them. Either way, we will have to pass a hard coded number in to pyproj.

What happens if you do the area_of_use approach above?
CRS.from_proj4("+proj=...").area_of_use.bounds

I agree we don't want to run this on every initialization and hardcoding some values would be useful. My overall comment is that you've chosen some values that are close to being a square, so do we want to just make those a square arbitrarily since you are somewhat choosing arbitrary initial points anyways.

I wonder if the values being different from epsg.io are due to what you mention with the lon_0, lat_0 not being 0 as is typical in many other projections.

I just get a None object if I pass in the proj4 string. We can make it a perfect square, but that either leave some empty space or make the projection of some points outside the bounds. If we were to make it a square, I suggest we use the largest absolute value (11802684.083372328). However, I still much prefer that we just leave it slightly imperfect.

It is unclear what the bounds in epsg.io comes from. In any case, the correct transformed value should look something like a square. This is what it looks like when I apply the bounds on epsg.io:

This is what it looks like when I set lon_0, and lat_0 to zero under the same bounds

lib/cartopy/tests/crs/test_spilhaus.py

Co-authored-by: Greg Lucas <greg.m.lucas@gmail.com>

philippemiron · 2025-05-12T15:55:03Z

Should we also add this projection to docs/source/reference/projections.rst ? Ps: thanks for taking the lead on that @MaceKuailv !!

MaceKuailv · 2025-05-12T16:03:21Z

Should we also add this projection to docs/source/reference/projections.rst ? Ps: thanks for taking the lead on that @MaceKuailv !!

Thank you! @philippemiron

Yeah, we can definitely do that. My only concern is that between we merge this PR and making a public release, users will see this projection in the documentation and no be able to use it without cloning the repo. I could open another PR just to add this to the documentation. Let me know what the preferred work flow is.

greglucas · 2025-05-12T19:21:38Z

Yeah, we can definitely do that. My only concern is that between we merge this PR and making a public release, users will see this projection in the documentation and no be able to use it without cloning the repo. I could open another PR just to add this to the documentation. Let me know what the preferred work flow is.

The documentation is only released when we make a code release, so no problem there. You can run that make_projections script here or in a follow-up PR. But, I agree it would be good to run that before the next release. 👍

greglucas

Overall looks good 👍 a few more doc updates.

lib/cartopy/tests/crs/test_spilhaus.py

lib/cartopy/tests/mpl/test_mpl_integration.py

Co-authored-by: Greg Lucas <greg.m.lucas@gmail.com>

greglucas

👍 Looks good to me. One more minor formatting change and a request to add a link to the ArcGIS story or something similar to explain the projection a little more.

docs/make_projection.py

greglucas · 2025-05-15T03:31:06Z

lib/cartopy/crs.py

@@ -3219,6 +3219,40 @@ def x_limits(self):
    def y_limits(self):
        return self._y_limits

+class Spilhaus(Projection):
+    """
+    Spilhaus World Ocean Map in a Square


I think we should add a link to the story about this projection which explains some of the ideas here too and is a nice historical reference.
https://storymaps.arcgis.com/stories/756bcae18d304a1eac140f19f4d5cb3d

In particular it explains why those two points you chose below are special. Maybe even put those antipodal points int he docstring here too?

Co-authored-by: Greg Lucas <greg.m.lucas@gmail.com>

greglucas

LGTM, thanks for incorporating this!

philippemiron · 2025-05-16T18:35:37Z

@greglucas do you know when/how new version of cartopy are released? Cheers and have a good weekend!

greglucas · 2025-05-22T02:27:40Z

@greglucas do you know when/how new version of cartopy are released?

We don't have any set schedule, but it has been every roughly 6 months - 1 year. It looks like we are about due for another release and have quite a few things ready to go. I opened this issue #2532 to track any features folks want to see get in. Feel free to leave a comment there if you have anything else you'd like to see in!

guidocioni · 2025-05-27T07:11:05Z

Hey guys,
I was having fun with this new projection today trying to plot some maps of SST but I noticed there is something weird going on...

Do you know what could cause those artefacts?

rcomer · 2025-05-27T09:50:20Z

@guidocioni is that a contour plot and did you use transform_first? If so it could be the issue discussed at #2233

guidocioni · 2025-05-27T09:59:04Z

@guidocioni is that a contour plot and did you use transform_first? If so it could be the issue discussed at #2233

That's the first thing I thought about, but no, it is a pcolormesh. At this resolution a contour plot would take ages.

I've seen this behaviour before and it happens on non-standard projections sometimes, but I'm not sure what it is exactly.

For reference, this is the code you can use to reproduce it (provided you have of course a cartopy version with this PR already inside). You can run it on your side because it automatically downloads (and cache) the data needed for the plot

import xarray as xr
import fsspec
import pandas as pd
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import time
from aiohttp.client_exceptions import ServerDisconnectedError

date = pd.Timestamp("2024-01-01")
url = f"https://www.star.nesdis.noaa.gov/pub/socd/mecb/crw/data/5km/v3.1_op/nc/v1.0/daily/sst/{date.year}/coraltemp_v3.1_{date.strftime('%Y%m%d')}.nc"
while True:
    try:
        file = fsspec.open_local(
            f"simplecache::{url}", simplecache={"cache_storage": "/tmp/"}
        )
        ds = xr.open_dataset(file, engine="netcdf4").squeeze()
        break
    except ServerDisconnectedError:
        time.sleep(0.5)

# Set up the figure and the projection for the plot
fig, ax = plt.subplots(
    figsize=(10, 10),
    subplot_kw={
        "projection": ccrs.Spilhaus()
    },
)

ax.add_feature(
    cfeature.LAND.with_scale("110m"), edgecolor="black", linewidth=0.25, zorder=2
)
lon2d, lat2d = np.meshgrid(ds["lon"].values, ds["lat"].values)

ax.pcolormesh(
    lon2d,
    lat2d,
    ds["analysed_sst"].values,
    transform=ccrs.PlateCarree(),
)

rcomer · 2025-05-27T10:09:40Z

@guidocioni I suggest posting your example as a new issue.

guidocioni · 2025-05-27T10:10:54Z

@guidocioni I suggest posting your example as a new issue.

I wanted to but does it make sense as this has not yet been published into a release? Otherwise I'm happy to do so :)

MaceKuailv · 2025-05-27T12:44:07Z

@guidocioni I ran into this issue as well. It is essential the same issue when you try to plot across datelines. Since the border here is quite a bit harder to define than “180W”, this makes the plotting tricky.

I have a solution that involves figuring out the boundaries and masking the boundary grid points with NaN. It is a pretty ad hoc solution, but it works OK for high resolution simulations. I can share an example in your new issue.

guidocioni · 2025-05-27T13:04:35Z

@guidocioni I ran into this issue as well. It is essential the same issue when you try to plot across datelines. Since the border here is quite a bit harder to define than “180W”, this makes the plotting tricky.

I have a solution that involves figuring out the boundaries and masking the boundary grid points with NaN. It is a pretty ad hoc solution, but it works OK for high resolution simulations. I can share an example in your new issue.

Yes, that would definitely help! We can move the discussion into #2542

MaceKuailv added 5 commits May 10, 2025 11:49

change import from numpy to libc just to pip install -e, p.s. there i…

b54b56a

…s a http error when I first run the test

create the projection, but it is really slow

e5ec3f0

verify that proj4 did not introduce extra error and implement options

9bf54f9

add test for mpl_integration

7db6e14

add test image, fix a very embarrassing typo

1c66194

MaceKuailv added 3 commits May 10, 2025 16:26

add new option for greenland_at

0f068c3

style fix

4bf4f9c

resolve conflict

3ecbd81

skip tests if pyproj too old

8f574de

greglucas reviewed May 10, 2025

View reviewed changes

also skip based on PROJ version

83ce837

rcomer linked an issue May 10, 2025 that may be closed by this pull request

Projection request: Spilhaus projection (world ocean map) #1376

Closed

MaceKuailv and others added 5 commits May 10, 2025 18:22

Styling of documentation

5551006

Co-authored-by: Greg Lucas <greg.m.lucas@gmail.com>

fix based on greg's comment

65f8096

not skip for pyproj

f85997d

conflict resolve

533df1b

add a empty line after doc string, don't know why

e554ca8

greglucas reviewed May 11, 2025

View reviewed changes

lib/cartopy/tests/crs/test_spilhaus.py Outdated Show resolved Hide resolved

MaceKuailv and others added 5 commits May 12, 2025 11:12

reduce proj4 argument passed in.

d446388

Co-authored-by: Greg Lucas <greg.m.lucas@gmail.com>

add globe argument

234808c

Co-authored-by: Greg Lucas <greg.m.lucas@gmail.com>

pass globe into __init__

d6cb858

Co-authored-by: Greg Lucas <greg.m.lucas@gmail.com>

fix skip reason

a848abf

Co-authored-by: Greg Lucas <greg.m.lucas@gmail.com>

fix test docstring "should have read what I coppy and paste"

bf30dcc

Co-authored-by: Greg Lucas <greg.m.lucas@gmail.com>

greglucas reviewed May 12, 2025

View reviewed changes

lib/cartopy/tests/crs/test_spilhaus.py Outdated Show resolved Hide resolved

lib/cartopy/tests/mpl/test_mpl_integration.py Outdated Show resolved Hide resolved

MaceKuailv and others added 4 commits May 12, 2025 15:44

fix skip reason

2cdcf84

Co-authored-by: Greg Lucas <greg.m.lucas@gmail.com>

fix skip reason

af8dfc0

Co-authored-by: Greg Lucas <greg.m.lucas@gmail.com>

run make_projections

092bc50

style

80b466f

greglucas approved these changes May 15, 2025

View reviewed changes

MaceKuailv and others added 3 commits May 15, 2025 11:08

Style in make projection

109d2fe

Co-authored-by: Greg Lucas <greg.m.lucas@gmail.com>

update spilhaus description

fa82706

typo and style

40dbfa6

greglucas approved these changes May 15, 2025

View reviewed changes

greglucas merged commit 7c66fd5 into SciTools:main May 15, 2025
23 checks passed

QuLogic added this to the Next Release milestone May 16, 2025

greglucas mentioned this pull request May 23, 2025

TST/FIX: Spilhaus doesn't contain units or datum by default, so remove #2535

Merged

Spilhaus projection #2529

Spilhaus projection #2529

Uh oh!

Conversation

MaceKuailv commented May 10, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Rationale

Example

Uh oh!

CLAassistant commented May 10, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

MaceKuailv commented May 10, 2025

Uh oh!

MaceKuailv commented May 10, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

rcomer commented May 10, 2025

Uh oh!

MaceKuailv commented May 10, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

greglucas left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

greglucas May 11, 2025

Choose a reason for hiding this comment

Uh oh!

MaceKuailv May 12, 2025

Choose a reason for hiding this comment

Uh oh!

greglucas May 12, 2025

Choose a reason for hiding this comment

Uh oh!

MaceKuailv May 12, 2025

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

philippemiron commented May 12, 2025

Uh oh!

MaceKuailv commented May 12, 2025

Uh oh!

greglucas commented May 12, 2025

Uh oh!

greglucas left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

greglucas left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

greglucas May 15, 2025

Choose a reason for hiding this comment

Uh oh!

greglucas left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

philippemiron commented May 16, 2025

Uh oh!

greglucas commented May 22, 2025

Uh oh!

guidocioni commented May 27, 2025

Uh oh!

rcomer commented May 27, 2025

Uh oh!

guidocioni commented May 27, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

MaceKuailv commented May 10, 2025 •

edited

Loading

CLAassistant commented May 10, 2025 •

edited

Loading

MaceKuailv commented May 10, 2025 •

edited

Loading

MaceKuailv commented May 10, 2025 •

edited

Loading

guidocioni commented May 27, 2025 •

edited

Loading