SciTools
`Spilhaus` projection not working correctly when plotting raster data
I'm trying to plot raster data over a Spilhaus projection but I'm getting some artefacts
This code is MWE to produce the image
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(),
)
when trying with imshow I get this instead, I guess the data is not properly rotated (but no artefacts).
I think that what might be going on is that for pcolormesh we only project each of the corners and draw straight lines between them. For this projection there is quite a bit of distortion and my assumption is that one of the corners may be transformed to the other side of the projection. So the projected square has one corner on the opposite side of the map. We try to get around that by identifying "long" cells, but the rotation here might prevent that from being caught / registered. It would be nice to have a more robust way of identifying these cases. https://github.com/SciTools/cartopy/blob/aedaad977556e0ad15b9256e9cafd15913108fd4/lib/cartopy/mpl/geoaxes.py#L1845-L1863
Here is a simpler reproducer without xarray/downloading showing what I mean with the one weird-looking point and a minimal number of cells.
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
lons = np.linspace(-160, 160, 2)
lats = np.linspace(-30, 30, 2)
lon2d, lat2d = np.meshgrid(lons, lats)
z = lon2d + lat2d
# Set up the figure and the projection for the plot
fig, ax = plt.subplots(
figsize=(10, 10),
subplot_kw={
"projection": ccrs.Spilhaus()
},
)
ax.coastlines()
ax.set_global()
ax.pcolormesh(
lon2d,
lat2d,
z,
transform=ccrs.PlateCarree(),
)
plt.show()
Here is how I created the mask for plotting. Obviously, this masks a bit more than the cartopy default, but it should be fine for high resolution situations.
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
def get_spilhaus_mask(lon2d,lat2d):
x = projection.transform_points(ccrs.PlateCarree(),lon2d,lat2d)
y = x[...,1]
x = x[...,0]
bound = np.mean(abs(np.array(projection.bounds)))
near_bound = np.logical_or(abs(x)>0.9*bound, abs(y)>0.9*bound)
sign_change = np.zeros_like(x)
change_sign0 = np.logical_or(x[1:, 1:]*x[:-1, :-1]<0,
y[1:, 1:]*y[:-1, :-1]<0)
change_sign1 = np.logical_or(x[1:, :-1]*x[:-1, 1:]<0,
y[1:, :-1]*y[:-1, 1:]<0)
change_sign = np.logical_or(change_sign0,change_sign1)
sign_change[1:,1:] = change_sign
sign_change[:-1,:-1] = np.logical_or(sign_change[:-1,:-1],change_sign)
sign_change[1:,:-1] = np.logical_or(sign_change[1:,:-1],change_sign)
sign_change[:-1,1:] = np.logical_or(sign_change[:-1,1:],change_sign)
return np.logical_and(near_bound,sign_change)
lons = np.linspace(-180, 180, 100)
lats = np.linspace(-90, 90, 100)
lon2d, lat2d = np.meshgrid(lons, lats)
mask = get_spilhaus_mask(lon2d,lat2d)
z = lon2d + lat2d
z[mask] = np.nan
# Set up the figure and the projection for the plot
fig, ax = plt.subplots(
figsize=(10, 10),
subplot_kw={
"projection": Spilhaus()
},
)
ax.coastlines()
ax.set_global()
ax.pcolormesh(
lon2d,
lat2d,
z,
transform=ccrs.PlateCarree(),
)
plt.show()
It is on my todo list to come up with a way to do this properly, which involves utilizing the tiling property of the projection and add ghost points just outside the domain. This has been on the back burner though, since I am trying to graduate...
Here is how I created the mask for plotting. Obviously, this masks a bit more than the cartopy default, but it should be fine for high resolution situations.
import cartopy.crs as ccrs import cartopy.feature as cfeature import matplotlib.pyplot as plt import numpy as np def get_spilhaus_mask(lon2d,lat2d): x = projection.transform_points(ccrs.PlateCarree(),lon2d,lat2d) y = x[...,1] x = x[...,0] bound = np.mean(abs(np.array(projection.bounds))) near_bound = np.logical_or(abs(x)>0.9*bound, abs(y)>0.9*bound) sign_change = np.zeros_like(x) change_sign0 = np.logical_or(x[1:, 1:]*x[:-1, :-1]<0, y[1:, 1:]*y[:-1, :-1]<0) change_sign1 = np.logical_or(x[1:, :-1]*x[:-1, 1:]<0, y[1:, :-1]*y[:-1, 1:]<0) change_sign = np.logical_or(change_sign0,change_sign1) sign_change[1:,1:] = change_sign sign_change[:-1,:-1] = np.logical_or(sign_change[:-1,:-1],change_sign) sign_change[1:,:-1] = np.logical_or(sign_change[1:,:-1],change_sign) sign_change[:-1,1:] = np.logical_or(sign_change[:-1,1:],change_sign) return np.logical_and(near_bound,sign_change) lons = np.linspace(-180, 180, 100) lats = np.linspace(-90, 90, 100) lon2d, lat2d = np.meshgrid(lons, lats) mask = get_spilhaus_mask(lon2d,lat2d) z = lon2d + lat2d z[mask] = np.nan # Set up the figure and the projection for the plot fig, ax = plt.subplots( figsize=(10, 10), subplot_kw={ "projection": Spilhaus() }, ) ax.coastlines() ax.set_global() ax.pcolormesh( lon2d, lat2d, z, transform=ccrs.PlateCarree(), ) plt.show()
It is on my todo list to come up with a way to do this properly, which involves utilizing the tiling property of the projection and add ghost points just outside the domain. This has been on the back burner though, since I am trying to graduate...
Hey @MaceKuailv, thanks a lot for this!
I tried on my data and it does indeed seem to make things better.
However, there is still some data outside of the bounds, and I think this was also present in your example
Luckily most of those lines are on land, so it's not a big deal for me, but some of those appear in the lower left and upper right corner, and I would like to get rid of those :)
In my completely ignorance of your logic I tried to decrease the value used here
near_bound = np.logical_or(abs(x)>0.9*bound, abs(y)>0.9*bound)
but it didn't make a difference.
Do you know how I could tweak the formula to make a more aggressive masking?
BTW, projection in your example is ccrs.Spilhaus(), right?
@guidocioni I didn't realize that I missed the very crucial step of rotating the domain. Sorry about that.
I also wrote some inline comments to explain what I have done.
def get_spilhaus_mask(lon2d,lat2d):
x = projection.transform_points(ccrs.PlateCarree(),lon2d,lat2d)
# Get the projected picture
y = x[...,1]
x = x[...,0]
bound = np.mean(abs(np.array(projection.bounds)))
# The center of the projected domain is (0,0),
# This is selecting points near the boundary.
# As you will see, not all of this is masked.
near_bound = np.logical_or(abs(x)>0.95*bound, abs(y)>0.95*bound)
# After rotating 45 degrees, the four edges of the domain are
# in four different quadrants.
theta = np.pi/4
rot_x = x*np.cos(theta)-y*np.sin(theta)
rot_y = y*np.cos(theta)+x*np.sin(theta)
# We now need to Figure out which points and their neighbors
# are at different quadrants.
# We first look at the corners of the grid and figure out
# which ones are at those boundaries.
change_sign0 = np.logical_or(rot_x[1:, 1:] *rot_x[:-1, :-1]<0,
rot_y[1:, 1:] *rot_y[:-1, :-1]<0)
change_sign1 = np.logical_or(rot_x[1:, :-1]*rot_x[:-1, 1:]<0,
rot_y[1:, :-1]*rot_y[:-1, 1:]<0)
change_sign = np.logical_or(change_sign0,change_sign1)
# If a grid cell touches one of the selected corner, we mask it.
# In other words, if the neighboring 3*3 grid points experience quadrants change,
# we will mask the center grid point.
sign_change = np.zeros_like(x)
sign_change[1:,1:] = change_sign
sign_change[:-1,:-1] = np.logical_or(sign_change[:-1,:-1],change_sign)
sign_change[1:,:-1] = np.logical_or(sign_change[1:,:-1],change_sign)
sign_change[:-1,1:] = np.logical_or(sign_change[:-1,1:],change_sign)
# Change of quandrants can happen at the four edges or
# at x=0 or y=0. We only select those sufficiently close to the boundary.
# This method will mask some points near the four corners, but those are all on land.
return np.logical_and(near_bound,sign_change)
@guidocioni Your SST map also works now:
I am using jet, which does not look as good as your color map...
Hey @MaceKuailv , first of all thanks a lot for getting back to me. In the meantime I used Claude to generate a new version of the masking function based on your example and it got me close to your second solution (unbelievable, sometimes AI can really be useful :D). I will try also your solution once I have some time.
The final result looks something like this. I wanted to make an animation to show equatorial waves in this projection so that I can lose 5 minutes trying to understand the geometry in Spilhaus :D
https://github.com/user-attachments/assets/2c5f46b1-929c-4b03-9e4c-4a6c64d24532
Oh wow, this is amazing. Can you tell me what color map are you using? It looks really good and have great contrast. Can you also show us what Claude came up with?
Oh wow, this is amazing. Can you tell me what color map are you using? It looks really good and have great contrast. Can you also show us what Claude came up with?
it's part of a series of colormaps I manually created based on figures found online or just personal experience. They're more tailored towards meteorological data.
I always wanted to make a repo just for those but was always too lazy so I'm still copy pasting them every time i need them. You can find them in one of my last projects where I'm trying to finally collect them all together: https://github.com/guidocioni/icon-2i in the colormaps folder. There's some details in the README on how to use them and if you check the plotting scripts there are also some levels that work well with them .Have fun!
@guidocioni I can say we'd be happy to have the colormaps in MetPy.
Hey @dopplershift , I'd be happy to draft a PR. Can you point me towards the part of the code where the colormaps are? I use metpy daily but never to import colormaps 😄
Love to see this collaboration in the open! ❤ Great stuff everyone.
A few comments:
- Maybe add an example to Cartopy's gallery to demonstrate the new projection and the work you're doing with a potentially helpful masking function for others to take advantage of.
- Can the masking function be generalized for pcolormesh more generally so that we would get that for all users on any projection?
- Can the masking function be generalized for pcolormesh more generally so that we would get that for all users on any projection?
@greglucas this is actually something that puzzled me since i started using cartopy.
Is the data outside of the projection footprint ever masked when plotting with pclormesh, contourf...? For large arrays I've always seen a big difference in plotting time when masking (or not) the input data, which always made me think that cartopy is trying to plot the data also where is not actually visible. But I always had to do it manually beforehand, and I usually employ a simple slicing on the coordinates of the input array using the extent of the target projection converted into plate carree, which is quite a raw approach.
Love to see this collaboration in the open! ❤ Great stuff everyone.
A few comments:
* Maybe add an example to Cartopy's gallery to demonstrate the new projection and the work you're doing with a potentially helpful masking function for others to take advantage of. * Can the masking function be generalized for pcolormesh more generally so that we would get that for all users on any projection?
@guidocioni Can you share the masking function that Claude came up with?
I think the general way of determining the mask involves calculating whether the discontinuity set (a great circle path in this case and many other cases) intersect each individual mesh. This may involve adding attribute to each projection manually, and maybe not very fast.
Love to see this collaboration in the open! ❤ Great stuff everyone. A few comments:
* Maybe add an example to Cartopy's gallery to demonstrate the new projection and the work you're doing with a potentially helpful masking function for others to take advantage of. * Can the masking function be generalized for pcolormesh more generally so that we would get that for all users on any projection?@guidocioni Can you share the masking function that Claude came up with?
I think the general way of determining the mask involves calculating whether the discontinuity set (a great circle path in this case and many other cases) intersect each individual mesh. This may involve adding attribute to each projection manually, and maybe not very fast.
I don't think Claude`s way is smartern than yours, but it worked for me once I reduced the boundary. I think it only cuts a rectangular box, but if you make it 'large' enough you don't end up seeing it. Again, it may be nonsense, but it worked for me and I didn't have enough time to debug, so here goes nothing :)
def get_spilhaus_mask(lon2d, lat2d):
# Transform coordinates to Spilhaus projection
x = ccrs.Spilhaus().transform_points(ccrs.PlateCarree(), lon2d, lat2d)
y = x[..., 1]
x = x[..., 0]
# Get projection bounds
bound = np.mean(abs(np.array(ccrs.Spilhaus().bounds)))
# Increase boundary threshold and add distance-based masking
near_bound = np.logical_or(abs(x) > 0.999 * bound, abs(y) > 0.999 * bound)
# Calculate distances to edges
edge_dist = np.minimum(
np.minimum(abs(x - bound), abs(x + bound)),
np.minimum(abs(y - bound), abs(y + bound)),
)
edge_mask = edge_dist < bound * 0.001
# Original sign change detection
sign_change = np.zeros_like(x)
change_sign0 = np.logical_or(
x[1:, 1:] * x[:-1, :-1] < 0, y[1:, 1:] * y[:-1, :-1] < 0
)
change_sign1 = np.logical_or(
x[1:, :-1] * x[:-1, 1:] < 0, y[1:, :-1] * y[:-1, 1:] < 0
)
change_sign = np.logical_or(change_sign0, change_sign1)
# Expand sign change mask
sign_change[1:, 1:] = change_sign
sign_change[:-1, :-1] = np.logical_or(sign_change[:-1, :-1], change_sign)
sign_change[1:, :-1] = np.logical_or(sign_change[1:, :-1], change_sign)
sign_change[:-1, 1:] = np.logical_or(sign_change[:-1, 1:], change_sign)
# Combine masks with additional edge detection
final_mask = np.logical_or(np.logical_and(near_bound, sign_change), edge_mask)
# Add additional buffer around masked regions
final_mask = np.maximum.reduce(
[
final_mask,
np.roll(final_mask, 1, axis=0),
np.roll(final_mask, -1, axis=0),
np.roll(final_mask, 1, axis=1),
np.roll(final_mask, -1, axis=1),
]
)
return final_mask
Is the data outside of the projection footprint ever masked when plotting with pclormesh, contourf...? For large arrays I've always seen a big difference in plotting time when masking (or not) the input data, which always made me think that cartopy is trying to plot the data also where is not actually visible. But I always had to do it manually beforehand, and I usually employ a simple slicing on the coordinates of the input array using the extent of the target projection converted into plate carree, which is quite a raw approach.
Sorry, sloppy choice of words on my part. The data isn't masked within Cartopy because we need to transform everything to see where it ends up in the map projection and enable setting smaller/larger limits. So what you are doing to speed things up makes sense when that is possible.
To clarify what I meant though, in pcolormesh we do "mask" some cells that cross the dateline of projections and wrap around the left to right side of the map. Those show up as these streaks that you see here and we try to hide them and instead draw those point with a pcolor() call for just those "masked" cells. That is something to try here, do a pcolor() call instead of pcolormesh() and see if that fixes those streaks or not.
Hi everyone, I'm trying to test the projection with the latest cartopy release. I'm not sure what I'm doing wrong, but I'm getting the following error message. Someone has any pointers? Cheers!
uv init spilhaus
cd spilhaus
uv add "cartopy>=0.25"
uv run python -c "import cartopy.crs as ccrs; ccrs.Spilhaus()"
Initialized project `spilhaus` at `~/repos/spilhaus`
Using CPython 3.12.3
Creating virtual environment at: .venv
Resolved 17 packages in 6ms
Installed 16 packages in 68ms
+ cartopy==0.25.0
+ certifi==2025.8.3
+ contourpy==1.3.3
+ cycler==0.12.1
+ fonttools==4.59.0
+ kiwisolver==1.4.8
+ matplotlib==3.10.5
+ numpy==2.3.2
+ packaging==25.0
+ pillow==11.3.0
+ pyparsing==3.2.3
+ pyproj==3.7.1
+ pyshp==2.3.1
+ python-dateutil==2.9.0.post0
+ shapely==2.1.1
+ six==1.17.0
Traceback (most recent call last):
File "<string>", line 1, in <module>
File "~/repos/spilhaus/.venv/lib/python3.12/site-packages/cartopy/crs.py", line 3257, in __init__
super().__init__(proj4_params, globe=globe)
File "~/repos/spilhaus/.venv/lib/python3.12/site-packages/cartopy/crs.py", line 677, in __init__
super().__init__(*args, **kwargs)
File "~/repos/spilhaus/.venv/lib/python3.12/site-packages/cartopy/crs.py", line 203, in __init__
super().__init__(self.proj4_init)
File "~/repos/spilhaus/.venv/lib/python3.12/site-packages/pyproj/crs/crs.py", line 350, in __init__
self._local.crs = _CRS(self.srs)
^^^^^^^^^^^^^^
File "pyproj/_crs.pyx", line 2364, in pyproj._crs._CRS.__init__
pyproj.exceptions.CRSError: Invalid projection: +proj=spilhaus +ellps=WGS84 +rot=45 +x_0=0.0 +y_0=0.0 +no_defs +type=crs: (Internal Proj Error: proj_create: Error 1027 (Invalid value for an argument): Unknown projection)
Hi everyone, I'm trying to test the projection with the latest cartopy release. I'm not sure what I'm doing wrong, but I'm getting the following error message. Someone has any pointers? Cheers!
uv init spilhaus cd spilhaus uv add "cartopy>=0.25" uv run python -c "import cartopy.crs as ccrs; ccrs.Spilhaus()"Initialized project `spilhaus` at `~/repos/spilhaus` Using CPython 3.12.3 Creating virtual environment at: .venv Resolved 17 packages in 6ms Installed 16 packages in 68ms + cartopy==0.25.0 + certifi==2025.8.3 + contourpy==1.3.3 + cycler==0.12.1 + fonttools==4.59.0 + kiwisolver==1.4.8 + matplotlib==3.10.5 + numpy==2.3.2 + packaging==25.0 + pillow==11.3.0 + pyparsing==3.2.3 + pyproj==3.7.1 + pyshp==2.3.1 + python-dateutil==2.9.0.post0 + shapely==2.1.1 + six==1.17.0 Traceback (most recent call last): File "<string>", line 1, in <module> File "~/repos/spilhaus/.venv/lib/python3.12/site-packages/cartopy/crs.py", line 3257, in __init__ super().__init__(proj4_params, globe=globe) File "~/repos/spilhaus/.venv/lib/python3.12/site-packages/cartopy/crs.py", line 677, in __init__ super().__init__(*args, **kwargs) File "~/repos/spilhaus/.venv/lib/python3.12/site-packages/cartopy/crs.py", line 203, in __init__ super().__init__(self.proj4_init) File "~/repos/spilhaus/.venv/lib/python3.12/site-packages/pyproj/crs/crs.py", line 350, in __init__ self._local.crs = _CRS(self.srs) ^^^^^^^^^^^^^^ File "pyproj/_crs.pyx", line 2364, in pyproj._crs._CRS.__init__ pyproj.exceptions.CRSError: Invalid projection: +proj=spilhaus +ellps=WGS84 +rot=45 +x_0=0.0 +y_0=0.0 +no_defs +type=crs: (Internal Proj Error: proj_create: Error 1027 (Invalid value for an argument): Unknown projection)
You need proj>= 9.6.0. See here https://github.com/SciTools/cartopy/issues/1376.
I'm going to close this issue as the original problem has been (kind of) resolved.