Extending longitude to beyond 360 degrees.

examples/scalar_data/extend_longitude_over_360.py

@@ -0,0 +1,133 @@
+"""
+Extending longitude beyond 360 degrees
+======================================
+
+This example demonstrates the use of the over parameter with
+Cartopy.  When using a cylindrical projection, setting over to True
+engages proj's "+over" switch and enables longitudes to be extended
+so that a single map can show more than 360 degrees of the globe.  This
+can be useful for ensuring that large structures can be shown in their
+entirety, unbroken by the edge of the map.
+
+The underlying data needs to be explicitly extended, for which the
+utility routines extend_lons_np and extend_lons_xr are included in
+cartopy.util.  This allows for the possibility of non-repeated data,
+such as an object's trajectory; a trivial example of this is shown
+by plotting Chicago twice with different coloured points.
+
+The script also applies Nightshade, which requires the "over"
+parameter to be set, and tissot and coastlines, which work transparently.
+
+Note that due a limitation in the underlying proj library, the longitudes
+are limited to [-572.95, 572.95].
+"""
+
+import datetime
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.ndimage import gaussian_filter
+
+import cartopy.crs as ccrs
+from cartopy.feature.nightshade import Nightshade
+from cartopy.util import extend_lons_np
+
+
+date = datetime.datetime.now()
+#date = datetime.datetime(2025, 3, 20, 9, 1, 0) # NH Spring equinox 2025
+#date = datetime.datetime(2025, 6, 21, 2, 42, 0) # NH Summer solstice 2025
+
+# Coordinates of a few places
+toulouse = 43.604500, 1.444000
+nyalesund = 78.9, 11.9
+chicago = 39.162, -84.45689
+
+### Load some Copernicus ocean altimetry data
+# ds = open_dataset("nrt_global_allsat_phy_l4_20241125_20241125.nc")
+# u = ds.ugosa[0,2:-1:4,1::3]
+# v = ds.vgosa[0,2:-1:4,1::3]
+# u, v = ds.ugosa[0, ...] , ds.vgosa[0,...]
+# speed = np.sqrt(u**2 + v**2)
+# x, y = speed.longitude, speed.latitude
+
+### Or just make up some random, pretend meteorological data
+incr = 0.25
+nlats, nlons = 720, 1440
+lat = np.arange(incr/2, 90, incr)
+lat = np.concat([-lat[-1::-1], lat])
+lon = np.arange(incr/2, 180, incr)
+lon = np.concat([-lon[-1::-1], lon])
+size = (nlats, nlons)
+rfield = np.random.normal(size=size)
+rfield = np.concat((rfield, rfield), axis=1)
+feature = gaussian_filter(rfield,
+            sigma=[nlats/12,nlons/4])[:,int(nlons/2):int(nlons*3/2)]
+
+# var = speed
+var = feature
+
+# Extrema:
+factor = 1
+vmin = np.nanmin(var)/factor
+vmax = np.nanmax(var)/factor
+# # Make symmetrical
+# vmin = -vmax
+
+extend_cbar = "both"
+
+def a_transform(arr):
+    """Some random transformation to differentiate data."""
+    amax = arr.max()
+    return amax/2 - arr
+
+# Design a few test configurations
+mapconf1 = dict(title="Standard Mercator", lonmin=-180, lonmax=180, over=False,
+                trans=ccrs.PlateCarree, proj=ccrs.Mercator,
+                central_longitude=0)
+mapconf2 = dict(title="Extended Mercator", lonmin=-390, lonmax=525, over=True,
+                trans=ccrs.PlateCarree, proj=ccrs.Mercator,
+                central_longitude=0)
+mapconf3 = dict(title="Extended Plate Carrée", lonmin=-390, lonmax=525, over=True,
+                trans=ccrs.PlateCarree, proj=ccrs.PlateCarree,
+                central_longitude=0)
+
+mapconfs = [mapconf1, mapconf2, mapconf3]
+
+for ind, mapconf in enumerate(mapconfs[1:2]):
+    print(ind, mapconf)
+    lonmin, lonmax = mapconf["lonmin"], mapconf["lonmax"]
+    over = mapconf["over"]
+    central_longitude = mapconf["central_longitude"]
+    proj = mapconf["proj"]
+    trans = mapconf["trans"]
+    projection = proj(over=over, central_longitude=central_longitude)
+    transform = trans(over=over, central_longitude=central_longitude)
+    title = mapconf["title"] + " [" + str(lonmin) + "," + str(lonmax) + "]"
+    lon_ext, lat_ext, var_ext = extend_lons_np(lon, lat, var, lonmin, lonmax)
+    ## Uncomment the following line to highlight the extra data (for Xarrays)
+    # var_ext = var_ext.where(
+    #     np.abs(var_ext.longitude)<180, a_transform(var_ext))
+
+    fig = plt.figure(figsize=(10,4), facecolor=(0,0,0,0))
+    ax = plt.axes(projection=projection)
+    ax.pcolormesh(lon_ext, lat_ext, var_ext,
+                  vmin=vmin, vmax = vmax, transform=transform)
+    ax.coastlines()
+    ax.gridlines(draw_labels=True, dms=True,
+                 x_inline=False, y_inline=False,
+                 crs=ccrs.PlateCarree(over=over))
+    ax.add_feature(Nightshade(date, alpha=0.2, delta=0.1, over=over))
+    ax.plot(toulouse[1], toulouse[0], "ro", transform=ccrs.Geodetic())
+    ax.plot(nyalesund[1], nyalesund[0], "ro", transform=ccrs.Geodetic())
+    ax.plot(chicago[1], chicago[0], "ro", transform=ccrs.Geodetic())
+    ax.plot(chicago[1]+360, chicago[0], "orange", marker="o", transform=ccrs.Geodetic())
+    ax.stock_img()
+    ax.tissot(
+        lons=np.arange(-590, 580,  200),
+        lats=[-75, -60, -45, -30, -10, 20, 50, 65, 80],
+        n_samples=40
+    )
+    ax.set_title(title)
+
+plt.show()
+#plt.savefig("eke_cartopy_extended_mercator.png")

lib/cartopy/crs.py

@@ -135,7 +135,7 @@ class is the very core of cartopy, all coordinate reference systems in cartopy
     #: Whether this projection can handle ellipses.
     _handles_ellipses = True
 
-    def __init__(self, proj4_params, globe=None):
+    def __init__(self, proj4_params, globe=None, over=False):
         """
         Parameters
         ----------
@@ -149,8 +149,23 @@ def __init__(self, proj4_params, globe=None):
         globe: :class:`~cartopy.crs.Globe` instance, optional
             If omitted, the default Globe instance will be created.
             See :class:`~cartopy.crs.Globe` for details.
+        over: boolean, optional.  Defaults to False.
+            If True, adds the +over parameter to the PROJ string,
+            and enables longitudes to be extended beyond 360
+            degrees.  This only makes sense with cylindrical
+            projections, and enables more than one rotation of
+            the globe to be visualised.
 
         """
+
+        self.over = over
+        if over is True:
+            if isinstance(proj4_params, list):
+                proj4_params.append(("over", None))
+            elif isinstance(proj4_params, str):
+                proj4_params+" +over"
+            else:
+                print("Error: proj4_params neither str nor list")
         self.input = (proj4_params, globe)
 
         # for compatibility with pyproj.CRS and rasterio.crs.CRS
@@ -419,8 +434,9 @@ def transform_points(self, src_crs, x, y, z=None, trap=False):
                     self.is_geodetic()):
                 # convert from [0,360] to [-180,180]
                 x = np.array(x, copy=True)
-                to_180 = (x > 180) | (x < -180)
-                x[to_180] = (((x[to_180] + 180) % 360) - 180)
+                if self.over is False:
+                    to_180 = (x > 180) | (x < -180)
+                    x[to_180] = (((x[to_180] + 180) % 360) - 180)
             try:
                 result[:, 0], result[:, 1], result[:, 2] = \
                     _safe_pj_transform(src_crs, self, x, y, z, trap=trap)
@@ -573,7 +589,7 @@ class Geodetic(CRS):
 
     """
 
-    def __init__(self, globe=None):
+    def __init__(self, globe=None, over=False):
         """
         Parameters
         ----------
@@ -583,7 +599,7 @@ def __init__(self, globe=None):
         """
         proj4_params = [('proj', 'lonlat')]
         globe = globe or Globe(datum='WGS84')
-        super().__init__(proj4_params, globe)
+        super().__init__(proj4_params, globe, over=over)
 
     # XXX Implement fwd such as Basemap's Geod.
     # Would be used in the tissot example.
@@ -699,7 +715,12 @@ def __init__(self, *args, **kwargs):
         elif self.is_geographic:
             # If the projection is geographic without an area of use, assume
             # the bounds are the full globe.
-            self.bounds = (-180, 180, -90, 90)
+            # So that feature_artist can handle +over, the
+            # bounds are extended beyond +/- 180 degrees.
+            if self.over:
+                self.bounds = (-572.95, 572.95, -90, 90)
+            else:
+                self.bounds = (-180, 180, -90, 90)
 
     @property
     def boundary(self):
@@ -1305,10 +1326,10 @@ class _RectangularProjection(Projection, metaclass=ABCMeta):
     """
     _wrappable = True
 
-    def __init__(self, proj4_params, half_width, half_height, globe=None):
+    def __init__(self, proj4_params, half_width, half_height, globe=None, over=False):
         self._half_width = half_width
         self._half_height = half_height
-        super().__init__(proj4_params, globe=globe)
+        super().__init__(proj4_params, globe=globe, over=over)
 
     @property
     def boundary(self):
@@ -1348,17 +1369,20 @@ def _ellipse_boundary(semimajor=2, semiminor=1, easting=0, northing=0, n=201):
 
 
 class PlateCarree(_CylindricalProjection):
-    def __init__(self, central_longitude=0.0, globe=None):
+    def __init__(self, central_longitude=0.0, globe=None, over=False):
         globe = globe or Globe(semimajor_axis=WGS84_SEMIMAJOR_AXIS)
         proj4_params = [('proj', 'eqc'), ('lon_0', central_longitude),
                         ('to_meter', math.radians(1) * (
                             globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS)),
                         ('vto_meter', 1)]
-        x_max = 180
+        if over is True:
+            x_max = 572.95 # Maximum allowed value in pyproj with +over
+        else:
+            x_max = 180
         y_max = 90
         # Set the threshold around 0.5 if the x max is 180.
-        self.threshold = x_max / 360
-        super().__init__(proj4_params, x_max, y_max, globe=globe)
+        self.threshold = 0.5
+        super().__init__(proj4_params, x_max, y_max, globe=globe, over=over)
 
     def _bbox_and_offset(self, other_plate_carree):
         """
@@ -1628,7 +1652,8 @@ class Mercator(Projection):
     def __init__(self, central_longitude=0.0,
                  min_latitude=-80.0, max_latitude=84.0,
                  globe=None, latitude_true_scale=None,
-                 false_easting=0.0, false_northing=0.0, scale_factor=None):
+                 false_easting=0.0, false_northing=0.0, scale_factor=None,
+                 over=False):
         """
         Parameters
         ----------
@@ -1650,6 +1675,8 @@ def __init__(self, central_longitude=0.0,
             Y offset from the planar origin in metres. Defaults to 0.
         scale_factor: optional
             Scale factor at natural origin. Defaults to unused.
+        over: optional
+            Extend map beyond 360 degrees. Defaults to False.
 
         Notes
         -----
@@ -1674,13 +1701,16 @@ def __init__(self, central_longitude=0.0,
             else:
                 proj4_params.append(('k_0', scale_factor))
 
-        super().__init__(proj4_params, globe=globe)
+        super().__init__(proj4_params, globe=globe, over=over)
 
         # Need to have x/y limits defined for the initial hash which
         # gets used within transform_points for caching
         self._x_limits = self._y_limits = None
         # Calculate limits.
-        minlon, maxlon = self._determine_longitude_bounds(central_longitude)
+        if over is False:
+            minlon, maxlon = self._determine_longitude_bounds(central_longitude)
+        else:
+            minlon, maxlon = -572.95, 572.95
         limits = self.transform_points(self.as_geodetic(),
                                        np.array([minlon, maxlon]),
                                        np.array([min_latitude, max_latitude]))

lib/cartopy/feature/__init__.py

@@ -17,6 +17,7 @@
 from abc import ABCMeta, abstractmethod
 
 import numpy as np
+import shapely.affinity as saffinity
 import shapely.geometry as sgeom
 
 import cartopy.crs
@@ -286,6 +287,7 @@ def geometries(self):
         Returns an iterator of (shapely) geometries for this feature.
 
         """
+
         key = (self.name, self.category, self.scale)
         if key not in _NATURAL_EARTH_GEOM_CACHE:
             path = shapereader.natural_earth(resolution=self.scale,
@@ -326,6 +328,91 @@ def with_scale(self, new_scale):
         return NaturalEarthFeature(self.category, self.name, new_scale,
                                    **self.kwargs)
 
+class NaturalEarthFeature_ext(NaturalEarthFeature):
+    def __init__(self, category, name, scale, extent=None, **kwargs):
+        """
+        Parameters
+        ----------
+        category
+            The category of the dataset, i.e. either 'cultural' or 'physical'.
+        name
+            The name of the dataset, e.g. 'admin_0_boundary_lines_land'.
+        scale
+            The dataset scale, i.e. one of '10m', '50m', or '110m',
+            or Scaler object. Dataset scales correspond to 1:10,000,000,
+            1:50,000,000, and 1:110,000,000 respectively.
+
+        Other Parameters
+        ----------------
+        **kwargs
+            Keyword arguments to be used when drawing this feature.
+
+        """
+
+        # Set over to True for all cases /maltron
+        super(NaturalEarthFeature, self).__init__(
+            cartopy.crs.PlateCarree(over=True), **kwargs
+        )
+        self.category = category
+        self.name = name
+
+        # Cast the given scale to a (constant) Scaler if a string is passed.
+        if isinstance(scale, str):
+            scale = Scaler(scale)
+
+        self.scaler = scale
+        # Make sure this is a valid resolution
+        self._validate_scale()
+        self.extent = extent
+
+    def geometries(self):
+        """
+        Returns an iterator of (shapely) geometries for this feature.
+
+        """
+        key = (self.name, self.category, self.scale)
+        if key not in _NATURAL_EARTH_GEOM_CACHE:
+            path = shapereader.natural_earth(resolution=self.scale,
+                                             category=self.category,
+                                             name=self.name)
+            reader = shapereader.Reader(path)
+            if reader.crs is not None:
+                self._crs = reader.crs
+            geometries = tuple(reader.geometries())
+            _NATURAL_EARTH_GEOM_CACHE[key] = geometries
+        else:
+            geometries = _NATURAL_EARTH_GEOM_CACHE[key]
+
+        if self.extent is not None:
+            if self.extent[1] - self.extent[0] > 360:
+                def extend_geoms(geoms, extent, xoffset=360):
+                    extent_geom = sgeom.box(extent[0], extent[2],
+                                            extent[1], extent[3])
+                    new_geoms = []
+                    for geom in geoms:
+                        geom_ext = saffinity.translate(geom, xoff=xoffset, yoff=0)
+                        if extent_geom.intersects(geom_ext):
+                            new_geoms.append(geom_ext)
+                    return new_geoms
+
+                geoms_left = []
+                geoms_right = []
+                if self.extent[1]-self.extent[0] > 360:
+                    if self.extent[0] < -180:
+                        for offset in np.arange(-360, self.extent[0]-360, -360):
+                            geoms_left += extend_geoms(
+                                geometries, self.extent, xoffset=offset
+                            )
+
+                    if self.extent[1] > 180:
+                        for offset in np.arange(360, self.extent[1]+360, 360):
+                            geoms_right += extend_geoms(
+                                geometries, self.extent, xoffset=offset
+                            )
+
+                geometries = tuple(geoms_left) + geometries + tuple(geoms_right)
+
+        return iter(geometries)
 
 class GSHHSFeature(Feature):
     """