diff --git a/coast/_utils/general_utils.py b/coast/_utils/general_utils.py index 7b40e974..2a590600 100644 --- a/coast/_utils/general_utils.py +++ b/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 diff --git a/coast/data/coast.py b/coast/data/coast.py index 737717ab..4db4627c 100644 --- a/coast/data/coast.py +++ b/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 diff --git a/coast/data/gridded.py b/coast/data/gridded.py index 5d58a3bd..687a644e 100644 --- a/coast/data/gridded.py +++ b/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: diff --git a/config/example_cmip5_grid_t.json b/config/example_cmip5_grid_t.json new file mode 100644 index 00000000..ced87225 --- /dev/null +++ b/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": [] +} diff --git a/example_scripts/find_models_nearest_point.py b/example_scripts/find_models_nearest_point.py new file mode 100644 index 00000000..f8ec52e4 --- /dev/null +++ b/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}") diff --git a/example_scripts/subsetting_by_lon_lat.py b/example_scripts/subsetting_by_lon_lat.py new file mode 100644 index 00000000..a2f7ca45 --- /dev/null +++ b/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])