updating Steve's changes to match with users updated requirement

coast/_utils/general_utils.py

@@ -170,10 +170,47 @@ def polar_to_cartesian(r, theta, degrees=True):
     return x, y
 
 
-def subset_indices_lonlat_box(array_lon, array_lat, lon_min, lon_max, lat_min, lat_max):
-    ind_lon = np.logical_and(array_lon <= lon_max, array_lon >= lon_min)
-    ind_lat = np.logical_and(array_lat <= lat_max, array_lat >= lat_min)
-    ind = np.where(np.logical_and(ind_lon, ind_lat))
+def subset_indices_lonlat_box(array_lon, array_lat, ww=np.NaN, ee=np.NaN, ss=np.NaN, nn=np.NaN):
+    """array_lon  longitudes (degrees)  - no particular order
+      array_lat  latitudes  (degrees)  - no particular order
+      ww   western boundary (degrees)  - identifies limiting meridian on the left (can be +/-'ve)
+                                         if NaN ee must also be NaN and no selection takes place
+                                         on longitude
+      ee   eastern boundary (degrees)  - identifies limiting meridian on the right (can be +/-'ve)
+      ss   southern boundary (degrees) - can be NaN (i.e. south pole)
+      nn   northern boundary (degrees) - can be NaN (i.e. north pole)
+
+    Note (1) Boundaries can be signed. E.g. -150 = 150W = 210E
+    Note (2) Vertices on Western and Southern boundaries are selected. Those
+             on Northern and Eastern boundaries are not. E.g. ww=1,ee=2
+             might select a single meridian of vertices if grid interval = 1
+    Note (3) Models sometimes duplicate vertices at the (model) boundary meridians.
+             These duplications are not removed by this function
+    """
+    if len(array_lat) != len(array_lon):
+        raise ValueError("The arrays must match on dimension")
+    if np.isnan(ss):
+        ss = -100
+    if np.isnan(nn):
+        nn = 100
+    if ss >= nn:
+        raise ValueError("Northern limit must be greater than Southern")
+    ind_lat = np.logical_and(array_lat >= ss, array_lat < nn)
+    if np.logical_xor(np.isnan(ww), np.isnan(ee)):
+        raise ValueError("both E/W boundaries must be non-NaN or omitted altogether")
+
+    if np.isnan(ww):
+        return ind_lat
+    #
+    #  Meridian limits specified. We reset the zero meridian to co-incide with ww
+    #
+    ww = np.mod(ww, 360)
+    ee = np.mod(ee, 360)
+    if ww == ee:
+        raise ValueError("Not allowed equal values for western and eastern boundaries")
+    array_lon = np.mod(array_lon - ww, 360)
+    ind_lon = array_lon < np.mod(ee - ww, 360)
+    ind = np.logical_and(ind_lon, ind_lat)
     return ind
 
 

coast/data/coast.py

@@ -1,6 +1,7 @@
 from dask import array
 import xarray as xr
 import numpy as np
+from .._utils import general_utils as gu
 from dask.distributed import Client
 import copy
 from .._utils.logging_util import get_slug, debug, info, warn, warning
@@ -211,10 +212,10 @@ def subset_indices_by_distance(self, centre_lon: float, centre_lat: float, radiu
         """
         This method returns a `tuple` of indices within the `radius` of the lon/lat point given by the user.
 
-        Distance is calculated as haversine - see `self.calculate_haversine_distance`
+        Distance is calculated as haversine - see `calculate_haversine_distance`
 
-        :param centre_lon: The longitude of the users central point
-        :param centre_lat: The latitude of the users central point
+        :param centre_lon: The longitude of the user's central point
+        :param centre_lat: The latitude of the user's central point
         :param radius: The haversine distance (in km) from the central point
         :return: All indices in a `tuple` with the haversine distance of the central point
         """
@@ -226,65 +227,34 @@ def subset_indices_by_distance(self, centre_lon: float, centre_lat: float, radiu
         # Calculate the distances between every model point and the specified
         # centre. Calls another routine dist_haversine.
 
-        dist = self.calculate_haversine_distance(centre_lon, centre_lat, lon, lat)
+        dist = gu.calculate_haversine_distance(centre_lon, centre_lat, lon, lat)
         indices_bool = dist < radius
         indices = np.where(indices_bool.compute())
 
         return xr.DataArray(indices[0]), xr.DataArray(indices[1])
 
-    def subset_indices_lonlat_box(self, lonbounds, latbounds):
-        """Generates array indices for data which lies in a given lon/lat box.
+    def subset_indices_lonlat_box(self, ww=np.NaN, ee=np.NaN, ss=np.NaN, nn=np.NaN):
+        """
+          ww   western boundary (degrees)  - identifies limiting meridian on the left (can be +/-'ve)
+                                             if NaN ee must also be NaN and no selection takes place
+                                             on longitude
+          ee   eastern boundary (degrees)  - identifies limiting meridian on the right (can be +/-'ve)
+          ss   southern boundary (degrees) - can be NaN (i.e. south pole)
+          nn   northern boundary (degrees) - can be NaN (i.e. north pole)
+        Note (1) boundaries can be signed. E.g. -150 is 150W = 210E
+        Note (2) models sometimes duplicate vertices at boundary meridians.
+                 These duplications are not removed by subsetting
+
 
-        Keyword arguments:
-        lon       -- Longitudes, 1D or 2D.
-        lat       -- Latitudes, 1D or 2D
-        lonbounds -- Array of form [min_longitude=-180, max_longitude=180]
-        latbounds -- Array of form [min_latitude, max_latitude]
 
-        return: Indices corresponding to datapoints inside specified box
         """
         debug(f"Subsetting {get_slug(self)} indices within lon/lat")
+
         lon_str = "longitude"
         lat_str = "latitude"
         lon = self.dataset[lon_str].copy()  # TODO Add a comment explaining why this needs to be copied
         lat = self.dataset[lat_str]
-        ff = lon > lonbounds[0]
-        ff *= lon < lonbounds[1]
-        ff *= lat > latbounds[0]
-        ff *= lat < latbounds[1]
-
-        return np.where(ff)
-
-    @staticmethod
-    def calculate_haversine_distance(lon1, lat1, lon2, lat2):  # TODO This could be a static method
-        """
-        # Estimation of geographical distance using the Haversine function.
-        # Input can be single values or 1D arrays of locations. This
-        # does NOT create a distance matrix but outputs another 1D array.
-        # This works for either location vectors of equal length OR a single loc
-        # and an arbitrary length location vector.
-        #
-        # lon1, lat1 :: Location(s) 1.
-        # lon2, lat2 :: Location(s) 2.
-        """
-
-        debug(f"Calculating haversine distance between {lon1},{lat1} and {lon2},{lat2}")
-
-        # Convert to radians for calculations
-        lon1 = np.deg2rad(lon1)
-        lat1 = np.deg2rad(lat1)
-        lon2 = np.deg2rad(lon2)
-        lat2 = np.deg2rad(lat2)
-
-        # Latitude and longitude differences
-        dlat = (lat2 - lat1) / 2
-        dlon = (lon2 - lon1) / 2
-
-        # Haversine function.
-        distance = np.sin(dlat) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon) ** 2
-        distance = 2 * 6371.007176 * np.arcsin(np.sqrt(distance))
-
-        return distance
+        return gu.subset_indices_lonlat_box(lon, lat, ww=ww, ee=ee, nn=nn, ss=ss)
 
     def get_subset_as_xarray(
         self, var: str, points_x: slice, points_y: slice, line_length: int = None, time_counter: int = 0

coast/data/gridded.py

@@ -269,8 +269,9 @@ def find_j_i(self, lat: float, lon: float):
         :return: the y and x coordinates for the NEMO object's grid_ref, i.e. t,u,v,f,w.
         """
         debug(f"Finding j,i for {lat},{lon} from {get_slug(self)}")
-        dist2 = np.square(self.dataset.latitude - lat) + np.square(self.dataset.longitude - lon)
-        [y, x] = np.unravel_index(dist2.argmin(), dist2.shape)
+        # dist2 = xr.ufuncs.square(self.dataset.latitude - lat) + xr.ufuncs.square(self.dataset.longitude - lon)
+        dist = general_utils.calculate_haversine_distance(self.dataset.longitude, self.dataset.latitude, lon, lat)
+        [y, x] = np.unravel_index(dist.argmin(), dist.shape)
         return [y, x]
 
     def find_j_i_domain(self, lat: float, lon: float, dataset_domain: xr.DataArray):
@@ -287,8 +288,8 @@ def find_j_i_domain(self, lat: float, lon: float, dataset_domain: xr.DataArray):
         debug(f"Finding j,i domain for {lat},{lon} from {get_slug(self)} using {get_slug(dataset_domain)}")
         internal_lat = dataset_domain[self.grid_vars[1]]  # [f"gphi{self.grid_ref.replace('-grid','')}"]
         internal_lon = dataset_domain[self.grid_vars[0]]  # [f"glam{self.grid_ref.replace('-grid','')}"]
-        dist2 = np.square(internal_lat - lat) + np.square(internal_lon - lon)
-        [_, y, x] = np.unravel_index(dist2.argmin(), dist2.shape)
+        dist = general_utils.calculate_haversine_distance(internal_lon, internal_lat, lon, lat)
+        [_, y, x] = np.unravel_index(dist.argmin(), dist.shape)
         return [y, x]
 
     def transect_indices(self, start: tuple, end: tuple) -> tuple:

config/example_cmip5_grid_t.json

@@ -0,0 +1,37 @@
+{
+  "type": "gridded",
+  "dimensionality": 4,
+  "chunks": {},
+  "grid_ref": {
+    "t-grid": []
+  },
+  "dataset": {
+    "dimension_map": {
+      "time": "t_dim",
+      "lev": "z_dim",
+      "i": "x_dim",
+      "j": "y_dim"
+    },
+    "variable_map": {
+      "time":"time",
+      "lat": "latitude",
+      "lon": "longitude",
+      "so": "salinity"
+    },
+    "coord_vars": [
+        "longitude",
+        "latitude",
+        "time",
+        "depth_0"
+        ]
+  },
+  "domain": {
+    "dimension_map": {},
+    "variable_map": {}
+  },
+  "static_variables": {
+    "not_grid_vars" : [],
+    "delete_vars": []
+  },
+  "processing_flags": []
+}

example_scripts/find_models_nearest_point.py

@@ -0,0 +1,22 @@
+"""
+Demonstration of wrapped longitude coordinates when finding the nearest
+ grid index for specified latitude and lonitude
+"""
+import coast
+import os.path as path
+
+# Global data files
+dn_files = "/projectsa/COAsT/GLOBAL_example_data/"
+filename = "so_Omon_IPSL-CM5A-LR_historical_r1i1p1_200001-200512.nc"
+fn_ispl_dat_t = path.join(dn_files, filename)
+# Configuration files describing the data files
+fn_config_t_grid = "./config/example_cmip5_grid_t.json"
+
+# Load the data as a Gridded object
+ispl = coast.Gridded(fn_ispl_dat_t, config=fn_config_t_grid)
+
+[j1, i1] = ispl.find_j_i(-9, 50)
+print(f"At (-9N,50E) nearest j,i indices: {j1,i1}")
+
+[j1, i1] = ispl.find_j_i(-9, -310)  # Same point on globe gives same result
+print(f"At (-9N,-310E) nearest j,i indices: {j1,i1}")

example_scripts/subsetting_by_lon_lat.py

@@ -0,0 +1,26 @@
+import numpy as np
+import coast
+
+lon = np.arange(10, 30)
+checkDim = False
+if checkDim:
+    lat = np.arange(20, 41)  # Testing array size check
+else:
+    lat = lon  # = lon-10 would leave some selections empty in the following
+print("The dataset is a set of vertices on a straight line in lon/lat space")
+print("Here's what we are starting with - for longitude (= latitude)")
+print(lon)
+ind = coast.general_utils.subset_indices_lonlat_box(lon, lat, ss=15, ww=20, ee=25)
+print("Here's the subset for [20,25] - note omission of 25")
+print(lon[ind])
+print(
+    "Interchange east and west. This should go from 25E all the way around to 20E. But for some reason it misses out 10-14"
+)
+ind = coast.general_utils.subset_indices_lonlat_box(lon, lat, ss=15, ww=25, ee=20)
+print(lon[ind])
+print("Use negative longitude")
+ind = coast.general_utils.subset_indices_lonlat_box(lon, lat, ss=15, ww=-340, ee=25)
+print(lon[ind])
+# print("Don't use same value for the two boundaries")
+# ind = cst.general_utils.subset_indices_lonlat_box(lat,lon,ss=-15,ww=25,ee=-335)
+# print(lon[ind])