diff --git a/tests/coordinates.py b/tests/coordinates.py new file mode 100644 index 0000000000000000000000000000000000000000..9781cf10bdb5730a25e029e3864ec83b1729eb81 --- /dev/null +++ b/tests/coordinates.py @@ -0,0 +1,39 @@ +#!/bin/python3 + +from toargridding.variables import Coordinate +from toargridding.metadata import Coordinates +import numpy as np +from pandas import Series + +def coord_to_index(coord : Series, x0_axis : float, d_axis : float, binWrap:int=None) -> np.array: + """converts a coordinate into a bin index on one axis + + Parameters: + ---------- + coord: + coordinate for conversion + x0_axis: + offset of the axis + d_axis: + resolution of the axis + """ + binIDs = (np.ceil((coord / d_axis) - 0.5) - x0_axis / d_axis).astype(int) + binIDs[binIDs >= binWrap] -=binWrap + return binIDs + + +print("Script for visual inspection of coordinate results for longitude and latitude") + +lon = Coordinate.from_resolution(Coordinates.longitude, 2.5,-180,180, True) +print(lon) + +testData = [x+0.1 for x in range(-180,180)] +testBins = coord_to_index(Series(testData), lon.min, lon.step, len(lon.data) ) +print(testBins) +print(len(lon.data)) +print(lon.data[144]) + +print("\n") + +#lat = Coordinate.from_resolution(Coordinates.latitude, 1.9,-90,90, False) +#print(lat) diff --git a/toargridding/grids.py b/toargridding/grids.py index 096b2fcdb955044ebf574c1b346b9fb889fec072..cfdd46263117f215ee3d1a73fea66ebd214afee6 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, max4binWrap : int=None) -> np.array: """converts a coordinate into a bin index on one axis Parameters: @@ -289,9 +289,13 @@ class RegularGrid(GridDefinition): offset of the axis d_axis: resolution of the axis + max4binWrap: + 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 max4binWrap is not None: + ids[ids >= max4binWrap] %= max4binWrap + 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 +323,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..4b3b712663b7481aded8fbdbd403ef31ece91c10 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,33 @@ 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 ? + 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) + #TODO need decision, if I need to remove the last point or not + #difference longitude and latitude! + 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