diff --git a/examples/scalar_data/extend_longitude_over_360.py b/examples/scalar_data/extend_longitude_over_360.py new file mode 100644 index 000000000..04f012bc7 --- /dev/null +++ b/examples/scalar_data/extend_longitude_over_360.py @@ -0,0 +1,135 @@ +""" +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], + mode="wrap")[:,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") diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py index 8a6dd4ccf..093c498ce 100644 --- a/lib/cartopy/crs.py +++ b/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])) diff --git a/lib/cartopy/feature/__init__.py b/lib/cartopy/feature/__init__.py index f021b6b16..3e18d6b7d 100644 --- a/lib/cartopy/feature/__init__.py +++ b/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): """ diff --git a/lib/cartopy/feature/nightshade.py b/lib/cartopy/feature/nightshade.py index 49c3c44f0..55076eac9 100644 --- a/lib/cartopy/feature/nightshade.py +++ b/lib/cartopy/feature/nightshade.py @@ -14,7 +14,7 @@ class Nightshade(ShapelyFeature): def __init__(self, date=None, delta=0.1, refraction=-0.83, - color="k", alpha=0.5, **kwargs): + color="k", alpha=0.5, over=False, **kwargs): """ Shade the darkside of the Earth, accounting for refraction. @@ -94,8 +94,37 @@ def __init__(self, date=None, delta=0.1, refraction=-0.83, kwargs.setdefault('alpha', alpha) geom = sgeom.Polygon(np.column_stack((x, y))) - return super().__init__( - [geom], rotated_pole, **kwargs) + + if over: + # Assume a cylindrical projection for this case + # and perform reprojection here + projected_geom = ccrs.PlateCarree(over=False).project_geometry( + geom, rotated_pole) + + pxs = [] + pys = [] + # Once projected, the result is a MultiPolygon + for item in range(0, len(projected_geom.geoms)): + pxl, pyl = np.array( + projected_geom.geoms._get_geom_item(item).boundary.coords.xy) + pxs.append(pxl) + pys.append(pyl) + + px = np.concatenate(pxs) + py = np.concatenate(pys) + geoms = [] + # We don't know the axes extent, so just go bigger than max lon. + # It will get cropped to the required extent in FeatureArtist. + for offset in np.arange(-720, 720+1, 360): + geoms.append(sgeom.Polygon(np.column_stack((px + offset, py)))) + + geom = sgeom.MultiPolygon(geoms) + + return super().__init__( + [geom], ccrs.PlateCarree(over=True), **kwargs) + else: + return super().__init__( + [geom], rotated_pole, **kwargs) def _julian_day(date): diff --git a/lib/cartopy/img_transform.py b/lib/cartopy/img_transform.py index d0e44b913..cc7e3bfb1 100644 --- a/lib/cartopy/img_transform.py +++ b/lib/cartopy/img_transform.py @@ -202,7 +202,7 @@ def _determine_bounds(x_coords, y_coords, source_cs): bounds = dict(x=[]) half_px = abs(np.diff(x_coords[:2])).max() / 2. - if (((hasattr(source_cs, 'is_geodetic') and + if source_cs.over is False and (((hasattr(source_cs, 'is_geodetic') and source_cs.is_geodetic()) or isinstance(source_cs, ccrs.PlateCarree)) and x_coords.max() > 180): if x_coords.min() < 180: diff --git a/lib/cartopy/mpl/geoaxes.py b/lib/cartopy/mpl/geoaxes.py index 86e08ad67..09d070deb 100644 --- a/lib/cartopy/mpl/geoaxes.py +++ b/lib/cartopy/mpl/geoaxes.py @@ -607,8 +607,24 @@ def coastlines(self, resolution='auto', color='black', **kwargs): """ kwargs['edgecolor'] = color kwargs['facecolor'] = 'none' - feature = cartopy.feature.COASTLINE + if self.projection.over is True: + extent = self.get_extent() + target_proj = ccrs.PlateCarree(over=True) + source_proj = self.projection + extent = self.get_extent() + lon0, lat0 = target_proj.transform_point(extent[0], + extent[2], source_proj) + lon1, lat1 = target_proj.transform_point(extent[1], + extent[3], source_proj) + extent = [lon0, lon1, lat0, lat1] + + feature = cartopy.feature.NaturalEarthFeature_ext( + 'physical', 'coastline', cartopy.feature.auto_scaler, extent=extent, + edgecolor='black', facecolor='never') + """Automatically scaled coastline, including major islands.""" + else: + feature = cartopy.feature.COASTLINE # The coastline feature is automatically scaled by default, but for # anything else, including custom scaler instances, create a new # feature which derives from the default one. @@ -617,6 +633,7 @@ def coastlines(self, resolution='auto', color='black', **kwargs): return self.add_feature(feature, **kwargs) + def tissot(self, rad_km=500, lons=None, lats=None, n_samples=80, **kwargs): """ Add Tissot's indicatrices to the axes. @@ -663,14 +680,28 @@ def tissot(self, rad_km=500, lons=None, lats=None, n_samples=80, **kwargs): if lons.shape != lats.shape: raise ValueError('lons and lats must have the same shape.') - for lon, lat in zip(lons, lats): - circle = geod.circle(lon, lat, rad_km * 1e3, n_samples=n_samples) - geoms.append(sgeom.Polygon(circle)) + if self.projection.over is True: + for lon, lat in zip(lons, lats): + # As this only applies to cylindrical projections, + # the circle will be the same no matter + # the longitude, so create every circle at 180 for convenience. + circle = geod.circle(180, lat, rad_km * 1e3, n_samples=n_samples) + circle[:,0] = (circle[:,0] % 360) -180 + lon + geoms.append(sgeom.Polygon(circle)) + + feature = cartopy.feature.ShapelyFeature( + geoms, ccrs.PlateCarree(over=True), **kwargs) - feature = cartopy.feature.ShapelyFeature(geoms, ccrs.Geodetic(), - **kwargs) + else: + for lon, lat in zip(lons, lats): + circle = geod.circle(lon, lat, rad_km * 1e3, n_samples=n_samples) + geoms.append(sgeom.Polygon(circle)) + + feature = cartopy.feature.ShapelyFeature(geoms, ccrs.Geodetic(), + **kwargs) return self.add_feature(feature) + def add_feature(self, feature, *, autolim=False, **kwargs): """ Add the given :class:`~cartopy.feature.Feature` instance to the axes. @@ -911,9 +942,9 @@ def set_xticks(self, ticks, minor=False, crs=None): if crs is not None and crs != self.projection: if not isinstance(crs, (ccrs._RectangularProjection, ccrs.Mercator)) or \ - not isinstance(self.projection, - (ccrs._RectangularProjection, - ccrs.Mercator)): + not isinstance(self.projection, + (ccrs._RectangularProjection, + ccrs.Mercator)): raise RuntimeError('Cannot handle non-rectangular coordinate ' 'systems.') proj_xyz = self.projection.transform_points(crs, @@ -958,9 +989,9 @@ def set_yticks(self, ticks, minor=False, crs=None): if crs is not None and crs != self.projection: if not isinstance(crs, (ccrs._RectangularProjection, ccrs.Mercator)) or \ - not isinstance(self.projection, - (ccrs._RectangularProjection, - ccrs.Mercator)): + not isinstance(self.projection, + (ccrs._RectangularProjection, + ccrs.Mercator)): raise RuntimeError('Cannot handle non-rectangular coordinate ' 'systems.') proj_xyz = self.projection.transform_points(crs, @@ -985,13 +1016,48 @@ def stock_img(self, name='ne_shaded', **kwargs): """ if name == 'ne_shaded': - source_proj = ccrs.PlateCarree() + target_proj = ccrs.PlateCarree(over=True) + source_proj = self.projection + extent = self.get_extent() + lon0, lat0 = target_proj.transform_point( + extent[0], extent[2], source_proj) + lon1, lat1 = target_proj.transform_point( + extent[1],extent[3], source_proj) + fname = (config["repo_data_dir"] / 'raster' / 'natural_earth' / '50-natural-earth-1-downsampled.png') + image = imread(fname) + + if lon1 - lon0 > 360: + factor = image.shape[1]/360 + negext = np.min([lon0 - -180, 0]) + posext = np.max([0, lon1 - 180]) + + new_image = image + + for offset in np.arange(-360, negext, -360): + new_image = np.concatenate( + (image, new_image), axis=1 + ) + for offset in np.arange(360, posext, 360): + new_image = np.concatenate( + (new_image, image), axis=1 + ) + + leftmost = image[:,int((360+np.mod(negext, -360))*factor):,:] + rightmost = image[:,:int(np.mod(posext, 360)*factor),:] + new_image = np.concatenate((leftmost, new_image, rightmost), + axis=1) + # Only longitudinal extent is different because we have + # modified the image only on this axis + extent = [lon0, lon1, -90, 90] + else: + new_image = image + extent = [-180, 180, -90, 90] - return self.imshow(imread(fname), origin='upper', - transform=source_proj, - extent=[-180, 180, -90, 90], + return self.imshow(new_image, origin='upper', + transform=target_proj, + extent=extent, **kwargs) else: raise ValueError(f'Unknown stock image {name!r}.') @@ -1035,6 +1101,14 @@ def background_img(self, name='ne_shaded', resolution='low', extent=None, extent is used. """ + over = False + try: + lonmin, lonmax, latmin, latmax = extent + if lonmax - lonmin > 360: + over = True + except Exception: + pass + # read in the user's background image directory: if len(_USER_BG_IMGS) == 0: self.read_user_background_images() @@ -1066,7 +1140,7 @@ def background_img(self, name='ne_shaded', resolution='low', extent=None, # now get the projection from the metadata: if _USER_BG_IMGS[name]['__projection__'] == 'PlateCarree': # currently only PlateCarree is defined: - source_proj = ccrs.PlateCarree() + source_proj = ccrs.PlateCarree(over=over) else: raise NotImplementedError('Background image projection undefined') diff --git a/lib/cartopy/tests/crs/test_mercator.py b/lib/cartopy/tests/crs/test_mercator.py index 54dcaacd6..9f1c4441d 100644 --- a/lib/cartopy/tests/crs/test_mercator.py +++ b/lib/cartopy/tests/crs/test_mercator.py @@ -92,3 +92,12 @@ def test_scale_factor(): assert_almost_equal(crs.boundary.bounds, [-18808021, -14585266, 18808021, 17653216], decimal=0) + + +def test_over(): + crs = ccrs.Mercator(over=True) + + other_args = {'ellps=WGS84', 'lon_0=0.0', 'x_0=0.0', 'y_0=0.0', 'units=m', 'over'} + check_proj_params('merc', crs, other_args) + assert_almost_equal(crs.boundary.bounds, + [-63780502, -15496571, 63780502, 18764656], decimal=0) diff --git a/lib/cartopy/util.py b/lib/cartopy/util.py index 28142de35..5277f849a 100644 --- a/lib/cartopy/util.py +++ b/lib/cartopy/util.py @@ -403,3 +403,146 @@ def add_cyclic(data, x=None, y=None, axis=-1, raise ValueError(estr) out_y = _add_cyclic_data(y, axis=yaxis) return out_data, out_x, out_y + + +def extend_lons_xr(dsx, lonmin, lonmax, lon_name="longitude"): + """ + Function to extend data to longitudes beyond 360 degrees in either + direction by duplicating data. + + Version for xarray. Use extend_lons_np for numpy arrays. + + Parameters + ---------- + dsx + an xarray dataset + lonmin + the lowest longitude desired in the extended xarray + lonmax + the highest longitude desired in the extended xarray + lon_name + the xarray coordinate name for longitude + + """ + import xarray as xr + + # Check reasonableness + if lonmin >= lonmax: + print( + "lonmin must be less than lonmax, returning" + " original dataset unchanched" + ) + return dsx + + first_lon = dsx[lon_name][0] + last_lon = dsx[lon_name][-1] + extent = last_lon - first_lon + # Reset to a 360 degree extent before proceeding, otherwise + # we end up extending again the already extended array. + if extent >= 360: + last_lon = first_lon + 360 + dsx = dsx.where( + (dsx[lon_name] < last_lon) & (dsx[lon_name] >= first_lon), + drop=True) + + # Create the name dsx_new to avoid an error if + # neither condition is met. + dsx_new = dsx + if lonmax > last_lon: + nshifts = int((lonmax - last_lon)/360) + 2 + if nshifts !=0: + for n in range(1, abs(nshifts)): + shift = np.sign(nshifts)*n*360 + dsx_new = xr.concat( + [dsx_new, + dsx.assign_coords({lon_name: dsx[lon_name]+shift})], + dim=lon_name) + dsx_new = dsx_new.where(dsx_new[lon_name] < lonmax, drop=True) + + if lonmin < first_lon: + nshifts = int((lonmin - first_lon)/360) - 2 + if nshifts !=0: + for n in range(1, abs(nshifts)): + shift = np.sign(nshifts)*n*360 + dsx_new = xr.concat( + [dsx.assign_coords({lon_name: dsx[lon_name]+shift}), + dsx_new], dim=lon_name) + dsx_new = dsx_new.where(dsx_new[lon_name] > lonmin, drop=True) + + return dsx_new + + +def extend_lons_np(lons, lats, arr, lonmin, lonmax): + """ + Function to extend data to longitudes beyond 360 degrees in either + direction by duplicating data. + + Version for numpy arrays. Use extend_lons_xr for xarrays. + + Parameters + ---------- + lons + a 1-d Numpy array of longitudes + lats + a 1-d Numpy array of latitudes + arr + a Numpy array of the variable to extend with dimension + (..., len(lons), len(lats)) + lonmin + the lowest longitude desired in the extended ndarray + lonmax + the highest longitude desired in the extended ndarray + + """ + # Check reasonableness + if lonmin >= lonmax: + print( + "lonmin must be less than lonmax, returning" + " original dataset unchanched" + ) + return lons, lats, arr + + first_lon = lons[0] + last_lon = lons[-1] + extent = last_lon - first_lon + # Reset to a 360 degree extent before proceeding, otherwise + # we end up extending again the already extended array. + if extent >= 360: + last_lon = first_lon + 360 + valid_lons = (lons < last_lon) & (lons >= first_lon) + arr = arr[..., valid_lons] + lons = lons[valid_lons] + + new_arr = arr + new_lons = lons + if lonmax > last_lon: + nshifts = int((lonmax - last_lon)/360) + 2 + if nshifts !=0: + for n in range(1, abs(nshifts)): + shift = np.sign(nshifts)*n*360 + new_lons = np.concatenate( + (new_lons, lons + shift)) + new_arr = np.concatenate( + (new_arr, arr), + axis=-1 + ) + valid_lons = new_lons < lonmax + new_lons = new_lons[valid_lons] + new_arr = new_arr[...,valid_lons] + + if lonmin < first_lon: + nshifts = int((lonmin - first_lon)/360) - 2 + if nshifts !=0: + for n in range(1, abs(nshifts)): + shift = np.sign(nshifts)*n*360 + new_lons = np.concatenate(( + lons + shift, new_lons)) + new_arr = np.concatenate( + (arr, new_arr), + axis=-1 + ) + valid_lons = new_lons > lonmin + new_lons = new_lons[valid_lons] + new_arr = new_arr[..., valid_lons] + + return new_lons, lats, new_arr