Merge pull request #303 from SuperDARN/ehn/terminator

EHN: Options to Plot a Terminator

docs/user/coordinates.md

@@ -71,4 +71,56 @@ Spatial plots have two options for projections.
     Some combinations of Projs/Coords/RangeEstimates are not designed to work. 
     For example, you cannot plot a fan plot using range gates, spatial plots require a value in kilometers. 
     At the moment, AACGM Coordinates do not plot on Geographic projections as it has not been developed yet. 
-    Convection maps only support polar projections due to lack of interest in requiring geographic projections.
+    Convection maps only support polar projections due to lack of interest in requiring geographic projections.
+
+# Including a Terminator
+
+Spatial plots have the option to include a terminator called `nightshade` at a given height in the ionosphere. This functions uses the *Cartopy* `nightshade` function.
+Nightshade is only available using the geographic projection and can be implemented by adding `nightshade=250` to the spacial plot call where 250 is the desired height in the
+ionosphere to be in the Earth's shadow. If you would like to plot your own terminator on any plot, the `terminator` function will return the anti-sub-solar position and the 
+great circle distance to the terminator in geographic coordinates:
+```python
+antisolarpsn, arc_length, angle_of_terminator = terminator(date, nightshade)
+```
+The `antisolarpsn` is given in degrees lon, lat. The `arc_length` is in kilometers and the `angle_of_terminator` is the angle from the subsolar point to the terminator (i.e. is 90 degrees at ground level).
+The terminator position can be calculated using `(lat, lon) = new_coordinate(lat, lon, arc_length, bearing, R=Re)` for any bearing from the antisolar position. This can be converted to magnetic coordinates using the
+AACGMv2 library. Unfortunately, matplotlib is unable to plot the terminator using `fill` consistently, hence we leave this option up to the user.
+An example of this is shown below:
+```python
+import pydarn
+import aacgmv2
+import datetime as dt
+import matplotlib.pyplot as plt
+import numpy as np
+
+# North Winter
+_,_,ax1,ccrs1=pydarn.Fan.plot_fov(66, dt.datetime(2023, 12, 21, 0, 0),
+    lowlat= 5, boundary=True, line_color='red', coastline=True, nightshade=250)
+
+# Test to plot terminator if ever required - plot line not fill!
+# Get antisolar point in geographic coords and radius of terminator
+# at given height
+date = dt.datetime(2023, 12, 21, 0, 0)
+antisolarpsn, arc, ang = pydarn.terminator(date, 250)
+# Convert position to magnetic coordinates
+mlat, lon_mag, _ =  aacgmv2.convert_latlon(antisolarpsn[1],
+                                           antisolarpsn[0],
+                                           250, date, method_code='G2A')
+# Shift to MLT
+shifted_mlts = lon_mag - (aacgmv2.convert_mlt(lon_mag, date) * 15)
+shifted_lons = lon_mag - shifted_mlts
+mlon = np.radians(shifted_lons)
+# Get positions at a distance from new position to plot terminator
+lats = []
+lons = []
+for b in range(-180,180,1):
+    (lat, lon) = pydarn.new_coordinate(mlat, shifted_lons, arc, b, R=pydarn.Re)
+    nlon =np.radians(lon)
+    lats.append(lat)
+    lons.append(nlon)
+lats2 = np.zeros(len(lats))
+plt.plot(np.squeeze(lons), np.squeeze(lats), color='b', zorder=2.0,
+         linewidth=3.0, linestyle='dashed')
+
+plt.show()
+```

pydarn/__init__.py

@@ -36,6 +36,7 @@
 from .exceptions.warning_formatting import partial_record_warning
 from .exceptions.warning_formatting import cartopy_warning
 from .exceptions.warning_formatting import cartopy_print_warning
+from .exceptions.warning_formatting import nightshade_warning
 
 # importing utils
 from .utils.constants import Re
@@ -56,6 +57,7 @@
 from .utils.scan import build_scan
 from .utils.geo import geocentric_coordinates
 from .utils.coordinates import Coords
+from .utils.terminator import terminator, new_coordinate
 from .utils.recalculate_elevation import recalculate_elevation
 from .utils.filters import Boxcar
 

pydarn/exceptions/warning_formatting.py

@@ -43,6 +43,19 @@ def cartopy_warning():
                  "If you do not need to use cartopy for your plotting, ignore "
                  "this message.")
 
+def nightshade_warning():
+    """
+    This warning prints when a user tries to use something that is not
+    implemented yet, however the plot will still plot without the
+    feature
+    """
+    warnings.warn("Nightshade is not implemented for plots in AACGMv2 "
+                  "coordinates. You can plot nightshade in geographic "
+                  "coordinates, or you can call the terminator function "
+                  "to return the antisolarposition and radius to "
+                  "decide how you would like to plot."
+                  "https://pydarn.readthedocs.io/en/main/user/coordinates/")
+
 def partial_record_warning():
     """
     prints a warning that the data chosen to be plotted is missing some

pydarn/plotting/projections.py

@@ -6,6 +6,7 @@
 # 2022-03-22 MTS removed coastline call and added grid lines to cartopy plotting
 # 2022-05-20 CJM added options to plot coastlines
 # 2022-06-13 Elliott Day don't create new ax if ax passed in to Projs
+# 2023-02-22 CJM added options for nightshade
 #
 # Disclaimer:
 # pyDARN is under the LGPL v3 license found in the root directory LICENSE.md
@@ -21,15 +22,19 @@
 import aacgmv2
 import enum
 import matplotlib.pyplot as plt
+import matplotlib.patches as patches
 import numpy as np
 from packaging import version
 
-from pydarn import Hemisphere, plot_exceptions
+from pydarn import (Hemisphere, plot_exceptions, terminator, Re,
+                    new_coordinate, nightshade_warning)
 try:
     import cartopy
     # from cartopy.mpl import geoaxes
     import cartopy.crs as ccrs
     import cartopy.feature as cfeature
+    from cartopy.feature.nightshade import Nightshade
+    from shapely.geometry import mapping, Polygon
     cartopyInstalled = True
     if version.parse(cartopy.__version__) < version.parse("0.19"):
         cartopyVersion = False
@@ -74,7 +79,7 @@ def convert_to_mag(coords, date, alt):
 def axis_polar(date, ax: object = None, lowlat: int = 30,
                hemisphere: Hemisphere = Hemisphere.North,
                grid_lines: bool = True, coastline: bool = False,
-               **kwargs):
+               nightshade: int = 0, **kwargs):
 
     """
     Sets up the polar plot matplotlib axis object, for use in
@@ -144,13 +149,18 @@ def axis_polar(date, ax: object = None, lowlat: int = 30,
         # Plot each geometry object
         for geom in cc_mag.geometries():
             plt.plot(*geom.coords.xy, color='k', linewidth=0.5, zorder=2.0)
+
+    if nightshade:
+        nightshade_warning()
+
     return ax, None
 
 
 def axis_geological(date, ax: object = None,
                     hemisphere: Hemisphere = Hemisphere.North,
                     lowlat: int = 30, grid_lines: bool = True,
-                    coastline: bool = False, **kwargs):
+                    coastline: bool = False, nightshade: int = 0,
+                    **kwargs):
     """
     
     Sets up the cartopy orthographic plot axis object, for use in
@@ -209,6 +219,12 @@ def axis_geological(date, ax: object = None,
 
     if coastline is True:
         ax.coastlines()
+
+    if nightshade:
+        refraction_value = -np.degrees(np.arccos(Re / (Re + nightshade)))
+        ns = Nightshade(date, refraction=refraction_value, alpha=0.1)
+        ax.add_feature(ns)
+
     return ax, ccrs