Feature/small fixes

polytope/datacube/backends/datacube.py

@@ -1,4 +1,5 @@
 import importlib
+import logging
 import math
 from abc import ABC, abstractmethod
 from typing import Any
@@ -105,6 +106,9 @@ def get_indices(self, path: DatacubePath, axis, lower, upper, method=None):
         if offset is not None:
             # Note that we can only do unique if not dealing with time values
             idx_between = unique(idx_between)
+
+        logging.info(f"For axis {axis.name} between {lower} and {upper}, found indices {idx_between}")
+
         return idx_between
 
     def _look_up_datacube(self, search_ranges, search_ranges_offset, indexes, axis, method):

polytope/datacube/backends/fdb.py

@@ -1,3 +1,4 @@
+import logging
 from copy import deepcopy
 
 import pygribjump as pygj
@@ -7,7 +8,14 @@
 
 
 class FDBDatacube(Datacube):
-    def __init__(self, config={}, axis_options={}):
+    def __init__(self, config=None, axis_options=None):
+        if config is None:
+            config = {}
+        if axis_options is None:
+            axis_options = {}
+
+        logging.info("Created an FDB datacube with options: " + str(axis_options))
+
         self.axis_options = axis_options
         self.axis_counter = 0
         self._axes = None
@@ -24,6 +32,9 @@ def __init__(self, config={}, axis_options={}):
 
         self.fdb = pygj.GribJump()
         self.fdb_coordinates = self.fdb.axes(partial_request)
+
+        logging.info("Axes returned from GribJump are: " + str(self.fdb_coordinates))
+
         self.fdb_coordinates["values"] = []
         for name, values in self.fdb_coordinates.items():
             values.sort()
@@ -39,16 +50,23 @@ def __init__(self, config={}, axis_options={}):
                 val = self._axes[name].type
                 self._check_and_add_axes(options, name, val)
 
+        logging.info("Polytope created axes for: " + str(self._axes.keys()))
+
     def get(self, requests: IndexTree):
         fdb_requests = []
         fdb_requests_decoding_info = []
         self.get_fdb_requests(requests, fdb_requests, fdb_requests_decoding_info)
         output_values = self.fdb.extract(fdb_requests)
         self.assign_fdb_output_to_nodes(output_values, fdb_requests_decoding_info)
 
-    def get_fdb_requests(self, requests: IndexTree, fdb_requests=[], fdb_requests_decoding_info=[], leaf_path={}):
+    def get_fdb_requests(self, requests: IndexTree, fdb_requests=[], fdb_requests_decoding_info=[], leaf_path=None):
+        if leaf_path is None:
+            leaf_path = {}
+
         # First when request node is root, go to its children
         if requests.axis.name == "root":
+            logging.info("Looking for data for the tree: " + str([leaf.flatten() for leaf in requests.leaves]))
+
             for c in requests.children:
                 self.get_fdb_requests(c, fdb_requests, fdb_requests_decoding_info)
         # If request node has no children, we have a leaf so need to assign fdb values to it
@@ -79,7 +97,9 @@ def get_fdb_requests(self, requests: IndexTree, fdb_requests=[], fdb_requests_de
                 for c in requests.children:
                     self.get_fdb_requests(c, fdb_requests, fdb_requests_decoding_info, leaf_path)
 
-    def get_2nd_last_values(self, requests, leaf_path={}):
+    def get_2nd_last_values(self, requests, leaf_path=None):
+        if leaf_path is None:
+            leaf_path = {}
         # In this function, we recursively loop over the last two layers of the tree and store the indices of the
         # request ranges in those layers
         # TODO: here find nearest point first before retrieving etc

polytope/datacube/backends/xarray.py

@@ -8,7 +8,9 @@
 class XArrayDatacube(Datacube):
     """Xarray arrays are labelled, axes can be defined as strings or integers (e.g. "time" or 0)."""
 
-    def __init__(self, dataarray: xr.DataArray, axis_options={}):
+    def __init__(self, dataarray: xr.DataArray, axis_options=None):
+        if axis_options is None:
+            axis_options = {}
         self.axis_options = axis_options
         self.axis_counter = 0
         self._axes = None

polytope/datacube/index_tree.py

@@ -1,6 +1,6 @@
 import json
-from typing import OrderedDict
 import logging
+from typing import OrderedDict
 
 from sortedcontainers import SortedList
 

polytope/datacube/transformations/__init__.py

@@ -0,0 +1 @@
+from ..transformations.datacube_transformations import *

polytope/datacube/transformations/datacube_mappers.py

@@ -1,4 +1,5 @@
 import bisect
+import logging
 import math
 from copy import deepcopy
 from importlib import import_module
@@ -126,6 +127,12 @@ def unmap(self, first_val, second_val):
         second_val = [i for i in self.second_axis_vals(first_val) if second_val - tol <= i <= second_val + tol][0]
         second_idx = self.second_axis_vals(first_val).index(second_val)
         final_index = self.axes_idx_to_regular_idx(first_idx, second_idx)
+
+        logging.info(
+            f"Mapped the values {first_val} on axis {self._mapped_axes[0]} \
+                and {second_val} on axis {self._mapped_axes[1]} to value {final_index} on axis {self._base_axis}"
+        )
+
         return final_index
 
 
@@ -1627,6 +1634,12 @@ def unmap(self, first_val, second_val):
         second_val = [i for i in self.second_axis_vals(first_val) if second_val - tol <= i <= second_val + tol][0]
         second_idx = self.second_axis_vals(first_val).index(second_val)
         reduced_ll_index = self.axes_idx_to_reduced_ll_idx(first_idx, second_idx)
+
+        logging.info(
+            f"Mapped the values {first_val} on axis {self._mapped_axes[0]} \
+                and {second_val} on axis {self._mapped_axes[1]} to value {reduced_ll_index} on axis {self._base_axis}"
+        )
+
         return reduced_ll_index
 
 
@@ -4494,6 +4507,12 @@ def find_second_axis_idx(self, first_val, second_val):
     def unmap(self, first_val, second_val):
         (first_idx, second_idx) = self.find_second_axis_idx(first_val, second_val)
         octahedral_index = self.axes_idx_to_octahedral_idx(first_idx, second_idx)
+
+        logging.info(
+            f"Mapped the values {first_val} on axis {self._mapped_axes[0]} \
+                and {second_val} on axis {self._mapped_axes[1]} to value {octahedral_index} on axis {self._base_axis}"
+        )
+
         return octahedral_index
 
 

polytope/datacube/transformations/datacube_merger.py

@@ -1,3 +1,5 @@
+import logging
+
 import numpy as np
 import pandas as pd
 
@@ -37,6 +39,10 @@ def merged_values(self, datacube):
                 val_to_add = val_to_add.astype("datetime64[s]")
                 merged_values.append(val_to_add)
         merged_values = np.array(merged_values)
+        logging.info(
+            f"Merged values {first_ax_vals} on axis {self.name} and \
+                     values {second_ax_vals} on axis {second_ax_name} to values {merged_values}"
+        )
         return merged_values
 
     def transformation_axes_final(self):
@@ -56,6 +62,10 @@ def unmerge(self, merged_val):
         # TODO: maybe replacing like this is too specific to time/dates?
         first_val = str(first_val).replace("-", "")
         second_val = second_val.replace(":", "")
+        logging.info(
+            f"Unmerged value {merged_val} to values {first_val} on axis {self.name} \
+                     and {second_val} on axis {self._second_axis}"
+        )
         return (first_val, second_val)
 
     def change_val_type(self, axis_name, values):

polytope/engine/hullslicer.py

@@ -1,3 +1,4 @@
+import logging
 import math
 from copy import copy
 from itertools import chain
@@ -74,6 +75,9 @@ def _build_sliceable_child(self, polytope, ax, node, datacube, lower, upper, nex
                 remapped_val_interm = ax.remap([value, value])[0]
                 remapped_val = (remapped_val_interm[0] + remapped_val_interm[1]) / 2
                 remapped_val = round(remapped_val, int(-math.log10(ax.tol)))
+
+            logging.info(f"Added index {remapped_val} on axis {ax.name} to the tree")
+
             child = node.create_child(ax, remapped_val)
             child["unsliced_polytopes"] = copy(node["unsliced_polytopes"])
             child["unsliced_polytopes"].remove(polytope)

polytope/polytope.py

@@ -37,10 +37,13 @@ def __repr__(self):
 
 
 class Polytope:
-    def __init__(self, datacube, engine=None, axis_options={}):
+    def __init__(self, datacube, engine=None, axis_options=None):
         from .datacube import Datacube
         from .engine import Engine
 
+        if axis_options is None:
+            axis_options = {}
+
         self.datacube = Datacube.create(datacube, axis_options)
         self.engine = engine if engine is not None else Engine.default()
 

requirements.txt

@@ -2,8 +2,8 @@ decorator==5.1.1
 numpy==1.23.5
 pandas==1.5.2
 pypi==2.1
-requests==2.28.1
-scipy==1.9.3
+requests==2.31.0
+scipy==1.11.4
 sortedcontainers==2.4.0
 tripy==1.0.0
 typing==3.7.4.3

tests/test_ecmwf_oper_data_fdb.py

@@ -3,7 +3,7 @@
 
 from polytope.engine.hullslicer import HullSlicer
 from polytope.polytope import Polytope, Request
-from polytope.shapes import Box, Select
+from polytope.shapes import Box, Point, Select
 
 
 class TestSlicingFDBDatacube:
@@ -39,3 +39,21 @@ def test_fdb_datacube(self):
         result = self.API.retrieve(request)
         result.pprint()
         assert len(result.leaves) == 9
+
+    @pytest.mark.fdb
+    def test_fdb_datacube_point(self):
+        request = Request(
+            Select("step", [0, 1]),
+            Select("levtype", ["sfc"]),
+            Select("date", [pd.Timestamp("20240103T0000")]),
+            Select("domain", ["g"]),
+            Select("expver", ["0001"]),
+            Select("param", ["167"]),
+            Select("class", ["od"]),
+            Select("stream", ["oper"]),
+            Select("type", ["fc"]),
+            Point(["latitude", "longitude"], [[0.035149384216, 0.0]], method="surrounding"),
+        )
+        result = self.API.retrieve(request)
+        result.pprint()
+        assert len(result.leaves) == 12

tests/test_point_shape.py

@@ -41,3 +41,8 @@ def test_point_surrounding_step(self):
         request = Request(Point(["step", "level"], [[2, 10]], method="surrounding"), Select("date", ["2000-01-01"]))
         result = self.API.retrieve(request)
         assert len(result.leaves) == 6
+
+    def test_point_surrounding_exact_step(self):
+        request = Request(Point(["step", "level"], [[3, 10]], method="surrounding"), Select("date", ["2000-01-01"]))
+        result = self.API.retrieve(request)
+        assert len(result.leaves) == 9

tests/test_slice_date_range_fdb.py

@@ -27,7 +27,6 @@ def setup_method(self, method):
     def test_fdb_datacube(self):
         request = Request(
             Select("step", [0]),
-            Select("number", [1]),
             Select("levtype", ["sfc"]),
             Span("date", pd.Timestamp("20230625T120000"), pd.Timestamp("20230626T120000")),
             Select("domain", ["g"]),
@@ -36,6 +35,7 @@ def test_fdb_datacube(self):
             Select("class", ["od"]),
             Select("stream", ["oper"]),
             Select("type", ["an"]),
+            Select("number", [1]),
             Box(["latitude", "longitude"], [0, 0], [0.2, 0.2]),
         )
         result = self.API.retrieve(request)