Inconsistent definition of inside/outside Nightshade feature geometry

[Eliam76]

Description

For a project, I want to plot the intersection between the Earth and specific planes facing the Sun. For this, I use the Nightshade feature : I take advantage of the refraction parameter to adjusts the day/night geometry boundary, which has the same result as shifting a plane facing the Sun. I use stereographic projection, so all the resulting projected polygons are discs or circular holes.

However, the behavior of the feature is inconsistent when the refraction value is positive (normally the value is zero or negative), especially at dates close to summer solstices (but not always at the exact solstice). The definition of the inside/outside of the polygon seems randomly wrong, as if the points weren't ordered correctly.

These are pretty extreme values for the refraction but it may underline an issue with the source code. As a suggestion, since the Nightshade feature compute the solar position
lat, lon = _solar_position(date)
it should be possible to determine the correct orientation as the projected polygon (representing the shaded part of the Earth) should be the one not containing this point.

Code to reproduce

from datetime import datetime

import matplotlib.pyplot as plt

import cartopy.crs as ccrs
from cartopy.feature.nightshade import Nightshade

date = datetime(2024, 6, 11)
strdate = date.strftime("%d/%m/%Y, %H:%M")

proj = ccrs.Stereographic(90., 0.0)

fig, axs = plt.subplots(3, 4, figsize=(12 , 9), subplot_kw={'projection': proj})

for i, ax in enumerate(axs.flatten()):
    nightshade = Nightshade(date, refraction=  i)

    ax.add_feature(nightshade)
    ax.coastlines()
    ax.set_global()
    ax.set_title(f"Refraction={i}")

fig.suptitle(f'{strdate}', fontsize = 20.)

plt.show()

Full environment definition

Operating system

Windows 10 Family, 64bits edition

Cartopy version

0.23.0

pip list

py -3.12 -m pip list

Package           Version
----------------- -----------
asttokens         2.4.1
Cartopy           0.23.0
certifi           2024.8.30
colorama          0.4.6
comm              0.2.2
contourpy         1.3.0
cycler            0.12.1
debugpy           1.8.5
decorator         5.1.1
executing         2.0.1
fonttools         4.53.1
ipykernel         6.29.5
ipython           8.26.0
jedi              0.19.1
jupyter_client    8.6.2
jupyter_core      5.7.2
kiwisolver        1.4.7
matplotlib        3.9.2
matplotlib-inline 0.1.7
nest-asyncio      1.6.0
numpy             2.1.1
packaging         24.1
parso             0.8.4
pillow            10.4.0
pip               24.2
platformdirs      4.2.2
prompt_toolkit    3.0.47
psutil            6.0.0
pure_eval         0.2.3
Pygments          2.18.0
pyparsing         3.1.4
pyproj            3.6.1
pyshp             2.3.1
python-dateutil   2.9.0.post0
pywin32           306
pyzmq             26.1.0
shapely           2.0.6
six               1.16.0
stack-data        0.6.3
tornado           6.4.1
traitlets         5.14.3
tzfpy             0.15.6
wcwidth           0.2.13

[rcomer]

Thanks for the clear report @Eliam76. Everything seems fine if we plot in the projection in which Nightshade is defined, so it looks like this is an instance of problems in the transforms.

dummy_nightshade = Nightshade()
proj = dummy_nightshade.crs

[Eliam76]

Thank you for your answer. Yes it seems that there are many projections in which the Nightshade feature behaves correctly. I tried to "manually" project the geometry using the following modification of the loop :

for i, ax in enumerate(axs.flatten()):
    nightshade = Nightshade(date, refraction=  i)
    # Extract feature geom in the native projection
    ns_geom = next(nightshade.geometries())
    # Project geometry in axe crs
    proj_ns_geom = ax.projection.project_geometry(ns_geom, nightshade.crs)
    print(proj_ns_geom.area)
    # Create feature from projected geometry
    proj_ns_feat = ShapelyFeature(proj_ns_geom, crs=proj)

    ax.add_feature(proj_ns_feat)
    ax.coastlines()
    ax.set_global()
    ax.set_title(f"Refraction={i}")

which gives the same result as in my initial post

The output from the print(proj_ns_geom.area) line confirms that the issue occurs in the project_geometry function since there are very significant reductions/increases in proj_ns_geom area before the creation of the ShapelyFeature (the expected behavior would be for the evolution of surface area to be a monotonic function)

5151773704494255.0
5246243033890791.0
5344798516527959.0
5448653067659229.0
5559983616304591.0
5684227535698096.0
2013242591740023.2 <- should be greater than the previous area
1872608307796799.8
1744912668814188.0
1628637143988004.0
6325209799254183.0
1425090085522831.8

Looking at the Projection class (#L954) it seems indeed that it may be the origin of the issue according to this comment

        # Determine orientation of polygon.
        # TODO: Consider checking the internal rings have the opposite
        # orientation to the external rings?
        if src_crs.is_geodetic():
            is_ccw = True
        else:
            is_ccw = polygon.exterior.is_ccw

[rcomer]

I have done a little more digging...

The decision whether we have the interior or exterior of a polygon is based on the Shapely is_ccw property.

cartopy/lib/cartopy/crs.py

Lines 1181 to 1185 in d5a07ab

	for ring in rings:
	if ring.is_ccw != is_ccw:
	interior_rings.append(ring)
	else:
	exterior_rings.append(ring)

I intercepted the relevant ring for the "Refraction = 8" case. The is_ccw property for this gives me False, however if I plot it with colours from purple to yellow it definitely looks counter clockwise.

From the Shapely docs, it appears is_ccw is only reliable if the ring is valid. Checking through, we do not have a valid ring regardless of refraction. I was hoping we might be able to use the Shapely make_valid function to work around this. That is not implemented for LinearRing but is for Polygon. So it would mean doing something like

poly = sgeom.Polygon(ring)
geoms = shapely.make_valid(poly)

[find useful polygons within geoms]

new_ring = new_poly.exterior

Unfortunately this process does not preserve the direction of the vertices.

So I still have no idea for a fix :(

[greglucas]

I'm actually somewhat surprised the positive refraction angles give valid geometries at all. We are creating initial geometries with -98 -> +98 latitudes in that situation.

cartopy/lib/cartopy/feature/nightshade.py

Lines 73 to 75 in d5a07ab

	# Fill latitudes up and then down
	y[:npts] = np.linspace(-(90 + refraction), 90 + refraction, npts)
	y[npts:] = y[:npts][::-1]

I'm not entirely clear whether we would want to do some validation on the refraction angle (raise if it isn't negative?), or immediately take the negative absolute value for users refraction = -abs(refraction).

[Eliam76]

I'm actually somewhat surprised the positive refraction angles give valid geometries at all. We are creating initial geometries with -98 -> +98 latitudes in that situation.

cartopy/lib/cartopy/feature/nightshade.py

Lines 73 to 75 in d5a07ab

	# Fill latitudes up and then down
	y[:npts] = np.linspace(-(90 + refraction), 90 + refraction, npts)
	y[npts:] = y[:npts][::-1]

I'm not entirely clear whether we would want to do some validation on the refraction angle (raise if it isn't negative?), or immediately take the negative absolute value for users refraction = -abs(refraction).

Wow I feel bad missing this obvious issue. Since the feature boundary is plane-symmetric it should in theory handle the negative and positive refraction with the following small change

        y[:npts] = np.linspace(-(90 - np.abs(refraction)), 90 - np.abs(refraction), npts)

With this modification there's no sudden switch of inside/outside definition, however there are weird "reflections" of the night shade along the vertical axis

I tried to plot the nightshade geometry points in rotated pole projection (which is used to create the geometry) along with the stereographic projection to understand the origin of these secondary circles.

I'm not sure but it seems that the projection algorithm interpolates the points in the source projection which adds points (not visible on my plot) zipping through the entire map near the south pole, which creates a second circle in stereographic projection.
However, what I don't understand is that I tried to re-derive the sunrise equation to be sure that angles definitions are correct and obtained what I thought to be a slightly different equation for arccos_tmp :

        arccos_tmp = np.clip(np.cos(np.deg2rad(90. - refraction)) /
                             np.cos(np.deg2rad(self.y)), -1, 1)

instead of

        arccos_tmp = np.clip(np.sin(np.deg2rad(refraction)) /
                             np.cos(np.deg2rad(self.y)), -1, 1)

and indeed it seems to work

But I quickly realized that it's in fact the same equation since cos(90° - a) = sin(a). So I don't quite understand why it works with the cosine version but not the sine.

[greglucas]

In the code you updated only the initial latitude range with the absolute value. If you also update the sin(-abs(refraction)) then I think it will be as expected. So, I'd suggest setting refraction = -abs(refraction) earlier in the code so we are referencing the same angle everywhere.

[Eliam76]

In the code you updated only the initial latitude range with the absolute value. If you also update the sin(-abs(refraction)) then I think it will be as expected. So, I'd suggest setting refraction = -abs(refraction) earlier in the code so we are referencing the same angle everywhere.

Yes, it is intentional. The solar terminator is a great circle with a plane pointing toward the Sun. In the Nightshade class the geometry is created in rotated pole coordinates with central longitude and latitude pointing toward the Sun, so that when the refraction is zero, the terminator goes from (rotated) pole to pole. Two small circles resulting from shifting the intersection plane with two opposite refraction value will share the same latitude range (hence the abs) but not the same longitude. And indeed, geometrically, the longitude of the intersection depends on the sine of the refraction angle. Using the absolute value would result in the same intersection geometry with positive and negative values, which is not (geometrically) the expected result.

Anyway I've discovered that the Geodesic class has a circle method that allows you to draw small and large circles without having to resort to the Nightshade class (which isn't its intended role), and therefore doesn't need to use positive refraction values that make no sense for a solar terminator. However, I would still suggest one of the following modification for the Nightshade class :

• using the negative of the absolute value everywhere since positive values have no physical sense concerning solar terminator,
• using the absolute value for the latitude but not the longitude, which can take any refraction angle. However there's still the issue of the interpolation in the rotated pole coordinates, the polygon should be closed passing through the projection boundary for positive refraction (see image below). One solution could be to define the rotated pole projection with opposite central lon/lat for positive refraction angles (and reverse the polygon), but it might be a bit tedious for nothing
• or, more interesting, using the circle method of the Geodesic class, which will (in theory) always produce a valid geometry with regularly spaced points and can be used with positive and negative refraction angles given the correct radius parameter

illustration : the polygon should close passing through the boundary of the projection instead of crossing the whole map as shown by the red line crossed with indigo

[greglucas]

@Eliam76 those are all good and interesting ideas. I'd be happy to review a PR for any of those updates. The geodesic version is especially interesting if it can be considered equivalent (I haven't thought enough about the math and the coordinate systems involved), but it sounds promising, nice thought!