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