diff --git a/polytope/datacube/datacube_axis.py b/polytope/datacube/datacube_axis.py index 1119be94..78a22807 100644 --- a/polytope/datacube/datacube_axis.py +++ b/polytope/datacube/datacube_axis.py @@ -1,5 +1,3 @@ -import bisect -import math from abc import ABC, abstractmethod from copy import deepcopy from typing import Any, List @@ -7,533 +5,11 @@ import numpy as np import pandas as pd -from ..utility.list_tools import bisect_left_cmp, bisect_right_cmp - - -def cyclic(cls): - if cls.is_cyclic: - from .transformations.datacube_cyclic import DatacubeAxisCyclic - - def update_range(): - for transform in cls.transformations: - if isinstance(transform, DatacubeAxisCyclic): - transformation = transform - cls.range = transformation.range - - def to_intervals(range): - update_range() - if range[0] == -math.inf: - range[0] = cls.range[0] - if range[1] == math.inf: - range[1] = cls.range[1] - axis_lower = cls.range[0] - axis_upper = cls.range[1] - axis_range = axis_upper - axis_lower - lower = range[0] - upper = range[1] - intervals = [] - if lower < axis_upper: - # In this case, we want to go from lower to the first remapped cyclic axis upper - # or the asked upper range value. - # For example, if we have cyclic range [0,360] and we want to break [-270,180] into intervals, - # we first want to obtain [-270, 0] as the first range, where 0 is the remapped cyclic axis upper - # but if we wanted to break [-270, -180] into intervals, we would want to get [-270,-180], - # where -180 is the asked upper range value. - loops = int((axis_upper - lower) / axis_range) - remapped_up = axis_upper - (loops) * axis_range - new_upper = min(upper, remapped_up) - else: - # In this case, since lower >= axis_upper, we need to either go to the asked upper range - # or we need to go to the first remapped cyclic axis upper which is higher than lower - new_upper = min(axis_upper + axis_range, upper) - while new_upper < lower: - new_upper = min(new_upper + axis_range, upper) - intervals.append([lower, new_upper]) - # Now that we have established what the first interval should be, we should just jump from cyclic range - # to cyclic range until we hit the asked upper range value. - new_up = deepcopy(new_upper) - while new_up < upper: - new_upper = new_up - new_up = min(upper, new_upper + axis_range) - intervals.append([new_upper, new_up]) - # Once we have added all the in-between ranges, we need to add the last interval - intervals.append([new_up, upper]) - return intervals - - def _remap_range_to_axis_range(range): - update_range() - axis_lower = cls.range[0] - axis_upper = cls.range[1] - axis_range = axis_upper - axis_lower - lower = range[0] - upper = range[1] - if lower < axis_lower: - # In this case we need to calculate the number of loops between the axis lower - # and the lower to recenter the lower - loops = int((axis_lower - lower - cls.tol) / axis_range) - return_lower = lower + (loops + 1) * axis_range - return_upper = upper + (loops + 1) * axis_range - elif lower >= axis_upper: - # In this case we need to calculate the number of loops between the axis upper - # and the lower to recenter the lower - loops = int((lower - axis_upper) / axis_range) - return_lower = lower - (loops + 1) * axis_range - return_upper = upper - (loops + 1) * axis_range - else: - # In this case, the lower value is already in the right range - return_lower = lower - return_upper = upper - return [return_lower, return_upper] - - def _remap_val_to_axis_range(value): - return_range = _remap_range_to_axis_range([value, value]) - return return_range[0] - - def remap(range: List): - update_range() - if cls.range[0] - cls.tol <= range[0] <= cls.range[1] + cls.tol: - if cls.range[0] - cls.tol <= range[1] <= cls.range[1] + cls.tol: - # If we are already in the cyclic range, return it - return [range] - elif abs(range[0] - range[1]) <= 2 * cls.tol: - # If we have a range that is just one point, then it should still be counted - # and so we should take a small interval around it to find values inbetween - range = [ - _remap_val_to_axis_range(range[0]) - cls.tol, - _remap_val_to_axis_range(range[0]) + cls.tol, - ] - return [range] - range_intervals = cls.to_intervals(range) - ranges = [] - for interval in range_intervals: - if abs(interval[0] - interval[1]) > 0: - # If the interval is not just a single point, we remap it to the axis range - range = _remap_range_to_axis_range([interval[0], interval[1]]) - up = range[1] - low = range[0] - if up < low: - # Make sure we remap in the right order - ranges.append([up - cls.tol, low + cls.tol]) - else: - ranges.append([low - cls.tol, up + cls.tol]) - return ranges - - old_find_indexes = cls.find_indexes - - def find_indexes(path, datacube): - return old_find_indexes(path, datacube) - - old_unmap_path_key = cls.unmap_path_key - - def unmap_path_key(key_value_path, leaf_path, unwanted_path): - value = key_value_path[cls.name] - for transform in cls.transformations: - if isinstance(transform, DatacubeAxisCyclic): - if cls.name == transform.name: - new_val = _remap_val_to_axis_range(value) - key_value_path[cls.name] = new_val - key_value_path, leaf_path, unwanted_path = old_unmap_path_key(key_value_path, leaf_path, unwanted_path) - return (key_value_path, leaf_path, unwanted_path) - - old_unmap_to_datacube = cls.unmap_to_datacube - - def unmap_to_datacube(path, unmapped_path): - (path, unmapped_path) = old_unmap_to_datacube(path, unmapped_path) - return (path, unmapped_path) - - old_find_indices_between = cls.find_indices_between - - def find_indices_between(index_ranges, low, up, datacube, method=None): - update_range() - indexes_between_ranges = [] - - if method != "surrounding" and method != "nearest": - return old_find_indices_between(index_ranges, low, up, datacube, method) - else: - for indexes in index_ranges: - if cls.name in datacube.complete_axes: - start = indexes.searchsorted(low, "left") - end = indexes.searchsorted(up, "right") - else: - start = bisect.bisect_left(indexes, low) - end = bisect.bisect_right(indexes, up) - - if start - 1 < 0: - index_val_found = indexes[-1:][0] - indexes_between_ranges.append([index_val_found]) - if end + 1 > len(indexes): - index_val_found = indexes[:2][0] - indexes_between_ranges.append([index_val_found]) - start = max(start - 1, 0) - end = min(end + 1, len(indexes)) - if cls.name in datacube.complete_axes: - indexes_between = indexes[start:end].to_list() - else: - indexes_between = indexes[start:end] - indexes_between_ranges.append(indexes_between) - return indexes_between_ranges - - def offset(range): - # We first unpad the range by the axis tolerance to make sure that - # we find the wanted range of the cyclic axis since we padded by the axis tolerance before. - # Also, it's safer that we find the offset of a value inside the range instead of on the border - unpadded_range = [range[0] + 1.5 * cls.tol, range[1] - 1.5 * cls.tol] - cyclic_range = _remap_range_to_axis_range(unpadded_range) - offset = unpadded_range[0] - cyclic_range[0] - return offset - - cls.to_intervals = to_intervals - cls.remap = remap - cls.offset = offset - cls.find_indexes = find_indexes - cls.unmap_to_datacube = unmap_to_datacube - cls.find_indices_between = find_indices_between - cls.unmap_path_key = unmap_path_key - - return cls - - -def mapper(cls): - from .transformations.datacube_mappers import DatacubeMapper - - if cls.has_mapper: - - def find_indexes(path, datacube): - # first, find the relevant transformation object that is a mapping in the cls.transformation dico - for transform in cls.transformations: - if isinstance(transform, DatacubeMapper): - transformation = transform - if cls.name == transformation._mapped_axes()[0]: - return transformation.first_axis_vals() - if cls.name == transformation._mapped_axes()[1]: - first_val = path[transformation._mapped_axes()[0]] - return transformation.second_axis_vals(first_val) - - old_unmap_to_datacube = cls.unmap_to_datacube - - def unmap_to_datacube(path, unmapped_path): - (path, unmapped_path) = old_unmap_to_datacube(path, unmapped_path) - for transform in cls.transformations: - if isinstance(transform, DatacubeMapper): - if cls.name == transform._mapped_axes()[0]: - # if we are on the first axis, then need to add the first val to unmapped_path - first_val = path.get(cls.name, None) - path.pop(cls.name, None) - if cls.name not in unmapped_path: - # if for some reason, the unmapped_path already has the first axis val, then don't update - unmapped_path[cls.name] = first_val - if cls.name == transform._mapped_axes()[1]: - # if we are on the second axis, then the val of the first axis is stored - # inside unmapped_path so can get it from there - second_val = path.get(cls.name, None) - path.pop(cls.name, None) - first_val = unmapped_path.get(transform._mapped_axes()[0], None) - unmapped_path.pop(transform._mapped_axes()[0], None) - # if the first_val was not in the unmapped_path, then it's still in path - if first_val is None: - first_val = path.get(transform._mapped_axes()[0], None) - path.pop(transform._mapped_axes()[0], None) - if first_val is not None and second_val is not None: - unmapped_idx = transform.unmap(first_val, second_val) - unmapped_path[transform.old_axis] = unmapped_idx - return (path, unmapped_path) - - old_unmap_path_key = cls.unmap_path_key - - def unmap_path_key(key_value_path, leaf_path, unwanted_path): - key_value_path, leaf_path, unwanted_path = old_unmap_path_key(key_value_path, leaf_path, unwanted_path) - value = key_value_path[cls.name] - for transform in cls.transformations: - if isinstance(transform, DatacubeMapper): - if cls.name == transform._mapped_axes()[0]: - unwanted_val = key_value_path[transform._mapped_axes()[0]] - unwanted_path[cls.name] = unwanted_val - if cls.name == transform._mapped_axes()[1]: - first_val = unwanted_path[transform._mapped_axes()[0]] - unmapped_idx = transform.unmap(first_val, value) - leaf_path.pop(transform._mapped_axes()[0], None) - key_value_path.pop(cls.name) - key_value_path[transform.old_axis] = unmapped_idx - return (key_value_path, leaf_path, unwanted_path) - - def find_indices_between(index_ranges, low, up, datacube, method=None): - # TODO: add method for snappping - indexes_between_ranges = [] - for transform in cls.transformations: - if isinstance(transform, DatacubeMapper): - transformation = transform - if cls.name in transformation._mapped_axes(): - for idxs in index_ranges: - if method == "surrounding" or method == "nearest": - axis_reversed = transform._axis_reversed[cls.name] - if not axis_reversed: - start = bisect.bisect_left(idxs, low) - end = bisect.bisect_right(idxs, up) - else: - # TODO: do the custom bisect - end = bisect_left_cmp(idxs, low, cmp=lambda x, y: x > y) + 1 - start = bisect_right_cmp(idxs, up, cmp=lambda x, y: x > y) - start = max(start - 1, 0) - end = min(end + 1, len(idxs)) - indexes_between = idxs[start:end] - indexes_between_ranges.append(indexes_between) - else: - axis_reversed = transform._axis_reversed[cls.name] - if not axis_reversed: - lower_idx = bisect.bisect_left(idxs, low) - upper_idx = bisect.bisect_right(idxs, up) - indexes_between = idxs[lower_idx:upper_idx] - else: - # TODO: do the custom bisect - end_idx = bisect_left_cmp(idxs, low, cmp=lambda x, y: x > y) + 1 - start_idx = bisect_right_cmp(idxs, up, cmp=lambda x, y: x > y) - indexes_between = idxs[start_idx:end_idx] - indexes_between_ranges.append(indexes_between) - return indexes_between_ranges - - cls.find_indexes = find_indexes - cls.unmap_to_datacube = unmap_to_datacube - cls.find_indices_between = find_indices_between - cls.unmap_path_key = unmap_path_key - - return cls - - -def merge(cls): - from .transformations.datacube_merger import DatacubeAxisMerger - - if cls.has_merger: - - def find_indexes(path, datacube): - # first, find the relevant transformation object that is a mapping in the cls.transformation dico - for transform in cls.transformations: - if isinstance(transform, DatacubeAxisMerger): - transformation = transform - if cls.name == transformation._first_axis: - return transformation.merged_values(datacube) - - old_unmap_path_key = cls.unmap_path_key - - def unmap_path_key(key_value_path, leaf_path, unwanted_path): - key_value_path, leaf_path, unwanted_path = old_unmap_path_key(key_value_path, leaf_path, unwanted_path) - new_key_value_path = {} - value = key_value_path[cls.name] - for transform in cls.transformations: - if isinstance(transform, DatacubeAxisMerger): - if cls.name == transform._first_axis: - (first_val, second_val) = transform.unmerge(value) - new_key_value_path[transform._first_axis] = first_val - new_key_value_path[transform._second_axis] = second_val - return (new_key_value_path, leaf_path, unwanted_path) - - old_unmap_to_datacube = cls.unmap_to_datacube - - def unmap_to_datacube(path, unmapped_path): - (path, unmapped_path) = old_unmap_to_datacube(path, unmapped_path) - for transform in cls.transformations: - if isinstance(transform, DatacubeAxisMerger): - transformation = transform - if cls.name == transformation._first_axis: - old_val = path.get(cls.name, None) - (first_val, second_val) = transformation.unmerge(old_val) - path.pop(cls.name, None) - path[transformation._first_axis] = first_val - path[transformation._second_axis] = second_val - return (path, unmapped_path) - - def find_indices_between(index_ranges, low, up, datacube, method=None): - # TODO: add method for snappping - indexes_between_ranges = [] - for transform in cls.transformations: - if isinstance(transform, DatacubeAxisMerger): - transformation = transform - if cls.name in transformation._mapped_axes(): - for indexes in index_ranges: - if method == "surrounding" or method == "nearest": - start = indexes.index(low) - end = indexes.index(up) - start = max(start - 1, 0) - end = min(end + 1, len(indexes)) - indexes_between = indexes[start:end] - indexes_between_ranges.append(indexes_between) - else: - lower_idx = bisect.bisect_left(indexes, low) - upper_idx = bisect.bisect_right(indexes, up) - indexes_between = indexes[lower_idx:upper_idx] - indexes_between_ranges.append(indexes_between) - return indexes_between_ranges - - def remap(range): - return [range] - - cls.remap = remap - cls.find_indexes = find_indexes - cls.unmap_to_datacube = unmap_to_datacube - cls.find_indices_between = find_indices_between - cls.unmap_path_key = unmap_path_key - - return cls - - -def reverse(cls): - from .transformations.datacube_reverse import DatacubeAxisReverse - - if cls.reorder: - - def find_indexes(path, datacube): - # first, find the relevant transformation object that is a mapping in the cls.transformation dico - subarray = datacube.dataarray.sel(path, method="nearest") - unordered_indices = datacube.datacube_natural_indexes(cls, subarray) - if cls.name in datacube.complete_axes: - ordered_indices = unordered_indices.sort_values() - else: - ordered_indices = unordered_indices - return ordered_indices - - def find_indices_between(index_ranges, low, up, datacube, method=None): - # TODO: add method for snappping - indexes_between_ranges = [] - for transform in cls.transformations: - if isinstance(transform, DatacubeAxisReverse): - transformation = transform - if cls.name == transformation.name: - for indexes in index_ranges: - if cls.name in datacube.complete_axes: - # Find the range of indexes between lower and upper - # https://pandas.pydata.org/docs/reference/api/pandas.Index.searchsorted.html - # Assumes the indexes are already sorted (could sort to be sure) and monotonically - # increasing - if method == "surrounding" or method == "nearest": - start = indexes.searchsorted(low, "left") - end = indexes.searchsorted(up, "right") - start = max(start - 1, 0) - end = min(end + 1, len(indexes)) - indexes_between = indexes[start:end].to_list() - indexes_between_ranges.append(indexes_between) - else: - start = indexes.searchsorted(low, "left") - end = indexes.searchsorted(up, "right") - indexes_between = indexes[start:end].to_list() - indexes_between_ranges.append(indexes_between) - else: - if method == "surrounding" or method == "nearest": - start = indexes.index(low) - end = indexes.index(up) - start = max(start - 1, 0) - end = min(end + 1, len(indexes)) - indexes_between = indexes[start:end] - indexes_between_ranges.append(indexes_between) - else: - lower_idx = bisect.bisect_left(indexes, low) - upper_idx = bisect.bisect_right(indexes, up) - indexes_between = indexes[lower_idx:upper_idx] - indexes_between_ranges.append(indexes_between) - return indexes_between_ranges - - def remap(range): - return [range] - - cls.remap = remap - cls.find_indexes = find_indexes - cls.find_indices_between = find_indices_between - - return cls - - -def type_change(cls): - from .transformations.datacube_type_change import DatacubeAxisTypeChange - - if cls.type_change: - old_find_indexes = cls.find_indexes - - def find_indexes(path, datacube): - for transform in cls.transformations: - if isinstance(transform, DatacubeAxisTypeChange): - transformation = transform - if cls.name == transformation.name: - original_vals = old_find_indexes(path, datacube) - return transformation.change_val_type(cls.name, original_vals) - - old_unmap_path_key = cls.unmap_path_key - - def unmap_path_key(key_value_path, leaf_path, unwanted_path): - key_value_path, leaf_path, unwanted_path = old_unmap_path_key(key_value_path, leaf_path, unwanted_path) - value = key_value_path[cls.name] - for transform in cls.transformations: - if isinstance(transform, DatacubeAxisTypeChange): - if cls.name == transform.name: - unchanged_val = transform.make_str(value) - key_value_path[cls.name] = unchanged_val - return (key_value_path, leaf_path, unwanted_path) - - def unmap_to_datacube(path, unmapped_path): - for transform in cls.transformations: - if isinstance(transform, DatacubeAxisTypeChange): - transformation = transform - if cls.name == transformation.name: - changed_val = path.get(cls.name, None) - unchanged_val = transformation.make_str(changed_val) - if cls.name in path: - path.pop(cls.name, None) - unmapped_path[cls.name] = unchanged_val - return (path, unmapped_path) - - def find_indices_between(index_ranges, low, up, datacube, method=None): - # TODO: add method for snappping - indexes_between_ranges = [] - for transform in cls.transformations: - if isinstance(transform, DatacubeAxisTypeChange): - transformation = transform - if cls.name == transformation.name: - for indexes in index_ranges: - if method == "surrounding" or method == "nearest": - start = indexes.index(low) - end = indexes.index(up) - start = max(start - 1, 0) - end = min(end + 1, len(indexes)) - indexes_between = indexes[start:end] - indexes_between_ranges.append(indexes_between) - else: - lower_idx = bisect.bisect_left(indexes, low) - upper_idx = bisect.bisect_right(indexes, up) - indexes_between = indexes[lower_idx:upper_idx] - indexes_between_ranges.append(indexes_between) - return indexes_between_ranges - - def remap(range): - return [range] - - cls.remap = remap - cls.find_indexes = find_indexes - cls.unmap_to_datacube = unmap_to_datacube - cls.find_indices_between = find_indices_between - cls.unmap_path_key = unmap_path_key - - return cls - - -def null(cls): - if cls.type_change: - old_find_indexes = cls.find_indexes - - def find_indexes(path, datacube): - return old_find_indexes(path, datacube) - - def find_indices_between(index_ranges, low, up, datacube, method=None): - indexes_between_ranges = [] - for indexes in index_ranges: - indexes_between = [i for i in indexes if low <= i <= up] - indexes_between_ranges.append(indexes_between) - return indexes_between_ranges - - def remap(range): - return [range] - - cls.remap = remap - cls.find_indexes = find_indexes - cls.find_indices_between = find_indices_between - - return cls +from .transformations.datacube_cyclic.cyclic_axis_decorator import cyclic +from .transformations.datacube_mappers.mapper_axis_decorator import mapper +from .transformations.datacube_merger.merger_axis_decorator import merge +from .transformations.datacube_reverse.reverse_axis_decorator import reverse +from .transformations.datacube_type_change.type_change_axis_decorator import type_change class DatacubeAxis(ABC): diff --git a/polytope/datacube/index_tree.py b/polytope/datacube/index_tree.py index b15e67fe..b4ce3bce 100644 --- a/polytope/datacube/index_tree.py +++ b/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 @@ -46,6 +46,13 @@ def leaves_with_ancestors(self): self._collect_leaf_nodes(leaves) return leaves + def pprint_2(self, level=0): + if self.axis.name == "root": + print("\n") + print("\t" * level + "\u21b3" + str(self)) + for child in self.children: + child.pprint_2(level + 1) + def _collect_leaf_nodes_old(self, leaves): if len(self.children) == 0: leaves.append(self) diff --git a/polytope/datacube/transformations/datacube_cyclic/cyclic_axis_decorator.py b/polytope/datacube/transformations/datacube_cyclic/cyclic_axis_decorator.py new file mode 100644 index 00000000..76a2aa42 --- /dev/null +++ b/polytope/datacube/transformations/datacube_cyclic/cyclic_axis_decorator.py @@ -0,0 +1,188 @@ +import bisect +import math +from copy import deepcopy +from typing import List + +from .datacube_cyclic import DatacubeAxisCyclic + + +def cyclic(cls): + if cls.is_cyclic: + + def update_range(): + for transform in cls.transformations: + if isinstance(transform, DatacubeAxisCyclic): + transformation = transform + cls.range = transformation.range + + def to_intervals(range): + update_range() + if range[0] == -math.inf: + range[0] = cls.range[0] + if range[1] == math.inf: + range[1] = cls.range[1] + axis_lower = cls.range[0] + axis_upper = cls.range[1] + axis_range = axis_upper - axis_lower + lower = range[0] + upper = range[1] + intervals = [] + if lower < axis_upper: + # In this case, we want to go from lower to the first remapped cyclic axis upper + # or the asked upper range value. + # For example, if we have cyclic range [0,360] and we want to break [-270,180] into intervals, + # we first want to obtain [-270, 0] as the first range, where 0 is the remapped cyclic axis upper + # but if we wanted to break [-270, -180] into intervals, we would want to get [-270,-180], + # where -180 is the asked upper range value. + loops = int((axis_upper - lower) / axis_range) + remapped_up = axis_upper - (loops) * axis_range + new_upper = min(upper, remapped_up) + else: + # In this case, since lower >= axis_upper, we need to either go to the asked upper range + # or we need to go to the first remapped cyclic axis upper which is higher than lower + new_upper = min(axis_upper + axis_range, upper) + while new_upper < lower: + new_upper = min(new_upper + axis_range, upper) + intervals.append([lower, new_upper]) + # Now that we have established what the first interval should be, we should just jump from cyclic range + # to cyclic range until we hit the asked upper range value. + new_up = deepcopy(new_upper) + while new_up < upper: + new_upper = new_up + new_up = min(upper, new_upper + axis_range) + intervals.append([new_upper, new_up]) + # Once we have added all the in-between ranges, we need to add the last interval + intervals.append([new_up, upper]) + return intervals + + def _remap_range_to_axis_range(range): + update_range() + axis_lower = cls.range[0] + axis_upper = cls.range[1] + axis_range = axis_upper - axis_lower + lower = range[0] + upper = range[1] + if lower < axis_lower: + # In this case we need to calculate the number of loops between the axis lower + # and the lower to recenter the lower + loops = int((axis_lower - lower - cls.tol) / axis_range) + return_lower = lower + (loops + 1) * axis_range + return_upper = upper + (loops + 1) * axis_range + elif lower >= axis_upper: + # In this case we need to calculate the number of loops between the axis upper + # and the lower to recenter the lower + loops = int((lower - axis_upper) / axis_range) + return_lower = lower - (loops + 1) * axis_range + return_upper = upper - (loops + 1) * axis_range + else: + # In this case, the lower value is already in the right range + return_lower = lower + return_upper = upper + return [return_lower, return_upper] + + def _remap_val_to_axis_range(value): + return_range = _remap_range_to_axis_range([value, value]) + return return_range[0] + + def remap(range: List): + update_range() + if cls.range[0] - cls.tol <= range[0] <= cls.range[1] + cls.tol: + if cls.range[0] - cls.tol <= range[1] <= cls.range[1] + cls.tol: + # If we are already in the cyclic range, return it + return [range] + elif abs(range[0] - range[1]) <= 2 * cls.tol: + # If we have a range that is just one point, then it should still be counted + # and so we should take a small interval around it to find values inbetween + range = [ + _remap_val_to_axis_range(range[0]) - cls.tol, + _remap_val_to_axis_range(range[0]) + cls.tol, + ] + return [range] + range_intervals = cls.to_intervals(range) + ranges = [] + for interval in range_intervals: + if abs(interval[0] - interval[1]) > 0: + # If the interval is not just a single point, we remap it to the axis range + range = _remap_range_to_axis_range([interval[0], interval[1]]) + up = range[1] + low = range[0] + if up < low: + # Make sure we remap in the right order + ranges.append([up - cls.tol, low + cls.tol]) + else: + ranges.append([low - cls.tol, up + cls.tol]) + return ranges + + old_find_indexes = cls.find_indexes + + def find_indexes(path, datacube): + return old_find_indexes(path, datacube) + + old_unmap_path_key = cls.unmap_path_key + + def unmap_path_key(key_value_path, leaf_path, unwanted_path): + value = key_value_path[cls.name] + for transform in cls.transformations: + if isinstance(transform, DatacubeAxisCyclic): + if cls.name == transform.name: + new_val = _remap_val_to_axis_range(value) + key_value_path[cls.name] = new_val + key_value_path, leaf_path, unwanted_path = old_unmap_path_key(key_value_path, leaf_path, unwanted_path) + return (key_value_path, leaf_path, unwanted_path) + + old_unmap_to_datacube = cls.unmap_to_datacube + + def unmap_to_datacube(path, unmapped_path): + (path, unmapped_path) = old_unmap_to_datacube(path, unmapped_path) + return (path, unmapped_path) + + old_find_indices_between = cls.find_indices_between + + def find_indices_between(index_ranges, low, up, datacube, method=None): + update_range() + indexes_between_ranges = [] + + if method != "surrounding" or method != "nearest": + return old_find_indices_between(index_ranges, low, up, datacube, method) + else: + for indexes in index_ranges: + if cls.name in datacube.complete_axes: + start = indexes.searchsorted(low, "left") + end = indexes.searchsorted(up, "right") + else: + start = bisect.bisect_left(indexes, low) + end = bisect.bisect_right(indexes, up) + + if start - 1 < 0: + index_val_found = indexes[-1:][0] + indexes_between_ranges.append([index_val_found]) + if end + 1 > len(indexes): + index_val_found = indexes[:2][0] + indexes_between_ranges.append([index_val_found]) + start = max(start - 1, 0) + end = min(end + 1, len(indexes)) + if cls.name in datacube.complete_axes: + indexes_between = indexes[start:end].to_list() + else: + indexes_between = indexes[start:end] + indexes_between_ranges.append(indexes_between) + return indexes_between_ranges + + def offset(range): + # We first unpad the range by the axis tolerance to make sure that + # we find the wanted range of the cyclic axis since we padded by the axis tolerance before. + # Also, it's safer that we find the offset of a value inside the range instead of on the border + unpadded_range = [range[0] + 1.5 * cls.tol, range[1] - 1.5 * cls.tol] + cyclic_range = _remap_range_to_axis_range(unpadded_range) + offset = unpadded_range[0] - cyclic_range[0] + return offset + + cls.to_intervals = to_intervals + cls.remap = remap + cls.offset = offset + cls.find_indexes = find_indexes + cls.unmap_to_datacube = unmap_to_datacube + cls.find_indices_between = find_indices_between + cls.unmap_path_key = unmap_path_key + + return cls diff --git a/polytope/datacube/transformations/datacube_cyclic.py b/polytope/datacube/transformations/datacube_cyclic/datacube_cyclic.py similarity index 90% rename from polytope/datacube/transformations/datacube_cyclic.py rename to polytope/datacube/transformations/datacube_cyclic/datacube_cyclic.py index 802285c3..86113aa2 100644 --- a/polytope/datacube/transformations/datacube_cyclic.py +++ b/polytope/datacube/transformations/datacube_cyclic/datacube_cyclic.py @@ -1,4 +1,4 @@ -from .datacube_transformations import DatacubeAxisTransformation +from ..datacube_transformations import DatacubeAxisTransformation class DatacubeAxisCyclic(DatacubeAxisTransformation): diff --git a/polytope/datacube/transformations/datacube_mappers/datacube_mappers.py b/polytope/datacube/transformations/datacube_mappers/datacube_mappers.py new file mode 100644 index 00000000..a2034bfb --- /dev/null +++ b/polytope/datacube/transformations/datacube_mappers/datacube_mappers.py @@ -0,0 +1,85 @@ +from copy import deepcopy +from importlib import import_module + +from ..datacube_transformations import DatacubeAxisTransformation + + +class DatacubeMapper(DatacubeAxisTransformation): + # Needs to implements DatacubeAxisTransformation methods + + def __init__(self, name, mapper_options): + self.transformation_options = mapper_options + self.grid_type = mapper_options["type"] + self.grid_resolution = mapper_options["resolution"] + self.grid_axes = mapper_options["axes"] + self.local_area = [] + if "local" in mapper_options.keys(): + self.local_area = mapper_options["local"] + self.old_axis = name + self._final_transformation = self.generate_final_transformation() + self._final_mapped_axes = self._final_transformation._mapped_axes + self._axis_reversed = self._final_transformation._axis_reversed + + def generate_final_transformation(self): + map_type = _type_to_datacube_mapper_lookup[self.grid_type] + module = import_module("polytope.datacube.transformations.datacube_mappers.mapper_types." + self.grid_type) + constructor = getattr(module, map_type) + transformation = deepcopy(constructor(self.old_axis, self.grid_axes, self.grid_resolution, self.local_area)) + return transformation + + def blocked_axes(self): + return [] + + def unwanted_axes(self): + return [self._final_mapped_axes[0]] + + def transformation_axes_final(self): + final_axes = self._final_mapped_axes + return final_axes + + # Needs to also implement its own methods + + def change_val_type(self, axis_name, values): + # the new axis_vals created will be floats + return [0.0] + + def _mapped_axes(self): + # NOTE: Each of the mapper method needs to call it's sub mapper method + final_axes = self._final_mapped_axes + return final_axes + + def _base_axis(self): + pass + + def _resolution(self): + pass + + def first_axis_vals(self): + return self._final_transformation.first_axis_vals() + + def second_axis_vals(self, first_val): + return self._final_transformation.second_axis_vals(first_val) + + def map_first_axis(self, lower, upper): + return self._final_transformation.map_first_axis(lower, upper) + + def map_second_axis(self, first_val, lower, upper): + return self._final_transformation.map_second_axis(first_val, lower, upper) + + def find_second_idx(self, first_val, second_val): + return self._final_transformation.find_second_idx(first_val, second_val) + + def unmap_first_val_to_start_line_idx(self, first_val): + return self._final_transformation.unmap_first_val_to_start_line_idx(first_val) + + def unmap(self, first_val, second_val): + return self._final_transformation.unmap(first_val, second_val) + + +_type_to_datacube_mapper_lookup = { + "octahedral": "OctahedralGridMapper", + "healpix": "HealpixGridMapper", + "regular": "RegularGridMapper", + "reduced_ll": "ReducedLatLonMapper", + "local_regular": "LocalRegularGridMapper", +} diff --git a/polytope/datacube/transformations/datacube_mappers/mapper_axis_decorator.py b/polytope/datacube/transformations/datacube_mappers/mapper_axis_decorator.py new file mode 100644 index 00000000..99432186 --- /dev/null +++ b/polytope/datacube/transformations/datacube_mappers/mapper_axis_decorator.py @@ -0,0 +1,108 @@ +import bisect + +from ....utility.list_tools import bisect_left_cmp, bisect_right_cmp +from .datacube_mappers import DatacubeMapper + + +def mapper(cls): + if cls.has_mapper: + + def find_indexes(path, datacube): + # first, find the relevant transformation object that is a mapping in the cls.transformation dico + for transform in cls.transformations: + if isinstance(transform, DatacubeMapper): + transformation = transform + if cls.name == transformation._mapped_axes()[0]: + return transformation.first_axis_vals() + if cls.name == transformation._mapped_axes()[1]: + first_val = path[transformation._mapped_axes()[0]] + return transformation.second_axis_vals(first_val) + + old_unmap_to_datacube = cls.unmap_to_datacube + + def unmap_to_datacube(path, unmapped_path): + (path, unmapped_path) = old_unmap_to_datacube(path, unmapped_path) + for transform in cls.transformations: + if isinstance(transform, DatacubeMapper): + if cls.name == transform._mapped_axes()[0]: + # if we are on the first axis, then need to add the first val to unmapped_path + first_val = path.get(cls.name, None) + path.pop(cls.name, None) + if cls.name not in unmapped_path: + # if for some reason, the unmapped_path already has the first axis val, then don't update + unmapped_path[cls.name] = first_val + if cls.name == transform._mapped_axes()[1]: + # if we are on the second axis, then the val of the first axis is stored + # inside unmapped_path so can get it from there + second_val = path.get(cls.name, None) + path.pop(cls.name, None) + first_val = unmapped_path.get(transform._mapped_axes()[0], None) + unmapped_path.pop(transform._mapped_axes()[0], None) + # if the first_val was not in the unmapped_path, then it's still in path + if first_val is None: + first_val = path.get(transform._mapped_axes()[0], None) + path.pop(transform._mapped_axes()[0], None) + if first_val is not None and second_val is not None: + unmapped_idx = transform.unmap(first_val, second_val) + unmapped_path[transform.old_axis] = unmapped_idx + return (path, unmapped_path) + + old_unmap_path_key = cls.unmap_path_key + + def unmap_path_key(key_value_path, leaf_path, unwanted_path): + key_value_path, leaf_path, unwanted_path = old_unmap_path_key(key_value_path, leaf_path, unwanted_path) + value = key_value_path[cls.name] + for transform in cls.transformations: + if isinstance(transform, DatacubeMapper): + if cls.name == transform._mapped_axes()[0]: + unwanted_val = key_value_path[transform._mapped_axes()[0]] + unwanted_path[cls.name] = unwanted_val + if cls.name == transform._mapped_axes()[1]: + first_val = unwanted_path[transform._mapped_axes()[0]] + unmapped_idx = transform.unmap(first_val, value) + leaf_path.pop(transform._mapped_axes()[0], None) + key_value_path.pop(cls.name) + key_value_path[transform.old_axis] = unmapped_idx + return (key_value_path, leaf_path, unwanted_path) + + def find_indices_between(index_ranges, low, up, datacube, method=None): + # TODO: add method for snappping + indexes_between_ranges = [] + for transform in cls.transformations: + if isinstance(transform, DatacubeMapper): + transformation = transform + if cls.name in transformation._mapped_axes(): + for idxs in index_ranges: + if method == "surrounding" or method == "nearest": + axis_reversed = transform._axis_reversed[cls.name] + if not axis_reversed: + start = bisect.bisect_left(idxs, low) + end = bisect.bisect_right(idxs, up) + else: + # TODO: do the custom bisect + end = bisect_left_cmp(idxs, low, cmp=lambda x, y: x > y) + 1 + start = bisect_right_cmp(idxs, up, cmp=lambda x, y: x > y) + start = max(start - 1, 0) + end = min(end + 1, len(idxs)) + indexes_between = idxs[start:end] + indexes_between_ranges.append(indexes_between) + else: + axis_reversed = transform._axis_reversed[cls.name] + if not axis_reversed: + lower_idx = bisect.bisect_left(idxs, low) + upper_idx = bisect.bisect_right(idxs, up) + indexes_between = idxs[lower_idx:upper_idx] + else: + # TODO: do the custom bisect + end_idx = bisect_left_cmp(idxs, low, cmp=lambda x, y: x > y) + 1 + start_idx = bisect_right_cmp(idxs, up, cmp=lambda x, y: x > y) + indexes_between = idxs[start_idx:end_idx] + indexes_between_ranges.append(indexes_between) + return indexes_between_ranges + + cls.find_indexes = find_indexes + cls.unmap_to_datacube = unmap_to_datacube + cls.find_indices_between = find_indices_between + cls.unmap_path_key = unmap_path_key + + return cls diff --git a/polytope/datacube/transformations/datacube_mappers/mapper_types/healpix.py b/polytope/datacube/transformations/datacube_mappers/mapper_types/healpix.py new file mode 100644 index 00000000..8589ec71 --- /dev/null +++ b/polytope/datacube/transformations/datacube_mappers/mapper_types/healpix.py @@ -0,0 +1,124 @@ +import bisect +import math + +from ..datacube_mappers import DatacubeMapper + + +class HealpixGridMapper(DatacubeMapper): + def __init__(self, base_axis, mapped_axes, resolution, local_area=[]): + # TODO: if local area is not empty list, raise NotImplemented + self._mapped_axes = mapped_axes + self._base_axis = base_axis + self._resolution = resolution + self._axis_reversed = {mapped_axes[0]: True, mapped_axes[1]: False} + self._first_axis_vals = self.first_axis_vals() + + def first_axis_vals(self): + rad2deg = 180 / math.pi + vals = [0] * (4 * self._resolution - 1) + + # Polar caps + for i in range(1, self._resolution): + val = 90 - (rad2deg * math.acos(1 - (i * i / (3 * self._resolution * self._resolution)))) + vals[i - 1] = val + vals[4 * self._resolution - 1 - i] = -val + # Equatorial belts + for i in range(self._resolution, 2 * self._resolution): + val = 90 - (rad2deg * math.acos((4 * self._resolution - 2 * i) / (3 * self._resolution))) + vals[i - 1] = val + vals[4 * self._resolution - 1 - i] = -val + # Equator + vals[2 * self._resolution - 1] = 0 + return vals + + def map_first_axis(self, lower, upper): + axis_lines = self._first_axis_vals + return_vals = [val for val in axis_lines if lower <= val <= upper] + return return_vals + + def second_axis_vals(self, first_val): + tol = 1e-8 + first_val = [i for i in self._first_axis_vals if first_val - tol <= i <= first_val + tol][0] + idx = self._first_axis_vals.index(first_val) + + # Polar caps + if idx < self._resolution - 1 or 3 * self._resolution - 1 < idx <= 4 * self._resolution - 2: + start = 45 / (idx + 1) + vals = [start + i * (360 / (4 * (idx + 1))) for i in range(4 * (idx + 1))] + return vals + # Equatorial belts + start = 45 / self._resolution + if self._resolution - 1 <= idx < 2 * self._resolution - 1 or 2 * self._resolution <= idx < 3 * self._resolution: + r_start = start * (2 - (((idx + 1) - self._resolution + 1) % 2)) + vals = [r_start + i * (360 / (4 * self._resolution)) for i in range(4 * self._resolution)] + if vals[-1] == 360: + vals[-1] = 0 + return vals + # Equator + temp_val = 1 if self._resolution % 2 else 0 + r_start = start * (1 - temp_val) + if idx == 2 * self._resolution - 1: + vals = [r_start + i * (360 / (4 * self._resolution)) for i in range(4 * self._resolution)] + return vals + + def map_second_axis(self, first_val, lower, upper): + axis_lines = self.second_axis_vals(first_val) + return_vals = [val for val in axis_lines if lower <= val <= upper] + return return_vals + + def axes_idx_to_healpix_idx(self, first_idx, second_idx): + idx = 0 + for i in range(self._resolution - 1): + if i != first_idx: + idx += 4 * (i + 1) + else: + idx += second_idx + return idx + for i in range(self._resolution - 1, 3 * self._resolution): + if i != first_idx: + idx += 4 * self._resolution + else: + idx += second_idx + return idx + for i in range(3 * self._resolution, 4 * self._resolution - 1): + if i != first_idx: + idx += 4 * (4 * self._resolution - 1 - i + 1) + else: + idx += second_idx + return idx + + def find_second_idx(self, first_val, second_val): + tol = 1e-10 + second_axis_vals = self.second_axis_vals(first_val) + second_idx = bisect.bisect_left(second_axis_vals, second_val - tol) + return second_idx + + def unmap_first_val_to_start_line_idx(self, first_val): + tol = 1e-8 + first_val = [i for i in self._first_axis_vals if first_val - tol <= i <= first_val + tol][0] + first_idx = self._first_axis_vals.index(first_val) + idx = 0 + for i in range(self._resolution - 1): + if i != first_idx: + idx += 4 * (i + 1) + else: + return idx + for i in range(self._resolution - 1, 3 * self._resolution): + if i != first_idx: + idx += 4 * self._resolution + else: + return idx + for i in range(3 * self._resolution, 4 * self._resolution - 1): + if i != first_idx: + idx += 4 * (4 * self._resolution - 1 - i + 1) + else: + return idx + + def unmap(self, first_val, second_val): + tol = 1e-8 + first_val = [i for i in self._first_axis_vals if first_val - tol <= i <= first_val + tol][0] + first_idx = self._first_axis_vals.index(first_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) + healpix_index = self.axes_idx_to_healpix_idx(first_idx, second_idx) + return healpix_index diff --git a/polytope/datacube/transformations/datacube_mappers/mapper_types/local_regular.py b/polytope/datacube/transformations/datacube_mappers/mapper_types/local_regular.py new file mode 100644 index 00000000..a1514778 --- /dev/null +++ b/polytope/datacube/transformations/datacube_mappers/mapper_types/local_regular.py @@ -0,0 +1,69 @@ +import bisect + +from ..datacube_mappers import DatacubeMapper + + +class LocalRegularGridMapper(DatacubeMapper): + def __init__(self, base_axis, mapped_axes, resolution, local_area=[]): + # TODO: if local area is not empty list, raise NotImplemented + self._mapped_axes = mapped_axes + self._base_axis = base_axis + self._first_axis_min = local_area[0] + self._first_axis_max = local_area[1] + self._second_axis_min = local_area[2] + self._second_axis_max = local_area[3] + if not isinstance(resolution, list): + self.first_resolution = resolution + self.second_resolution = resolution + else: + self.first_resolution = resolution[0] + self.second_resolution = resolution[1] + self._first_deg_increment = (local_area[1] - local_area[0]) / self.first_resolution + self._second_deg_increment = (local_area[3] - local_area[2]) / self.second_resolution + self._axis_reversed = {mapped_axes[0]: True, mapped_axes[1]: False} + self._first_axis_vals = self.first_axis_vals() + + def first_axis_vals(self): + first_ax_vals = [self._first_axis_max - i * self._first_deg_increment for i in range(self.first_resolution + 1)] + return first_ax_vals + + def map_first_axis(self, lower, upper): + axis_lines = self._first_axis_vals + return_vals = [val for val in axis_lines if lower <= val <= upper] + return return_vals + + def second_axis_vals(self, first_val): + second_ax_vals = [ + self._second_axis_min + i * self._second_deg_increment for i in range(self.second_resolution + 1) + ] + return second_ax_vals + + def map_second_axis(self, first_val, lower, upper): + axis_lines = self.second_axis_vals(first_val) + return_vals = [val for val in axis_lines if lower <= val <= upper] + return return_vals + + def axes_idx_to_regular_idx(self, first_idx, second_idx): + final_idx = first_idx * (self.second_resolution + 1) + second_idx + return final_idx + + def find_second_idx(self, first_val, second_val): + tol = 1e-10 + second_axis_vals = self.second_axis_vals(first_val) + second_idx = bisect.bisect_left(second_axis_vals, second_val - tol) + return second_idx + + def unmap_first_val_to_start_line_idx(self, first_val): + tol = 1e-8 + first_val = [i for i in self._first_axis_vals if first_val - tol <= i <= first_val + tol][0] + first_idx = self._first_axis_vals.index(first_val) + return first_idx * self.second_resolution + + def unmap(self, first_val, second_val): + tol = 1e-8 + first_val = [i for i in self._first_axis_vals if first_val - tol <= i <= first_val + tol][0] + first_idx = self._first_axis_vals.index(first_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) + return final_index diff --git a/polytope/datacube/transformations/datacube_mappers.py b/polytope/datacube/transformations/datacube_mappers/mapper_types/octahedral.py similarity index 71% rename from polytope/datacube/transformations/datacube_mappers.py rename to polytope/datacube/transformations/datacube_mappers/mapper_types/octahedral.py index 0a7a3f5e..730ac959 100644 --- a/polytope/datacube/transformations/datacube_mappers.py +++ b/polytope/datacube/transformations/datacube_mappers/mapper_types/octahedral.py @@ -1,1756 +1,12 @@ -import bisect import math -from copy import deepcopy -from importlib import import_module -from ...utility.list_tools import bisect_left_cmp, bisect_right_cmp -from .datacube_transformations import DatacubeAxisTransformation - - -class DatacubeMapper(DatacubeAxisTransformation): - # Needs to implements DatacubeAxisTransformation methods - - def __init__(self, name, mapper_options): - self.transformation_options = mapper_options - self.grid_type = mapper_options["type"] - self.grid_resolution = mapper_options["resolution"] - self.grid_axes = mapper_options["axes"] - self.old_axis = name - self._final_transformation = self.generate_final_transformation() - self._final_mapped_axes = self._final_transformation._mapped_axes - self._axis_reversed = self._final_transformation._axis_reversed - - def generate_final_transformation(self): - map_type = _type_to_datacube_mapper_lookup[self.grid_type] - module = import_module("polytope.datacube.transformations.datacube_mappers") - constructor = getattr(module, map_type) - transformation = deepcopy(constructor(self.old_axis, self.grid_axes, self.grid_resolution)) - return transformation - - def blocked_axes(self): - return [] - - def unwanted_axes(self): - return [self._final_mapped_axes[0]] - - def transformation_axes_final(self): - final_axes = self._final_mapped_axes - return final_axes - - # Needs to also implement its own methods - - def change_val_type(self, axis_name, values): - # the new axis_vals created will be floats - return [0.0] - - def _mapped_axes(self): - # NOTE: Each of the mapper method needs to call it's sub mapper method - final_axes = self._final_mapped_axes - return final_axes - - def _base_axis(self): - pass - - def _resolution(self): - pass - - def first_axis_vals(self): - return self._final_transformation.first_axis_vals() - - def second_axis_vals(self, first_val): - return self._final_transformation.second_axis_vals(first_val) - - def map_first_axis(self, lower, upper): - return self._final_transformation.map_first_axis(lower, upper) - - def map_second_axis(self, first_val, lower, upper): - return self._final_transformation.map_second_axis(first_val, lower, upper) - - def find_second_idx(self, first_val, second_val): - return self._final_transformation.find_second_idx(first_val, second_val) - - def unmap_first_val_to_start_line_idx(self, first_val): - return self._final_transformation.unmap_first_val_to_start_line_idx(first_val) - - def unmap(self, first_val, second_val): - return self._final_transformation.unmap(first_val, second_val) - - -class RegularGridMapper(DatacubeMapper): - def __init__(self, base_axis, mapped_axes, resolution): - self._mapped_axes = mapped_axes - self._base_axis = base_axis - self._resolution = resolution - self.deg_increment = 90 / self._resolution - self._axis_reversed = {mapped_axes[0]: True, mapped_axes[1]: False} - self._first_axis_vals = self.first_axis_vals() - - def first_axis_vals(self): - first_ax_vals = [90 - i * self.deg_increment for i in range(2 * self._resolution)] - return first_ax_vals - - def map_first_axis(self, lower, upper): - axis_lines = self._first_axis_vals - return_vals = [val for val in axis_lines if lower <= val <= upper] - return return_vals - - def second_axis_vals(self, first_val): - second_ax_vals = [i * self.deg_increment for i in range(4 * self._resolution)] - return second_ax_vals - - def map_second_axis(self, first_val, lower, upper): - axis_lines = self.second_axis_vals(first_val) - return_vals = [val for val in axis_lines if lower <= val <= upper] - return return_vals - - def axes_idx_to_regular_idx(self, first_idx, second_idx): - final_idx = first_idx * 4 * self._resolution + second_idx - return final_idx - - def find_second_idx(self, first_val, second_val): - tol = 1e-10 - second_axis_vals = self.second_axis_vals(first_val) - second_idx = bisect.bisect_left(second_axis_vals, second_val - tol) - return second_idx - - def unmap_first_val_to_start_line_idx(self, first_val): - tol = 1e-8 - first_val = [i for i in self._first_axis_vals if first_val - tol <= i <= first_val + tol][0] - first_idx = self._first_axis_vals.index(first_val) - return first_idx * 4 * self._resolution - - def unmap(self, first_val, second_val): - tol = 1e-8 - first_val = [i for i in self._first_axis_vals if first_val - tol <= i <= first_val + tol][0] - first_idx = self._first_axis_vals.index(first_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) - return final_index - - -class ReducedLatLonMapper(DatacubeMapper): - def __init__(self, base_axis, mapped_axes, resolution): - self._mapped_axes = mapped_axes - self._base_axis = base_axis - self._resolution = resolution - self._axis_reversed = {mapped_axes[0]: False, mapped_axes[1]: False} - self._first_axis_vals = self.first_axis_vals() - - def first_axis_vals(self): - resolution = 180 / (self._resolution - 1) - vals = [-90 + i * resolution for i in range(self._resolution)] - return vals - - def map_first_axis(self, lower, upper): - axis_lines = self._first_axis_vals - return_vals = [val for val in axis_lines if lower <= val <= upper] - return return_vals - - def lon_spacing(self): - if self._resolution == 1441: - return [ - 2, - 6, - 14, - 20, - 26, - 32, - 38, - 44, - 50, - 58, - 64, - 70, - 76, - 82, - 88, - 94, - 102, - 108, - 114, - 120, - 126, - 132, - 138, - 144, - 152, - 158, - 164, - 170, - 176, - 182, - 188, - 196, - 202, - 208, - 214, - 220, - 226, - 232, - 238, - 246, - 252, - 258, - 264, - 270, - 276, - 282, - 290, - 296, - 302, - 308, - 314, - 320, - 326, - 332, - 340, - 346, - 352, - 358, - 364, - 370, - 376, - 382, - 388, - 396, - 402, - 408, - 414, - 420, - 426, - 432, - 438, - 444, - 452, - 458, - 464, - 470, - 476, - 482, - 488, - 494, - 500, - 506, - 512, - 520, - 526, - 532, - 538, - 544, - 550, - 556, - 562, - 568, - 574, - 580, - 586, - 594, - 600, - 606, - 612, - 618, - 624, - 630, - 636, - 642, - 648, - 654, - 660, - 666, - 672, - 678, - 686, - 692, - 698, - 704, - 710, - 716, - 722, - 728, - 734, - 740, - 746, - 752, - 758, - 764, - 770, - 776, - 782, - 788, - 794, - 800, - 806, - 812, - 818, - 824, - 830, - 836, - 842, - 848, - 854, - 860, - 866, - 872, - 878, - 884, - 890, - 896, - 902, - 908, - 914, - 920, - 926, - 932, - 938, - 944, - 950, - 956, - 962, - 968, - 974, - 980, - 986, - 992, - 998, - 1004, - 1010, - 1014, - 1020, - 1026, - 1032, - 1038, - 1044, - 1050, - 1056, - 1062, - 1068, - 1074, - 1080, - 1086, - 1092, - 1096, - 1102, - 1108, - 1114, - 1120, - 1126, - 1132, - 1138, - 1144, - 1148, - 1154, - 1160, - 1166, - 1172, - 1178, - 1184, - 1190, - 1194, - 1200, - 1206, - 1212, - 1218, - 1224, - 1230, - 1234, - 1240, - 1246, - 1252, - 1258, - 1264, - 1268, - 1274, - 1280, - 1286, - 1292, - 1296, - 1302, - 1308, - 1314, - 1320, - 1324, - 1330, - 1336, - 1342, - 1348, - 1352, - 1358, - 1364, - 1370, - 1374, - 1380, - 1386, - 1392, - 1396, - 1402, - 1408, - 1414, - 1418, - 1424, - 1430, - 1436, - 1440, - 1446, - 1452, - 1456, - 1462, - 1468, - 1474, - 1478, - 1484, - 1490, - 1494, - 1500, - 1506, - 1510, - 1516, - 1522, - 1526, - 1532, - 1538, - 1542, - 1548, - 1554, - 1558, - 1564, - 1570, - 1574, - 1580, - 1584, - 1590, - 1596, - 1600, - 1606, - 1610, - 1616, - 1622, - 1626, - 1632, - 1636, - 1642, - 1648, - 1652, - 1658, - 1662, - 1668, - 1672, - 1678, - 1684, - 1688, - 1694, - 1698, - 1704, - 1708, - 1714, - 1718, - 1724, - 1728, - 1734, - 1738, - 1744, - 1748, - 1754, - 1758, - 1764, - 1768, - 1774, - 1778, - 1784, - 1788, - 1794, - 1798, - 1804, - 1808, - 1812, - 1818, - 1822, - 1828, - 1832, - 1838, - 1842, - 1846, - 1852, - 1856, - 1862, - 1866, - 1870, - 1876, - 1880, - 1886, - 1890, - 1894, - 1900, - 1904, - 1908, - 1914, - 1918, - 1922, - 1928, - 1932, - 1936, - 1942, - 1946, - 1950, - 1956, - 1960, - 1964, - 1970, - 1974, - 1978, - 1982, - 1988, - 1992, - 1996, - 2002, - 2006, - 2010, - 2014, - 2020, - 2024, - 2028, - 2032, - 2036, - 2042, - 2046, - 2050, - 2054, - 2060, - 2064, - 2068, - 2072, - 2076, - 2080, - 2086, - 2090, - 2094, - 2098, - 2102, - 2106, - 2112, - 2116, - 2120, - 2124, - 2128, - 2132, - 2136, - 2140, - 2144, - 2150, - 2154, - 2158, - 2162, - 2166, - 2170, - 2174, - 2178, - 2182, - 2186, - 2190, - 2194, - 2198, - 2202, - 2206, - 2210, - 2214, - 2218, - 2222, - 2226, - 2230, - 2234, - 2238, - 2242, - 2246, - 2250, - 2254, - 2258, - 2262, - 2266, - 2270, - 2274, - 2278, - 2282, - 2286, - 2290, - 2292, - 2296, - 2300, - 2304, - 2308, - 2312, - 2316, - 2320, - 2324, - 2326, - 2330, - 2334, - 2338, - 2342, - 2346, - 2348, - 2352, - 2356, - 2360, - 2364, - 2366, - 2370, - 2374, - 2378, - 2382, - 2384, - 2388, - 2392, - 2396, - 2398, - 2402, - 2406, - 2410, - 2412, - 2416, - 2420, - 2422, - 2426, - 2430, - 2432, - 2436, - 2440, - 2442, - 2446, - 2450, - 2452, - 2456, - 2460, - 2462, - 2466, - 2470, - 2472, - 2476, - 2478, - 2482, - 2486, - 2488, - 2492, - 2494, - 2498, - 2500, - 2504, - 2508, - 2510, - 2514, - 2516, - 2520, - 2522, - 2526, - 2528, - 2532, - 2534, - 2538, - 2540, - 2544, - 2546, - 2550, - 2552, - 2556, - 2558, - 2560, - 2564, - 2566, - 2570, - 2572, - 2576, - 2578, - 2580, - 2584, - 2586, - 2590, - 2592, - 2594, - 2598, - 2600, - 2602, - 2606, - 2608, - 2610, - 2614, - 2616, - 2618, - 2622, - 2624, - 2626, - 2628, - 2632, - 2634, - 2636, - 2640, - 2642, - 2644, - 2646, - 2650, - 2652, - 2654, - 2656, - 2658, - 2662, - 2664, - 2666, - 2668, - 2670, - 2674, - 2676, - 2678, - 2680, - 2682, - 2684, - 2686, - 2690, - 2692, - 2694, - 2696, - 2698, - 2700, - 2702, - 2704, - 2706, - 2708, - 2712, - 2714, - 2716, - 2718, - 2720, - 2722, - 2724, - 2726, - 2728, - 2730, - 2732, - 2734, - 2736, - 2738, - 2740, - 2742, - 2744, - 2746, - 2748, - 2750, - 2750, - 2752, - 2754, - 2756, - 2758, - 2760, - 2762, - 2764, - 2766, - 2768, - 2768, - 2770, - 2772, - 2774, - 2776, - 2778, - 2780, - 2780, - 2782, - 2784, - 2786, - 2788, - 2788, - 2790, - 2792, - 2794, - 2794, - 2796, - 2798, - 2800, - 2800, - 2802, - 2804, - 2806, - 2806, - 2808, - 2810, - 2810, - 2812, - 2814, - 2814, - 2816, - 2818, - 2818, - 2820, - 2822, - 2822, - 2824, - 2826, - 2826, - 2828, - 2828, - 2830, - 2832, - 2832, - 2834, - 2834, - 2836, - 2836, - 2838, - 2838, - 2840, - 2842, - 2842, - 2844, - 2844, - 2846, - 2846, - 2846, - 2848, - 2848, - 2850, - 2850, - 2852, - 2852, - 2854, - 2854, - 2856, - 2856, - 2856, - 2858, - 2858, - 2860, - 2860, - 2860, - 2862, - 2862, - 2862, - 2864, - 2864, - 2864, - 2866, - 2866, - 2866, - 2868, - 2868, - 2868, - 2868, - 2870, - 2870, - 2870, - 2872, - 2872, - 2872, - 2872, - 2874, - 2874, - 2874, - 2874, - 2874, - 2876, - 2876, - 2876, - 2876, - 2876, - 2876, - 2878, - 2878, - 2878, - 2878, - 2878, - 2878, - 2878, - 2878, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2880, - 2878, - 2878, - 2878, - 2878, - 2878, - 2878, - 2878, - 2878, - 2876, - 2876, - 2876, - 2876, - 2876, - 2876, - 2874, - 2874, - 2874, - 2874, - 2874, - 2872, - 2872, - 2872, - 2872, - 2870, - 2870, - 2870, - 2868, - 2868, - 2868, - 2868, - 2866, - 2866, - 2866, - 2864, - 2864, - 2864, - 2862, - 2862, - 2862, - 2860, - 2860, - 2860, - 2858, - 2858, - 2856, - 2856, - 2856, - 2854, - 2854, - 2852, - 2852, - 2850, - 2850, - 2848, - 2848, - 2846, - 2846, - 2846, - 2844, - 2844, - 2842, - 2842, - 2840, - 2838, - 2838, - 2836, - 2836, - 2834, - 2834, - 2832, - 2832, - 2830, - 2828, - 2828, - 2826, - 2826, - 2824, - 2822, - 2822, - 2820, - 2818, - 2818, - 2816, - 2814, - 2814, - 2812, - 2810, - 2810, - 2808, - 2806, - 2806, - 2804, - 2802, - 2800, - 2800, - 2798, - 2796, - 2794, - 2794, - 2792, - 2790, - 2788, - 2788, - 2786, - 2784, - 2782, - 2780, - 2780, - 2778, - 2776, - 2774, - 2772, - 2770, - 2768, - 2768, - 2766, - 2764, - 2762, - 2760, - 2758, - 2756, - 2754, - 2752, - 2750, - 2750, - 2748, - 2746, - 2744, - 2742, - 2740, - 2738, - 2736, - 2734, - 2732, - 2730, - 2728, - 2726, - 2724, - 2722, - 2720, - 2718, - 2716, - 2714, - 2712, - 2708, - 2706, - 2704, - 2702, - 2700, - 2698, - 2696, - 2694, - 2692, - 2690, - 2686, - 2684, - 2682, - 2680, - 2678, - 2676, - 2674, - 2670, - 2668, - 2666, - 2664, - 2662, - 2658, - 2656, - 2654, - 2652, - 2650, - 2646, - 2644, - 2642, - 2640, - 2636, - 2634, - 2632, - 2628, - 2626, - 2624, - 2622, - 2618, - 2616, - 2614, - 2610, - 2608, - 2606, - 2602, - 2600, - 2598, - 2594, - 2592, - 2590, - 2586, - 2584, - 2580, - 2578, - 2576, - 2572, - 2570, - 2566, - 2564, - 2560, - 2558, - 2556, - 2552, - 2550, - 2546, - 2544, - 2540, - 2538, - 2534, - 2532, - 2528, - 2526, - 2522, - 2520, - 2516, - 2514, - 2510, - 2508, - 2504, - 2500, - 2498, - 2494, - 2492, - 2488, - 2486, - 2482, - 2478, - 2476, - 2472, - 2470, - 2466, - 2462, - 2460, - 2456, - 2452, - 2450, - 2446, - 2442, - 2440, - 2436, - 2432, - 2430, - 2426, - 2422, - 2420, - 2416, - 2412, - 2410, - 2406, - 2402, - 2398, - 2396, - 2392, - 2388, - 2384, - 2382, - 2378, - 2374, - 2370, - 2366, - 2364, - 2360, - 2356, - 2352, - 2348, - 2346, - 2342, - 2338, - 2334, - 2330, - 2326, - 2324, - 2320, - 2316, - 2312, - 2308, - 2304, - 2300, - 2296, - 2292, - 2290, - 2286, - 2282, - 2278, - 2274, - 2270, - 2266, - 2262, - 2258, - 2254, - 2250, - 2246, - 2242, - 2238, - 2234, - 2230, - 2226, - 2222, - 2218, - 2214, - 2210, - 2206, - 2202, - 2198, - 2194, - 2190, - 2186, - 2182, - 2178, - 2174, - 2170, - 2166, - 2162, - 2158, - 2154, - 2150, - 2144, - 2140, - 2136, - 2132, - 2128, - 2124, - 2120, - 2116, - 2112, - 2106, - 2102, - 2098, - 2094, - 2090, - 2086, - 2080, - 2076, - 2072, - 2068, - 2064, - 2060, - 2054, - 2050, - 2046, - 2042, - 2036, - 2032, - 2028, - 2024, - 2020, - 2014, - 2010, - 2006, - 2002, - 1996, - 1992, - 1988, - 1982, - 1978, - 1974, - 1970, - 1964, - 1960, - 1956, - 1950, - 1946, - 1942, - 1936, - 1932, - 1928, - 1922, - 1918, - 1914, - 1908, - 1904, - 1900, - 1894, - 1890, - 1886, - 1880, - 1876, - 1870, - 1866, - 1862, - 1856, - 1852, - 1846, - 1842, - 1838, - 1832, - 1828, - 1822, - 1818, - 1812, - 1808, - 1804, - 1798, - 1794, - 1788, - 1784, - 1778, - 1774, - 1768, - 1764, - 1758, - 1754, - 1748, - 1744, - 1738, - 1734, - 1728, - 1724, - 1718, - 1714, - 1708, - 1704, - 1698, - 1694, - 1688, - 1684, - 1678, - 1672, - 1668, - 1662, - 1658, - 1652, - 1648, - 1642, - 1636, - 1632, - 1626, - 1622, - 1616, - 1610, - 1606, - 1600, - 1596, - 1590, - 1584, - 1580, - 1574, - 1570, - 1564, - 1558, - 1554, - 1548, - 1542, - 1538, - 1532, - 1526, - 1522, - 1516, - 1510, - 1506, - 1500, - 1494, - 1490, - 1484, - 1478, - 1474, - 1468, - 1462, - 1456, - 1452, - 1446, - 1440, - 1436, - 1430, - 1424, - 1418, - 1414, - 1408, - 1402, - 1396, - 1392, - 1386, - 1380, - 1374, - 1370, - 1364, - 1358, - 1352, - 1348, - 1342, - 1336, - 1330, - 1324, - 1320, - 1314, - 1308, - 1302, - 1296, - 1292, - 1286, - 1280, - 1274, - 1268, - 1264, - 1258, - 1252, - 1246, - 1240, - 1234, - 1230, - 1224, - 1218, - 1212, - 1206, - 1200, - 1194, - 1190, - 1184, - 1178, - 1172, - 1166, - 1160, - 1154, - 1148, - 1144, - 1138, - 1132, - 1126, - 1120, - 1114, - 1108, - 1102, - 1096, - 1092, - 1086, - 1080, - 1074, - 1068, - 1062, - 1056, - 1050, - 1044, - 1038, - 1032, - 1026, - 1020, - 1014, - 1010, - 1004, - 998, - 992, - 986, - 980, - 974, - 968, - 962, - 956, - 950, - 944, - 938, - 932, - 926, - 920, - 914, - 908, - 902, - 896, - 890, - 884, - 878, - 872, - 866, - 860, - 854, - 848, - 842, - 836, - 830, - 824, - 818, - 812, - 806, - 800, - 794, - 788, - 782, - 776, - 770, - 764, - 758, - 752, - 746, - 740, - 734, - 728, - 722, - 716, - 710, - 704, - 698, - 692, - 686, - 678, - 672, - 666, - 660, - 654, - 648, - 642, - 636, - 630, - 624, - 618, - 612, - 606, - 600, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - ] - - def second_axis_vals(self, first_val): - first_idx = self._first_axis_vals.index(first_val) - Ny = self.lon_spacing()[first_idx] - second_spacing = 360 / Ny - return [i * second_spacing for i in range(Ny)] - - def map_second_axis(self, first_val, lower, upper): - axis_lines = self.second_axis_vals(first_val) - return_vals = [val for val in axis_lines if lower <= val <= upper] - return return_vals - - def axes_idx_to_reduced_ll_idx(self, first_idx, second_idx): - Ny_array = self.lon_spacing() - idx = 0 - for i in range(self._resolution): - if i != first_idx: - idx += Ny_array[i] - else: - idx += second_idx - return idx - - def find_second_idx(self, first_val, second_val): - tol = 1e-10 - second_axis_vals = self.second_axis_vals(first_val) - second_idx = bisect.bisect_left(second_axis_vals, second_val - tol) - return second_idx - - def unmap(self, first_val, second_val): - tol = 1e-8 - first_val = [i for i in self._first_axis_vals if first_val - tol <= i <= first_val + tol][0] - first_idx = self._first_axis_vals.index(first_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) - return reduced_ll_index - - -class HealpixGridMapper(DatacubeMapper): - def __init__(self, base_axis, mapped_axes, resolution): - self._mapped_axes = mapped_axes - self._base_axis = base_axis - self._resolution = resolution - self._axis_reversed = {mapped_axes[0]: True, mapped_axes[1]: False} - self._first_axis_vals = self.first_axis_vals() - - def first_axis_vals(self): - rad2deg = 180 / math.pi - vals = [0] * (4 * self._resolution - 1) - - # Polar caps - for i in range(1, self._resolution): - val = 90 - (rad2deg * math.acos(1 - (i * i / (3 * self._resolution * self._resolution)))) - vals[i - 1] = val - vals[4 * self._resolution - 1 - i] = -val - # Equatorial belts - for i in range(self._resolution, 2 * self._resolution): - val = 90 - (rad2deg * math.acos((4 * self._resolution - 2 * i) / (3 * self._resolution))) - vals[i - 1] = val - vals[4 * self._resolution - 1 - i] = -val - # Equator - vals[2 * self._resolution - 1] = 0 - return vals - - def map_first_axis(self, lower, upper): - axis_lines = self._first_axis_vals - return_vals = [val for val in axis_lines if lower <= val <= upper] - return return_vals - - def second_axis_vals(self, first_val): - tol = 1e-8 - first_val = [i for i in self._first_axis_vals if first_val - tol <= i <= first_val + tol][0] - idx = self._first_axis_vals.index(first_val) - - # Polar caps - if idx < self._resolution - 1 or 3 * self._resolution - 1 < idx <= 4 * self._resolution - 2: - start = 45 / (idx + 1) - vals = [start + i * (360 / (4 * (idx + 1))) for i in range(4 * (idx + 1))] - return vals - # Equatorial belts - start = 45 / self._resolution - if self._resolution - 1 <= idx < 2 * self._resolution - 1 or 2 * self._resolution <= idx < 3 * self._resolution: - r_start = start * (2 - (((idx + 1) - self._resolution + 1) % 2)) - vals = [r_start + i * (360 / (4 * self._resolution)) for i in range(4 * self._resolution)] - if vals[-1] == 360: - vals[-1] = 0 - return vals - # Equator - temp_val = 1 if self._resolution % 2 else 0 - r_start = start * (1 - temp_val) - if idx == 2 * self._resolution - 1: - vals = [r_start + i * (360 / (4 * self._resolution)) for i in range(4 * self._resolution)] - return vals - - def map_second_axis(self, first_val, lower, upper): - axis_lines = self.second_axis_vals(first_val) - return_vals = [val for val in axis_lines if lower <= val <= upper] - return return_vals - - def axes_idx_to_healpix_idx(self, first_idx, second_idx): - idx = 0 - for i in range(self._resolution - 1): - if i != first_idx: - idx += 4 * (i + 1) - else: - idx += second_idx - return idx - for i in range(self._resolution - 1, 3 * self._resolution): - if i != first_idx: - idx += 4 * self._resolution - else: - idx += second_idx - return idx - for i in range(3 * self._resolution, 4 * self._resolution - 1): - if i != first_idx: - idx += 4 * (4 * self._resolution - 1 - i + 1) - else: - idx += second_idx - return idx - - def find_second_idx(self, first_val, second_val): - tol = 1e-10 - second_axis_vals = self.second_axis_vals(first_val) - second_idx = bisect.bisect_left(second_axis_vals, second_val - tol) - return second_idx - - def unmap_first_val_to_start_line_idx(self, first_val): - tol = 1e-8 - first_val = [i for i in self._first_axis_vals if first_val - tol <= i <= first_val + tol][0] - first_idx = self._first_axis_vals.index(first_val) - idx = 0 - for i in range(self._resolution - 1): - if i != first_idx: - idx += 4 * (i + 1) - else: - return idx - for i in range(self._resolution - 1, 3 * self._resolution): - if i != first_idx: - idx += 4 * self._resolution - else: - return idx - for i in range(3 * self._resolution, 4 * self._resolution - 1): - if i != first_idx: - idx += 4 * (4 * self._resolution - 1 - i + 1) - else: - return idx - - def unmap(self, first_val, second_val): - tol = 1e-8 - first_val = [i for i in self._first_axis_vals if first_val - tol <= i <= first_val + tol][0] - first_idx = self._first_axis_vals.index(first_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) - healpix_index = self.axes_idx_to_healpix_idx(first_idx, second_idx) - return healpix_index +from .....utility.list_tools import bisect_left_cmp, bisect_right_cmp +from ..datacube_mappers import DatacubeMapper class OctahedralGridMapper(DatacubeMapper): - def __init__(self, base_axis, mapped_axes, resolution): + def __init__(self, base_axis, mapped_axes, resolution, local_area=[]): + # TODO: if local area is not empty list, raise NotImplemented self._mapped_axes = mapped_axes self._base_axis = base_axis self._resolution = resolution @@ -4495,11 +2751,3 @@ 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) return octahedral_index - - -_type_to_datacube_mapper_lookup = { - "octahedral": "OctahedralGridMapper", - "healpix": "HealpixGridMapper", - "regular": "RegularGridMapper", - "reduced_ll": "ReducedLatLonMapper", -} diff --git a/polytope/datacube/transformations/datacube_mappers/mapper_types/reduced_ll.py b/polytope/datacube/transformations/datacube_mappers/mapper_types/reduced_ll.py new file mode 100644 index 00000000..ece09b4a --- /dev/null +++ b/polytope/datacube/transformations/datacube_mappers/mapper_types/reduced_ll.py @@ -0,0 +1,1505 @@ +import bisect + +from ..datacube_mappers import DatacubeMapper + + +class ReducedLatLonMapper(DatacubeMapper): + def __init__(self, base_axis, mapped_axes, resolution, local_area=[]): + # TODO: if local area is not empty list, raise NotImplemented + self._mapped_axes = mapped_axes + self._base_axis = base_axis + self._resolution = resolution + self._axis_reversed = {mapped_axes[0]: False, mapped_axes[1]: False} + self._first_axis_vals = self.first_axis_vals() + + def first_axis_vals(self): + resolution = 180 / (self._resolution - 1) + vals = [-90 + i * resolution for i in range(self._resolution)] + return vals + + def map_first_axis(self, lower, upper): + axis_lines = self._first_axis_vals + return_vals = [val for val in axis_lines if lower <= val <= upper] + return return_vals + + def lon_spacing(self): + if self._resolution == 1441: + return [ + 2, + 6, + 14, + 20, + 26, + 32, + 38, + 44, + 50, + 58, + 64, + 70, + 76, + 82, + 88, + 94, + 102, + 108, + 114, + 120, + 126, + 132, + 138, + 144, + 152, + 158, + 164, + 170, + 176, + 182, + 188, + 196, + 202, + 208, + 214, + 220, + 226, + 232, + 238, + 246, + 252, + 258, + 264, + 270, + 276, + 282, + 290, + 296, + 302, + 308, + 314, + 320, + 326, + 332, + 340, + 346, + 352, + 358, + 364, + 370, + 376, + 382, + 388, + 396, + 402, + 408, + 414, + 420, + 426, + 432, + 438, + 444, + 452, + 458, + 464, + 470, + 476, + 482, + 488, + 494, + 500, + 506, + 512, + 520, + 526, + 532, + 538, + 544, + 550, + 556, + 562, + 568, + 574, + 580, + 586, + 594, + 600, + 606, + 612, + 618, + 624, + 630, + 636, + 642, + 648, + 654, + 660, + 666, + 672, + 678, + 686, + 692, + 698, + 704, + 710, + 716, + 722, + 728, + 734, + 740, + 746, + 752, + 758, + 764, + 770, + 776, + 782, + 788, + 794, + 800, + 806, + 812, + 818, + 824, + 830, + 836, + 842, + 848, + 854, + 860, + 866, + 872, + 878, + 884, + 890, + 896, + 902, + 908, + 914, + 920, + 926, + 932, + 938, + 944, + 950, + 956, + 962, + 968, + 974, + 980, + 986, + 992, + 998, + 1004, + 1010, + 1014, + 1020, + 1026, + 1032, + 1038, + 1044, + 1050, + 1056, + 1062, + 1068, + 1074, + 1080, + 1086, + 1092, + 1096, + 1102, + 1108, + 1114, + 1120, + 1126, + 1132, + 1138, + 1144, + 1148, + 1154, + 1160, + 1166, + 1172, + 1178, + 1184, + 1190, + 1194, + 1200, + 1206, + 1212, + 1218, + 1224, + 1230, + 1234, + 1240, + 1246, + 1252, + 1258, + 1264, + 1268, + 1274, + 1280, + 1286, + 1292, + 1296, + 1302, + 1308, + 1314, + 1320, + 1324, + 1330, + 1336, + 1342, + 1348, + 1352, + 1358, + 1364, + 1370, + 1374, + 1380, + 1386, + 1392, + 1396, + 1402, + 1408, + 1414, + 1418, + 1424, + 1430, + 1436, + 1440, + 1446, + 1452, + 1456, + 1462, + 1468, + 1474, + 1478, + 1484, + 1490, + 1494, + 1500, + 1506, + 1510, + 1516, + 1522, + 1526, + 1532, + 1538, + 1542, + 1548, + 1554, + 1558, + 1564, + 1570, + 1574, + 1580, + 1584, + 1590, + 1596, + 1600, + 1606, + 1610, + 1616, + 1622, + 1626, + 1632, + 1636, + 1642, + 1648, + 1652, + 1658, + 1662, + 1668, + 1672, + 1678, + 1684, + 1688, + 1694, + 1698, + 1704, + 1708, + 1714, + 1718, + 1724, + 1728, + 1734, + 1738, + 1744, + 1748, + 1754, + 1758, + 1764, + 1768, + 1774, + 1778, + 1784, + 1788, + 1794, + 1798, + 1804, + 1808, + 1812, + 1818, + 1822, + 1828, + 1832, + 1838, + 1842, + 1846, + 1852, + 1856, + 1862, + 1866, + 1870, + 1876, + 1880, + 1886, + 1890, + 1894, + 1900, + 1904, + 1908, + 1914, + 1918, + 1922, + 1928, + 1932, + 1936, + 1942, + 1946, + 1950, + 1956, + 1960, + 1964, + 1970, + 1974, + 1978, + 1982, + 1988, + 1992, + 1996, + 2002, + 2006, + 2010, + 2014, + 2020, + 2024, + 2028, + 2032, + 2036, + 2042, + 2046, + 2050, + 2054, + 2060, + 2064, + 2068, + 2072, + 2076, + 2080, + 2086, + 2090, + 2094, + 2098, + 2102, + 2106, + 2112, + 2116, + 2120, + 2124, + 2128, + 2132, + 2136, + 2140, + 2144, + 2150, + 2154, + 2158, + 2162, + 2166, + 2170, + 2174, + 2178, + 2182, + 2186, + 2190, + 2194, + 2198, + 2202, + 2206, + 2210, + 2214, + 2218, + 2222, + 2226, + 2230, + 2234, + 2238, + 2242, + 2246, + 2250, + 2254, + 2258, + 2262, + 2266, + 2270, + 2274, + 2278, + 2282, + 2286, + 2290, + 2292, + 2296, + 2300, + 2304, + 2308, + 2312, + 2316, + 2320, + 2324, + 2326, + 2330, + 2334, + 2338, + 2342, + 2346, + 2348, + 2352, + 2356, + 2360, + 2364, + 2366, + 2370, + 2374, + 2378, + 2382, + 2384, + 2388, + 2392, + 2396, + 2398, + 2402, + 2406, + 2410, + 2412, + 2416, + 2420, + 2422, + 2426, + 2430, + 2432, + 2436, + 2440, + 2442, + 2446, + 2450, + 2452, + 2456, + 2460, + 2462, + 2466, + 2470, + 2472, + 2476, + 2478, + 2482, + 2486, + 2488, + 2492, + 2494, + 2498, + 2500, + 2504, + 2508, + 2510, + 2514, + 2516, + 2520, + 2522, + 2526, + 2528, + 2532, + 2534, + 2538, + 2540, + 2544, + 2546, + 2550, + 2552, + 2556, + 2558, + 2560, + 2564, + 2566, + 2570, + 2572, + 2576, + 2578, + 2580, + 2584, + 2586, + 2590, + 2592, + 2594, + 2598, + 2600, + 2602, + 2606, + 2608, + 2610, + 2614, + 2616, + 2618, + 2622, + 2624, + 2626, + 2628, + 2632, + 2634, + 2636, + 2640, + 2642, + 2644, + 2646, + 2650, + 2652, + 2654, + 2656, + 2658, + 2662, + 2664, + 2666, + 2668, + 2670, + 2674, + 2676, + 2678, + 2680, + 2682, + 2684, + 2686, + 2690, + 2692, + 2694, + 2696, + 2698, + 2700, + 2702, + 2704, + 2706, + 2708, + 2712, + 2714, + 2716, + 2718, + 2720, + 2722, + 2724, + 2726, + 2728, + 2730, + 2732, + 2734, + 2736, + 2738, + 2740, + 2742, + 2744, + 2746, + 2748, + 2750, + 2750, + 2752, + 2754, + 2756, + 2758, + 2760, + 2762, + 2764, + 2766, + 2768, + 2768, + 2770, + 2772, + 2774, + 2776, + 2778, + 2780, + 2780, + 2782, + 2784, + 2786, + 2788, + 2788, + 2790, + 2792, + 2794, + 2794, + 2796, + 2798, + 2800, + 2800, + 2802, + 2804, + 2806, + 2806, + 2808, + 2810, + 2810, + 2812, + 2814, + 2814, + 2816, + 2818, + 2818, + 2820, + 2822, + 2822, + 2824, + 2826, + 2826, + 2828, + 2828, + 2830, + 2832, + 2832, + 2834, + 2834, + 2836, + 2836, + 2838, + 2838, + 2840, + 2842, + 2842, + 2844, + 2844, + 2846, + 2846, + 2846, + 2848, + 2848, + 2850, + 2850, + 2852, + 2852, + 2854, + 2854, + 2856, + 2856, + 2856, + 2858, + 2858, + 2860, + 2860, + 2860, + 2862, + 2862, + 2862, + 2864, + 2864, + 2864, + 2866, + 2866, + 2866, + 2868, + 2868, + 2868, + 2868, + 2870, + 2870, + 2870, + 2872, + 2872, + 2872, + 2872, + 2874, + 2874, + 2874, + 2874, + 2874, + 2876, + 2876, + 2876, + 2876, + 2876, + 2876, + 2878, + 2878, + 2878, + 2878, + 2878, + 2878, + 2878, + 2878, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2880, + 2878, + 2878, + 2878, + 2878, + 2878, + 2878, + 2878, + 2878, + 2876, + 2876, + 2876, + 2876, + 2876, + 2876, + 2874, + 2874, + 2874, + 2874, + 2874, + 2872, + 2872, + 2872, + 2872, + 2870, + 2870, + 2870, + 2868, + 2868, + 2868, + 2868, + 2866, + 2866, + 2866, + 2864, + 2864, + 2864, + 2862, + 2862, + 2862, + 2860, + 2860, + 2860, + 2858, + 2858, + 2856, + 2856, + 2856, + 2854, + 2854, + 2852, + 2852, + 2850, + 2850, + 2848, + 2848, + 2846, + 2846, + 2846, + 2844, + 2844, + 2842, + 2842, + 2840, + 2838, + 2838, + 2836, + 2836, + 2834, + 2834, + 2832, + 2832, + 2830, + 2828, + 2828, + 2826, + 2826, + 2824, + 2822, + 2822, + 2820, + 2818, + 2818, + 2816, + 2814, + 2814, + 2812, + 2810, + 2810, + 2808, + 2806, + 2806, + 2804, + 2802, + 2800, + 2800, + 2798, + 2796, + 2794, + 2794, + 2792, + 2790, + 2788, + 2788, + 2786, + 2784, + 2782, + 2780, + 2780, + 2778, + 2776, + 2774, + 2772, + 2770, + 2768, + 2768, + 2766, + 2764, + 2762, + 2760, + 2758, + 2756, + 2754, + 2752, + 2750, + 2750, + 2748, + 2746, + 2744, + 2742, + 2740, + 2738, + 2736, + 2734, + 2732, + 2730, + 2728, + 2726, + 2724, + 2722, + 2720, + 2718, + 2716, + 2714, + 2712, + 2708, + 2706, + 2704, + 2702, + 2700, + 2698, + 2696, + 2694, + 2692, + 2690, + 2686, + 2684, + 2682, + 2680, + 2678, + 2676, + 2674, + 2670, + 2668, + 2666, + 2664, + 2662, + 2658, + 2656, + 2654, + 2652, + 2650, + 2646, + 2644, + 2642, + 2640, + 2636, + 2634, + 2632, + 2628, + 2626, + 2624, + 2622, + 2618, + 2616, + 2614, + 2610, + 2608, + 2606, + 2602, + 2600, + 2598, + 2594, + 2592, + 2590, + 2586, + 2584, + 2580, + 2578, + 2576, + 2572, + 2570, + 2566, + 2564, + 2560, + 2558, + 2556, + 2552, + 2550, + 2546, + 2544, + 2540, + 2538, + 2534, + 2532, + 2528, + 2526, + 2522, + 2520, + 2516, + 2514, + 2510, + 2508, + 2504, + 2500, + 2498, + 2494, + 2492, + 2488, + 2486, + 2482, + 2478, + 2476, + 2472, + 2470, + 2466, + 2462, + 2460, + 2456, + 2452, + 2450, + 2446, + 2442, + 2440, + 2436, + 2432, + 2430, + 2426, + 2422, + 2420, + 2416, + 2412, + 2410, + 2406, + 2402, + 2398, + 2396, + 2392, + 2388, + 2384, + 2382, + 2378, + 2374, + 2370, + 2366, + 2364, + 2360, + 2356, + 2352, + 2348, + 2346, + 2342, + 2338, + 2334, + 2330, + 2326, + 2324, + 2320, + 2316, + 2312, + 2308, + 2304, + 2300, + 2296, + 2292, + 2290, + 2286, + 2282, + 2278, + 2274, + 2270, + 2266, + 2262, + 2258, + 2254, + 2250, + 2246, + 2242, + 2238, + 2234, + 2230, + 2226, + 2222, + 2218, + 2214, + 2210, + 2206, + 2202, + 2198, + 2194, + 2190, + 2186, + 2182, + 2178, + 2174, + 2170, + 2166, + 2162, + 2158, + 2154, + 2150, + 2144, + 2140, + 2136, + 2132, + 2128, + 2124, + 2120, + 2116, + 2112, + 2106, + 2102, + 2098, + 2094, + 2090, + 2086, + 2080, + 2076, + 2072, + 2068, + 2064, + 2060, + 2054, + 2050, + 2046, + 2042, + 2036, + 2032, + 2028, + 2024, + 2020, + 2014, + 2010, + 2006, + 2002, + 1996, + 1992, + 1988, + 1982, + 1978, + 1974, + 1970, + 1964, + 1960, + 1956, + 1950, + 1946, + 1942, + 1936, + 1932, + 1928, + 1922, + 1918, + 1914, + 1908, + 1904, + 1900, + 1894, + 1890, + 1886, + 1880, + 1876, + 1870, + 1866, + 1862, + 1856, + 1852, + 1846, + 1842, + 1838, + 1832, + 1828, + 1822, + 1818, + 1812, + 1808, + 1804, + 1798, + 1794, + 1788, + 1784, + 1778, + 1774, + 1768, + 1764, + 1758, + 1754, + 1748, + 1744, + 1738, + 1734, + 1728, + 1724, + 1718, + 1714, + 1708, + 1704, + 1698, + 1694, + 1688, + 1684, + 1678, + 1672, + 1668, + 1662, + 1658, + 1652, + 1648, + 1642, + 1636, + 1632, + 1626, + 1622, + 1616, + 1610, + 1606, + 1600, + 1596, + 1590, + 1584, + 1580, + 1574, + 1570, + 1564, + 1558, + 1554, + 1548, + 1542, + 1538, + 1532, + 1526, + 1522, + 1516, + 1510, + 1506, + 1500, + 1494, + 1490, + 1484, + 1478, + 1474, + 1468, + 1462, + 1456, + 1452, + 1446, + 1440, + 1436, + 1430, + 1424, + 1418, + 1414, + 1408, + 1402, + 1396, + 1392, + 1386, + 1380, + 1374, + 1370, + 1364, + 1358, + 1352, + 1348, + 1342, + 1336, + 1330, + 1324, + 1320, + 1314, + 1308, + 1302, + 1296, + 1292, + 1286, + 1280, + 1274, + 1268, + 1264, + 1258, + 1252, + 1246, + 1240, + 1234, + 1230, + 1224, + 1218, + 1212, + 1206, + 1200, + 1194, + 1190, + 1184, + 1178, + 1172, + 1166, + 1160, + 1154, + 1148, + 1144, + 1138, + 1132, + 1126, + 1120, + 1114, + 1108, + 1102, + 1096, + 1092, + 1086, + 1080, + 1074, + 1068, + 1062, + 1056, + 1050, + 1044, + 1038, + 1032, + 1026, + 1020, + 1014, + 1010, + 1004, + 998, + 992, + 986, + 980, + 974, + 968, + 962, + 956, + 950, + 944, + 938, + 932, + 926, + 920, + 914, + 908, + 902, + 896, + 890, + 884, + 878, + 872, + 866, + 860, + 854, + 848, + 842, + 836, + 830, + 824, + 818, + 812, + 806, + 800, + 794, + 788, + 782, + 776, + 770, + 764, + 758, + 752, + 746, + 740, + 734, + 728, + 722, + 716, + 710, + 704, + 698, + 692, + 686, + 678, + 672, + 666, + 660, + 654, + 648, + 642, + 636, + 630, + 624, + 618, + 612, + 606, + 600, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + ] + + def second_axis_vals(self, first_val): + first_idx = self._first_axis_vals.index(first_val) + Ny = self.lon_spacing()[first_idx] + second_spacing = 360 / Ny + return [i * second_spacing for i in range(Ny)] + + def map_second_axis(self, first_val, lower, upper): + axis_lines = self.second_axis_vals(first_val) + return_vals = [val for val in axis_lines if lower <= val <= upper] + return return_vals + + def axes_idx_to_reduced_ll_idx(self, first_idx, second_idx): + Ny_array = self.lon_spacing() + idx = 0 + for i in range(self._resolution): + if i != first_idx: + idx += Ny_array[i] + else: + idx += second_idx + return idx + + def find_second_idx(self, first_val, second_val): + tol = 1e-10 + second_axis_vals = self.second_axis_vals(first_val) + second_idx = bisect.bisect_left(second_axis_vals, second_val - tol) + return second_idx + + def unmap(self, first_val, second_val): + tol = 1e-8 + first_val = [i for i in self._first_axis_vals if first_val - tol <= i <= first_val + tol][0] + first_idx = self._first_axis_vals.index(first_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) + return reduced_ll_index diff --git a/polytope/datacube/transformations/datacube_mappers/mapper_types/regular.py b/polytope/datacube/transformations/datacube_mappers/mapper_types/regular.py new file mode 100644 index 00000000..c8f207fc --- /dev/null +++ b/polytope/datacube/transformations/datacube_mappers/mapper_types/regular.py @@ -0,0 +1,57 @@ +import bisect + +from ..datacube_mappers import DatacubeMapper + + +class RegularGridMapper(DatacubeMapper): + def __init__(self, base_axis, mapped_axes, resolution, local_area=[]): + # TODO: if local area is not empty list, raise NotImplemented + self._mapped_axes = mapped_axes + self._base_axis = base_axis + self._resolution = resolution + self.deg_increment = 90 / self._resolution + self._axis_reversed = {mapped_axes[0]: True, mapped_axes[1]: False} + self._first_axis_vals = self.first_axis_vals() + + def first_axis_vals(self): + first_ax_vals = [90 - i * self.deg_increment for i in range(2 * self._resolution)] + return first_ax_vals + + def map_first_axis(self, lower, upper): + axis_lines = self._first_axis_vals + return_vals = [val for val in axis_lines if lower <= val <= upper] + return return_vals + + def second_axis_vals(self, first_val): + second_ax_vals = [i * self.deg_increment for i in range(4 * self._resolution)] + return second_ax_vals + + def map_second_axis(self, first_val, lower, upper): + axis_lines = self.second_axis_vals(first_val) + return_vals = [val for val in axis_lines if lower <= val <= upper] + return return_vals + + def axes_idx_to_regular_idx(self, first_idx, second_idx): + final_idx = first_idx * 4 * self._resolution + second_idx + return final_idx + + def find_second_idx(self, first_val, second_val): + tol = 1e-10 + second_axis_vals = self.second_axis_vals(first_val) + second_idx = bisect.bisect_left(second_axis_vals, second_val - tol) + return second_idx + + def unmap_first_val_to_start_line_idx(self, first_val): + tol = 1e-8 + first_val = [i for i in self._first_axis_vals if first_val - tol <= i <= first_val + tol][0] + first_idx = self._first_axis_vals.index(first_val) + return first_idx * 4 * self._resolution + + def unmap(self, first_val, second_val): + tol = 1e-8 + first_val = [i for i in self._first_axis_vals if first_val - tol <= i <= first_val + tol][0] + first_idx = self._first_axis_vals.index(first_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) + return final_index diff --git a/polytope/datacube/transformations/datacube_merger.py b/polytope/datacube/transformations/datacube_merger/datacube_merger.py similarity index 97% rename from polytope/datacube/transformations/datacube_merger.py rename to polytope/datacube/transformations/datacube_merger/datacube_merger.py index d6027867..ecd53372 100644 --- a/polytope/datacube/transformations/datacube_merger.py +++ b/polytope/datacube/transformations/datacube_merger/datacube_merger.py @@ -1,7 +1,7 @@ import numpy as np import pandas as pd -from .datacube_transformations import DatacubeAxisTransformation +from ..datacube_transformations import DatacubeAxisTransformation class DatacubeAxisMerger(DatacubeAxisTransformation): diff --git a/polytope/datacube/transformations/datacube_merger/merger_axis_decorator.py b/polytope/datacube/transformations/datacube_merger/merger_axis_decorator.py new file mode 100644 index 00000000..16d81639 --- /dev/null +++ b/polytope/datacube/transformations/datacube_merger/merger_axis_decorator.py @@ -0,0 +1,77 @@ +import bisect + +from .datacube_merger import DatacubeAxisMerger + + +def merge(cls): + if cls.has_merger: + + def find_indexes(path, datacube): + # first, find the relevant transformation object that is a mapping in the cls.transformation dico + for transform in cls.transformations: + if isinstance(transform, DatacubeAxisMerger): + transformation = transform + if cls.name == transformation._first_axis: + return transformation.merged_values(datacube) + + old_unmap_path_key = cls.unmap_path_key + + def unmap_path_key(key_value_path, leaf_path, unwanted_path): + key_value_path, leaf_path, unwanted_path = old_unmap_path_key(key_value_path, leaf_path, unwanted_path) + new_key_value_path = {} + value = key_value_path[cls.name] + for transform in cls.transformations: + if isinstance(transform, DatacubeAxisMerger): + if cls.name == transform._first_axis: + (first_val, second_val) = transform.unmerge(value) + new_key_value_path[transform._first_axis] = first_val + new_key_value_path[transform._second_axis] = second_val + return (new_key_value_path, leaf_path, unwanted_path) + + old_unmap_to_datacube = cls.unmap_to_datacube + + def unmap_to_datacube(path, unmapped_path): + (path, unmapped_path) = old_unmap_to_datacube(path, unmapped_path) + for transform in cls.transformations: + if isinstance(transform, DatacubeAxisMerger): + transformation = transform + if cls.name == transformation._first_axis: + old_val = path.get(cls.name, None) + (first_val, second_val) = transformation.unmerge(old_val) + path.pop(cls.name, None) + path[transformation._first_axis] = first_val + path[transformation._second_axis] = second_val + return (path, unmapped_path) + + def find_indices_between(index_ranges, low, up, datacube, method=None): + # TODO: add method for snappping + indexes_between_ranges = [] + for transform in cls.transformations: + if isinstance(transform, DatacubeAxisMerger): + transformation = transform + if cls.name in transformation._mapped_axes(): + for indexes in index_ranges: + if method == "surrounding" or method == "nearest": + start = indexes.index(low) + end = indexes.index(up) + start = max(start - 1, 0) + end = min(end + 1, len(indexes)) + indexes_between = indexes[start:end] + indexes_between_ranges.append(indexes_between) + else: + lower_idx = bisect.bisect_left(indexes, low) + upper_idx = bisect.bisect_right(indexes, up) + indexes_between = indexes[lower_idx:upper_idx] + indexes_between_ranges.append(indexes_between) + return indexes_between_ranges + + def remap(range): + return [range] + + cls.remap = remap + cls.find_indexes = find_indexes + cls.unmap_to_datacube = unmap_to_datacube + cls.find_indices_between = find_indices_between + cls.unmap_path_key = unmap_path_key + + return cls diff --git a/polytope/datacube/transformations/datacube_null_transformation.py b/polytope/datacube/transformations/datacube_null_transformation/datacube_null_transformation.py similarity index 88% rename from polytope/datacube/transformations/datacube_null_transformation.py rename to polytope/datacube/transformations/datacube_null_transformation/datacube_null_transformation.py index 55c94277..43dccbbe 100644 --- a/polytope/datacube/transformations/datacube_null_transformation.py +++ b/polytope/datacube/transformations/datacube_null_transformation/datacube_null_transformation.py @@ -1,4 +1,4 @@ -from .datacube_transformations import DatacubeAxisTransformation +from ..datacube_transformations import DatacubeAxisTransformation class DatacubeNullTransformation(DatacubeAxisTransformation): diff --git a/polytope/datacube/transformations/datacube_null_transformation/null_axis_decorator.py b/polytope/datacube/transformations/datacube_null_transformation/null_axis_decorator.py new file mode 100644 index 00000000..10e2644d --- /dev/null +++ b/polytope/datacube/transformations/datacube_null_transformation/null_axis_decorator.py @@ -0,0 +1,22 @@ +def null(cls): + if cls.type_change: + old_find_indexes = cls.find_indexes + + def find_indexes(path, datacube): + return old_find_indexes(path, datacube) + + def find_indices_between(index_ranges, low, up, datacube, method=None): + indexes_between_ranges = [] + for indexes in index_ranges: + indexes_between = [i for i in indexes if low <= i <= up] + indexes_between_ranges.append(indexes_between) + return indexes_between_ranges + + def remap(range): + return [range] + + cls.remap = remap + cls.find_indexes = find_indexes + cls.find_indices_between = find_indices_between + + return cls diff --git a/polytope/datacube/transformations/datacube_reverse.py b/polytope/datacube/transformations/datacube_reverse/datacube_reverse.py similarity index 88% rename from polytope/datacube/transformations/datacube_reverse.py rename to polytope/datacube/transformations/datacube_reverse/datacube_reverse.py index 6a556907..6e60de87 100644 --- a/polytope/datacube/transformations/datacube_reverse.py +++ b/polytope/datacube/transformations/datacube_reverse/datacube_reverse.py @@ -1,4 +1,4 @@ -from .datacube_transformations import DatacubeAxisTransformation +from ..datacube_transformations import DatacubeAxisTransformation class DatacubeAxisReverse(DatacubeAxisTransformation): diff --git a/polytope/datacube/transformations/datacube_reverse/reverse_axis_decorator.py b/polytope/datacube/transformations/datacube_reverse/reverse_axis_decorator.py new file mode 100644 index 00000000..18bc8bd6 --- /dev/null +++ b/polytope/datacube/transformations/datacube_reverse/reverse_axis_decorator.py @@ -0,0 +1,66 @@ +import bisect + +from .datacube_reverse import DatacubeAxisReverse + + +def reverse(cls): + if cls.reorder: + + def find_indexes(path, datacube): + # first, find the relevant transformation object that is a mapping in the cls.transformation dico + subarray = datacube.dataarray.sel(path, method="nearest") + unordered_indices = datacube.datacube_natural_indexes(cls, subarray) + if cls.name in datacube.complete_axes: + ordered_indices = unordered_indices.sort_values() + else: + ordered_indices = unordered_indices + return ordered_indices + + def find_indices_between(index_ranges, low, up, datacube, method=None): + # TODO: add method for snappping + indexes_between_ranges = [] + for transform in cls.transformations: + if isinstance(transform, DatacubeAxisReverse): + transformation = transform + if cls.name == transformation.name: + for indexes in index_ranges: + if cls.name in datacube.complete_axes: + # Find the range of indexes between lower and upper + # https://pandas.pydata.org/docs/reference/api/pandas.Index.searchsorted.html + # Assumes the indexes are already sorted (could sort to be sure) and monotonically + # increasing + if method == "surrounding" or method == "nearest": + start = indexes.searchsorted(low, "left") + end = indexes.searchsorted(up, "right") + start = max(start - 1, 0) + end = min(end + 1, len(indexes)) + indexes_between = indexes[start:end].to_list() + indexes_between_ranges.append(indexes_between) + else: + start = indexes.searchsorted(low, "left") + end = indexes.searchsorted(up, "right") + indexes_between = indexes[start:end].to_list() + indexes_between_ranges.append(indexes_between) + else: + if method == "surrounding" or method == "nearest": + start = indexes.index(low) + end = indexes.index(up) + start = max(start - 1, 0) + end = min(end + 1, len(indexes)) + indexes_between = indexes[start:end] + indexes_between_ranges.append(indexes_between) + else: + lower_idx = bisect.bisect_left(indexes, low) + upper_idx = bisect.bisect_right(indexes, up) + indexes_between = indexes[lower_idx:upper_idx] + indexes_between_ranges.append(indexes_between) + return indexes_between_ranges + + def remap(range): + return [range] + + cls.remap = remap + cls.find_indexes = find_indexes + cls.find_indices_between = find_indices_between + + return cls diff --git a/polytope/datacube/transformations/datacube_transformations.py b/polytope/datacube/transformations/datacube_transformations.py index 900ad16b..2077f346 100644 --- a/polytope/datacube/transformations/datacube_transformations.py +++ b/polytope/datacube/transformations/datacube_transformations.py @@ -8,8 +8,8 @@ class DatacubeAxisTransformation(ABC): def create_transform(name, transformation_type_key, transformation_options): transformation_type = _type_to_datacube_transformation_lookup[transformation_type_key] transformation_file_name = _type_to_transformation_file_lookup[transformation_type_key] - - module = import_module("polytope.datacube.transformations.datacube_" + transformation_file_name) + file_name = ".datacube_" + transformation_file_name + module = import_module("polytope.datacube.transformations" + file_name + file_name) constructor = getattr(module, transformation_type) transformation_type_option = transformation_options[transformation_type_key] new_transformation = deepcopy(constructor(name, transformation_type_option)) diff --git a/polytope/datacube/transformations/datacube_type_change.py b/polytope/datacube/transformations/datacube_type_change/datacube_type_change.py similarity index 93% rename from polytope/datacube/transformations/datacube_type_change.py rename to polytope/datacube/transformations/datacube_type_change/datacube_type_change.py index cdc046b7..8ba2ef0f 100644 --- a/polytope/datacube/transformations/datacube_type_change.py +++ b/polytope/datacube/transformations/datacube_type_change/datacube_type_change.py @@ -1,7 +1,7 @@ from copy import deepcopy from importlib import import_module -from .datacube_transformations import DatacubeAxisTransformation +from ..datacube_transformations import DatacubeAxisTransformation class DatacubeAxisTypeChange(DatacubeAxisTransformation): @@ -15,7 +15,7 @@ def __init__(self, name, type_options): def generate_final_transformation(self): map_type = _type_to_datacube_type_change_lookup[self.new_type] - module = import_module("polytope.datacube.transformations.datacube_type_change") + module = import_module("polytope.datacube.transformations.datacube_type_change.datacube_type_change") constructor = getattr(module, map_type) transformation = deepcopy(constructor(self.name, self.new_type)) return transformation diff --git a/polytope/datacube/transformations/datacube_type_change/type_change_axis_decorator.py b/polytope/datacube/transformations/datacube_type_change/type_change_axis_decorator.py new file mode 100644 index 00000000..0e366982 --- /dev/null +++ b/polytope/datacube/transformations/datacube_type_change/type_change_axis_decorator.py @@ -0,0 +1,73 @@ +import bisect + +from .datacube_type_change import DatacubeAxisTypeChange + + +def type_change(cls): + if cls.type_change: + old_find_indexes = cls.find_indexes + + def find_indexes(path, datacube): + for transform in cls.transformations: + if isinstance(transform, DatacubeAxisTypeChange): + transformation = transform + if cls.name == transformation.name: + original_vals = old_find_indexes(path, datacube) + return transformation.change_val_type(cls.name, original_vals) + + old_unmap_path_key = cls.unmap_path_key + + def unmap_path_key(key_value_path, leaf_path, unwanted_path): + key_value_path, leaf_path, unwanted_path = old_unmap_path_key(key_value_path, leaf_path, unwanted_path) + value = key_value_path[cls.name] + for transform in cls.transformations: + if isinstance(transform, DatacubeAxisTypeChange): + if cls.name == transform.name: + unchanged_val = transform.make_str(value) + key_value_path[cls.name] = unchanged_val + return (key_value_path, leaf_path, unwanted_path) + + def unmap_to_datacube(path, unmapped_path): + for transform in cls.transformations: + if isinstance(transform, DatacubeAxisTypeChange): + transformation = transform + if cls.name == transformation.name: + changed_val = path.get(cls.name, None) + unchanged_val = transformation.make_str(changed_val) + if cls.name in path: + path.pop(cls.name, None) + unmapped_path[cls.name] = unchanged_val + return (path, unmapped_path) + + def find_indices_between(index_ranges, low, up, datacube, method=None): + # TODO: add method for snappping + indexes_between_ranges = [] + for transform in cls.transformations: + if isinstance(transform, DatacubeAxisTypeChange): + transformation = transform + if cls.name == transformation.name: + for indexes in index_ranges: + if method == "surrounding" or method == "nearest": + start = indexes.index(low) + end = indexes.index(up) + start = max(start - 1, 0) + end = min(end + 1, len(indexes)) + indexes_between = indexes[start:end] + indexes_between_ranges.append(indexes_between) + else: + lower_idx = bisect.bisect_left(indexes, low) + upper_idx = bisect.bisect_right(indexes, up) + indexes_between = indexes[lower_idx:upper_idx] + indexes_between_ranges.append(indexes_between) + return indexes_between_ranges + + def remap(range): + return [range] + + cls.remap = remap + cls.find_indexes = find_indexes + cls.unmap_to_datacube = unmap_to_datacube + cls.find_indices_between = find_indices_between + cls.unmap_path_key = unmap_path_key + + return cls diff --git a/tests/test_ecmwf_oper_data_fdb.py b/tests/test_ecmwf_oper_data_fdb.py index 553407ea..ebea9d99 100644 --- a/tests/test_ecmwf_oper_data_fdb.py +++ b/tests/test_ecmwf_oper_data_fdb.py @@ -16,7 +16,7 @@ def setup_method(self, method): "date": {"merge": {"with": "time", "linkers": ["T", "00"]}}, "step": {"type_change": "int"}, } - self.config = {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper", "type": "fc"} + self.config = {"class": "od", "expver": "0001", "levtype": "sfc", "type": "fc", "stream": "oper"} self.fdbdatacube = FDBDatacube(self.config, axis_options=self.options) self.slicer = HullSlicer() self.API = Polytope(datacube=self.fdbdatacube, engine=self.slicer, axis_options=self.options) diff --git a/tests/test_local_grid_cyclic.py b/tests/test_local_grid_cyclic.py new file mode 100644 index 00000000..02598632 --- /dev/null +++ b/tests/test_local_grid_cyclic.py @@ -0,0 +1,72 @@ +import pandas as pd +import pytest + +from polytope.engine.hullslicer import HullSlicer +from polytope.polytope import Polytope, Request +from polytope.shapes import Point, Select + + +class TestSlicingFDBDatacube: + def setup_method(self, method): + from polytope.datacube.backends.fdb import FDBDatacube + + # Create a dataarray with 3 labelled axes using different index types + self.options = { + "values": { + "mapper": { + "type": "local_regular", + "resolution": [80, 80], + "axes": ["latitude", "longitude"], + "local": [-40, 40, -20, 60], + } + }, + "date": {"merge": {"with": "time", "linkers": ["T", "00"]}}, + "step": {"type_change": "int"}, + "number": {"type_change": "int"}, + "longitude": {"cyclic": [-180, 180]}, + } + self.config = {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"} + self.fdbdatacube = FDBDatacube(self.config, axis_options=self.options) + self.slicer = HullSlicer() + self.API = Polytope(datacube=self.fdbdatacube, engine=self.slicer, axis_options=self.options) + + # Testing different shapes + @pytest.mark.fdb + def test_fdb_datacube(self): + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Select("date", [pd.Timestamp("20240129T000000")]), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["fc"]), + Point(["latitude", "longitude"], [[-20, -20]]), + ) + result = self.API.retrieve(request) + result.pprint_2() + assert len(result.leaves) == 1 + assert result.leaves[0].flatten()["latitude"] == -20 + assert result.leaves[0].flatten()["longitude"] == -20 + + @pytest.mark.fdb + def test_fdb_datacube_2(self): + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Select("date", [pd.Timestamp("20240129T000000")]), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["fc"]), + Point(["latitude", "longitude"], [[-20, 50 + 360]]), + ) + result = self.API.retrieve(request) + result.pprint_2() + assert len(result.leaves) == 1 + assert result.leaves[0].flatten()["latitude"] == -20 + assert result.leaves[0].flatten()["longitude"] == 50 diff --git a/tests/test_local_regular_grid.py b/tests/test_local_regular_grid.py new file mode 100644 index 00000000..dfac410e --- /dev/null +++ b/tests/test_local_regular_grid.py @@ -0,0 +1,227 @@ +import pandas as pd +import pytest + +from polytope.engine.hullslicer import HullSlicer +from polytope.polytope import Polytope, Request +from polytope.shapes import Point, Select + + +class TestSlicingFDBDatacube: + def setup_method(self, method): + from polytope.datacube.backends.fdb import FDBDatacube + + # Create a dataarray with 3 labelled axes using different index types + self.options = { + "values": { + "mapper": { + "type": "local_regular", + "resolution": [80, 80], + "axes": ["latitude", "longitude"], + "local": [-40, 40, -20, 60], + } + }, + "date": {"merge": {"with": "time", "linkers": ["T", "00"]}}, + "step": {"type_change": "int"}, + "number": {"type_change": "int"}, + } + self.config = {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"} + self.fdbdatacube = FDBDatacube(self.config, axis_options=self.options) + self.slicer = HullSlicer() + self.API = Polytope(datacube=self.fdbdatacube, engine=self.slicer, axis_options=self.options) + + # Testing different shapes + @pytest.mark.fdb + def test_fdb_datacube(self): + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Select("date", [pd.Timestamp("20240129T000000")]), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["fc"]), + Point(["latitude", "longitude"], [[0.16, 0.176]], method="nearest"), + ) + result = self.API.retrieve(request) + result.pprint_2() + assert len(result.leaves) == 1 + assert result.leaves[0].flatten()["latitude"] == 0 + assert result.leaves[0].flatten()["longitude"] == 0 + + @pytest.mark.fdb + def test_point_outside_local_region(self): + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Select("date", [pd.Timestamp("20240129T000000")]), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["fc"]), + Point(["latitude", "longitude"], [[0.16, 61]], method="nearest"), + ) + result = self.API.retrieve(request) + result.pprint_2() + assert len(result.leaves) == 1 + assert result.leaves[0].flatten()["latitude"] == 0 + assert result.leaves[0].flatten()["longitude"] == 60 + + @pytest.mark.fdb + def test_point_outside_local_region_2(self): + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Select("date", [pd.Timestamp("20240129T000000")]), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["fc"]), + Point(["latitude", "longitude"], [[41, 1]], method="nearest"), + ) + result = self.API.retrieve(request) + result.pprint_2() + assert len(result.leaves) == 1 + assert result.leaves[0].flatten()["latitude"] == 40 + assert result.leaves[0].flatten()["longitude"] == 1 + + @pytest.mark.fdb + def test_point_outside_local_region_3(self): + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Select("date", [pd.Timestamp("20240129T000000")]), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["fc"]), + Point(["latitude", "longitude"], [[1, 61]]), + ) + result = self.API.retrieve(request) + result.pprint_2() + assert len(result.leaves) == 1 + assert result.is_root() + + @pytest.mark.fdb + def test_point_outside_local_region_4(self): + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Select("date", [pd.Timestamp("20240129T000000")]), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["fc"]), + Point(["latitude", "longitude"], [[41, 1]]), + ) + result = self.API.retrieve(request) + result.pprint_2() + assert len(result.leaves) == 1 + assert result.is_root() + + @pytest.mark.fdb + def test_point_outside_local_region_5(self): + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Select("date", [pd.Timestamp("20240129T000000")]), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["fc"]), + Point(["latitude", "longitude"], [[-41, 1]]), + ) + result = self.API.retrieve(request) + result.pprint_2() + assert len(result.leaves) == 1 + assert result.is_root() + + @pytest.mark.fdb + def test_point_outside_local_region_6(self): + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Select("date", [pd.Timestamp("20240129T000000")]), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["fc"]), + Point(["latitude", "longitude"], [[-30, -21]]), + ) + result = self.API.retrieve(request) + result.pprint_2() + assert len(result.leaves) == 1 + assert result.is_root() + + @pytest.mark.fdb + def test_point_outside_local_region_7(self): + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Select("date", [pd.Timestamp("20240129T000000")]), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["fc"]), + Point(["latitude", "longitude"], [[-41, 1]], method="nearest"), + ) + result = self.API.retrieve(request) + result.pprint_2() + assert len(result.leaves) == 1 + assert result.leaves[0].flatten()["latitude"] == -40 + assert result.leaves[0].flatten()["longitude"] == 1 + + @pytest.mark.fdb + def test_point_outside_local_region_8(self): + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Select("date", [pd.Timestamp("20240129T000000")]), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["fc"]), + Point(["latitude", "longitude"], [[-30, -21]], method="nearest"), + ) + result = self.API.retrieve(request) + result.pprint_2() + assert len(result.leaves) == 1 + assert result.leaves[0].flatten()["latitude"] == -30 + assert result.leaves[0].flatten()["longitude"] == -20 + + @pytest.mark.fdb + def test_point_outside_local_region_9(self): + request = Request( + Select("step", [0]), + Select("levtype", ["sfc"]), + Select("date", [pd.Timestamp("20240129T000000")]), + Select("domain", ["g"]), + Select("expver", ["0001"]), + Select("param", ["167"]), + Select("class", ["od"]), + Select("stream", ["oper"]), + Select("type", ["fc"]), + Point(["latitude", "longitude"], [[-30, -21]], method="surrounding"), + ) + result = self.API.retrieve(request) + result.pprint_2() + assert len(result.leaves) == 3 + assert result.leaves[0].flatten()["latitude"] == -31 + assert result.leaves[0].flatten()["longitude"] == -20 diff --git a/tests/test_mappers.py b/tests/test_mappers.py index b9a39477..1231f19d 100644 --- a/tests/test_mappers.py +++ b/tests/test_mappers.py @@ -1,4 +1,6 @@ -from polytope.datacube.transformations.datacube_mappers import OctahedralGridMapper +from polytope.datacube.transformations.datacube_mappers.mapper_types.octahedral import ( + OctahedralGridMapper, +) class TestMapper: diff --git a/tests/test_point_nearest.py b/tests/test_point_nearest.py index 081050da..fa53b980 100644 --- a/tests/test_point_nearest.py +++ b/tests/test_point_nearest.py @@ -95,7 +95,7 @@ def test_fdb_datacube_true_point_3(self): result = self.API.retrieve(request) result.pprint() assert len(result.leaves) == 1 - assert result.leaves[0].value == 0 + assert result.leaves[0].value == 359.929906542056 assert result.leaves[0].axis.name == "longitude" @pytest.mark.fdb diff --git a/tests/test_slice_date_range_fdb.py b/tests/test_slice_date_range_fdb.py index ec14ba48..319bc2c2 100644 --- a/tests/test_slice_date_range_fdb.py +++ b/tests/test_slice_date_range_fdb.py @@ -17,7 +17,7 @@ def setup_method(self, method): "step": {"type_change": "int"}, "number": {"type_change": "int"}, } - self.config = {"class": "od", "expver": "0001", "levtype": "sfc", "step": "0"} + self.config = {"class": "od", "expver": "0001", "levtype": "sfc", "stream": "oper"} self.fdbdatacube = FDBDatacube(self.config, axis_options=self.options) self.slicer = HullSlicer() self.API = Polytope(datacube=self.fdbdatacube, engine=self.slicer, axis_options=self.options) @@ -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"]),