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).