diff --git a/tests/coordinates.py b/tests/coordinates.py new file mode 100644 index 0000000000000000000000000000000000000000..86a3ec6fbb7355da7211c5925d49d4fa96f26f36 --- /dev/null +++ b/tests/coordinates.py @@ -0,0 +1,30 @@ +#!/bin/python3 + +from toargridding.variables import Coordinate +from toargridding.metadata import Coordinates +import numpy as np +from pandas import Series +from toargridding.grids import RegularGrid +import matplotlib.pyplot as plt + +print("Script for visual inspection of coordinate results for longitude and latitude") + +lon = Coordinate.from_resolution(Coordinates.longitude, 1,-180,180, True) +print(lon) + +print("\n") + +lat = Coordinate.from_resolution(Coordinates.latitude, 1.9,-90,90, False) +print(lat) + +testGrid = RegularGrid(lat_resolution=1.9, lon_resolution=2.5) + +testData = [x+0.1 for x in np.arange(-180,181, .5)] +testBins = testGrid.coord_to_index(Series(testData), lon.min, lon.step, len(lon.data) ) + +for val, bin in zip(testData, testBins): + print(val, "->", bin) + + +plt.plot(testData, testBins) +plt.show() \ No newline at end of file diff --git a/toargridding/grids.py b/toargridding/grids.py index 096b2fcdb955044ebf574c1b346b9fb889fec072..fc320300baaa027c2a30e4b9c78e48c07768b1c7 100644 --- a/toargridding/grids.py +++ b/toargridding/grids.py @@ -110,10 +110,10 @@ class RegularGrid(GridDefinition): # TODO make sure only sensible resolutions self.lat = Coordinate.from_resolution( - Coordinates.latitude, lat_resolution, min=-90, max=90 + Coordinates.latitude, lat_resolution, min=-90, max=90, wraps=False ) self.lon = Coordinate.from_resolution( - Coordinates.longitude, lon_resolution, min=-180, max=180 + Coordinates.longitude, lon_resolution, min=-180, max=180, wraps=True ) spatial_shape = (self.lon.size, self.lat.size) spatial_size = self.lon.size * self.lat.size @@ -272,13 +272,13 @@ class RegularGrid(GridDefinition): """ id_x = self.coord_to_index(coords[self.lat.name], self.lat.min, self.lat.step) - id_y = self.coord_to_index(coords[self.lon.name], self.lon.min, self.lon.step) + id_y = self.coord_to_index(coords[self.lon.name], self.lon.min, self.lon.step, len(self.lon.data)) id_i = self._as_i_index[id_x, id_y] return pd.Series(id_i, index=id_x.index) - def coord_to_index(self, coord : pd.Series, x0_axis : float, d_axis : float) -> np.array: + def coord_to_index(self, coord : pd.Series, x0_axis : float, d_axis : float, maxBin4Wrap : int=None) -> np.array: """converts a coordinate into a bin index on one axis Parameters: @@ -289,9 +289,14 @@ class RegularGrid(GridDefinition): offset of the axis d_axis: resolution of the axis + maxBin4Wrap: + this is the maximum bin number. If it is set, thw bin numbers wrap around, like for the latitude """ - - return (np.ceil((coord / d_axis) - 0.5) - x0_axis / d_axis).astype(int) + ids = (np.ceil((coord / d_axis) - 0.5) - x0_axis / d_axis).astype(int) + if maxBin4Wrap is not None: + ids[ids >= maxBin4Wrap] %= maxBin4Wrap + ids[ids < 0] += maxBin4Wrap + return ids def get_empty_grid(self, time: Variable, metadata: Metadata) -> xr.Dataset: # TODO make CF-compliant => docs """creation of an empty dataset without data @@ -319,4 +324,4 @@ class RegularGrid(GridDefinition): def get_id(self): """provide an ID for this grid and its resolution. """ - return f"_regular{self.lat.step}x{self.lon.step}" + return f"_regular{len(self.lat.data)}x{len(self.lon.data)}" diff --git a/toargridding/metadata.py b/toargridding/metadata.py index 0e7370ab3c09ba7dc1275307f61d284485306ff4..3689e04f4093b8e30bb74c6531aa426ccaaca717 100644 --- a/toargridding/metadata.py +++ b/toargridding/metadata.py @@ -234,7 +234,7 @@ def get_cf_metadata(variable: Variables, metadata: Metadata | None)-> Dict: "units": "degrees_north", "coverage_type": "coordinate", }, - "encoding": {"dtype": "int32"}, + "encoding": {"dtype": "float32"}, } case Variables.longitude.name: cf_metadata = { @@ -244,7 +244,7 @@ def get_cf_metadata(variable: Variables, metadata: Metadata | None)-> Dict: "units": "degrees_east", "coverage_type": "coordinate", }, - "encoding": {"dtype": "int32"}, + "encoding": {"dtype": "float32"}, } case Variables.time.name: cf_metadata = { diff --git a/toargridding/variables.py b/toargridding/variables.py index 5445cf543363e2f1228a9bebc0c452c9b22f2422..0440e0bf5de35a5cc14575943e0f419eb1d2344f 100644 --- a/toargridding/variables.py +++ b/toargridding/variables.py @@ -35,7 +35,7 @@ class Variable: """construction from analysis results """ cf_metadata = get_cf_metadata(variable, metadata=metadata) - # print(variable.name, cf_metadata) + #print(variable.name, cf_metadata) return cls(variable, data, **cf_metadata, **kwargs) def as_data_array(self, dims) -> xr.DataArray: @@ -84,7 +84,7 @@ class Coordinate(Variable): @classmethod def from_resolution( - cls, variable: Variables, resolution: float, min: float, max: float + cls, variable: Variables, resolution: float, min: float, max: float, wraps : bool ): """construction from a data range and resolution @@ -94,18 +94,31 @@ class Coordinate(Variable): ---------- resolution: width of a bin; actual size will be selected to obtain equidistant steps between steps + If the resolution does not result in a set of equidistant bins, a bin is added and the resolution is slightly decreased. min: lowest value on coordinate axis max: highest value on coordinate axis - + wraps: + if true, the last grid point is removed as it is the same as the first one, e.g. for longitude [-180, -180+res, 180-res]. + if false, the last grid point is kept like for the latitude [-90, -90+res... , 90-res,90] """ span = max - min - n = int(span / resolution) # raise error if invalid inputs ? + n = int(span / resolution) #TODO: raise error if invalid inputs ? + if n*resolution != span: + print(f"[DEBUG:] Resolution {resolution} does not provide an equidistant division of the span [{min},{max}]") + n+=1 + step = span / n + print(f"[DEBUG:] Adoption resolution {resolution} to {step}") + else: + step = resolution data = np.linspace(min, max, n + 1) - return cls.from_data(data, variable, None, step=resolution) + if wraps: + return cls.from_data(data[:-1], variable, None, step=step) + else: + return cls.from_data(data, variable, None, step=step) def as_data_array(self): """conversion to data array