Merge branch '16-coordinate-ranges-do-not-match' into 'wip_tests_and_notebooks'

Resolve "coordinate ranges do not match"

See merge request !4

tests/coordinates.py

+#!/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

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)}"

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 = {

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