diff --git a/semeio/workflows/localisation/local_script_lib.py b/semeio/workflows/localisation/local_script_lib.py index ffa6c715a..2c4a2aab2 100644 --- a/semeio/workflows/localisation/local_script_lib.py +++ b/semeio/workflows/localisation/local_script_lib.py @@ -92,6 +92,218 @@ def from_list(cls, input_list): return cls([Parameter(key, val) for key, val in result.items()]) +@dataclass +class Decay: + obs_pos: list + main_range: float + perp_range: float + azimuth: float + grid: object + + def __post_init__(self): + angle = (90.0 - self.azimuth) * math.pi / 180.0 + self.cosangle = math.cos(angle) + self.sinangle = math.sin(angle) + + def get_dx_dy(self, data_index): + try: + # Assume the grid is 3D Grid + x, y, _ = self.grid.get_xyz(active_index=data_index) + except AttributeError: + # Assume the grid is a 2D Surface grid + # pylint: disable=no-member + x, y = self.grid.getXY(data_index) + x_unrotated = x - self.obs_pos[0] + y_unrotated = y - self.obs_pos[1] + + dx = ( + x_unrotated * self.cosangle + y_unrotated * self.sinangle + ) / self.main_range + dy = ( + -x_unrotated * self.sinangle + y_unrotated * self.cosangle + ) / self.perp_range + return dx, dy + + def norm_dist_square(self, data_index): + dx, dy = self.get_dx_dy(data_index) + d2 = dx**2 + dy**2 + return d2 + + +@dataclass +class GaussianDecay(Decay): + cutoff: bool + + def __call__(self, data_index): + d2 = super().norm_dist_square(data_index) + if self.cutoff and d2 > 1.0: + return 0.0 + exp_arg = -3.0 * d2 + return math.exp(exp_arg) + + +@dataclass +class ConstGaussianDecay(Decay): + normalised_tapering_range: float + cutoff: bool + + def __call__(self, data_index): + d2 = super().norm_dist_square(data_index) + d = math.sqrt(d2) + if d <= 1.0: + return 1.0 + if self.cutoff and d > self.normalised_tapering_range: + return 0.0 + + distance_from_inner_ellipse = (d - 1) / (self.normalised_tapering_range - 1) + exp_arg = -3 * distance_from_inner_ellipse**2 + return math.exp(exp_arg) + + +@dataclass +class ExponentialDecay(Decay): + cutoff: bool + + def __call__(self, data_index): + d2 = super().norm_dist_square(data_index) + d = math.sqrt(d2) + if self.cutoff and d > 1.0: + return 0.0 + exp_arg = -3.0 * d + return math.exp(exp_arg) + + +@dataclass +class ConstExponentialDecay(Decay): + normalised_tapering_range: float + cutoff: bool + + def __call__(self, data_index): + d2 = super().norm_dist_square(data_index) + d = math.sqrt(d2) + if d <= 1.0: + return 1.0 + if self.cutoff and d > self.normalised_tapering_range: + return 0.0 + + distance_from_inner_ellipse = (d - 1) / (self.normalised_tapering_range - 1) + exp_arg = -3 * distance_from_inner_ellipse + return math.exp(exp_arg) + + +@dataclass +class ConstantScalingFactor: + value: float + + def __call__(self, _): + return self.value + + +class ScalingValues: + scaling_param_number = 1 + corr_name = None + + @classmethod + def initialize(cls): + cls.scaling_param_number = 1 + cls.corr_name = None + + @classmethod + def write_qc_parameter_field( + cls, + node_name, + corr_name, + field_scale, + grid, + param_for_field, + log_level=LogLevel.OFF, + ): + # pylint: disable=too-many-arguments + + if param_for_field is None or field_scale is None: + return + + # Write scaling parameter once per corr_name + if corr_name == cls.corr_name: + return + + cls.corr_name = corr_name + + # Need a parameter name <= 8 character long for GRDECL format + scaling_kw_name = "S_" + str(cls.scaling_param_number) + filename = cls.corr_name + "_" + node_name + "_" + scaling_kw_name + ".roff" + + scaling_values = np.reshape( + param_for_field, (grid.getNX(), grid.getNY(), grid.getNZ()), "C" + ) + save_field(scaling_values, scaling_kw_name, filename, FieldFileFormat.ROFF) + + print( + "Write calculated scaling factor with name: " + f"{scaling_kw_name} to file: {filename}" + ) + debug_print( + f"Write calculated scaling factor with name: " + f"{scaling_kw_name} to file: {filename}", + LogLevel.LEVEL3, + log_level, + ) + + # Increase parameter number to define unique parameter name + cls.scaling_param_number = cls.scaling_param_number + 1 + + @classmethod + def write_qc_parameter_surface( + # pylint: disable=too-many-locals + cls, + node_name: str, + corr_name: str, + surface_scale: bool, + reference_surface_file: str, + param_for_surface: object, + log_level=LogLevel.OFF, + ): + # pylint: disable=too-many-arguments + + if param_for_surface is None or surface_scale is None: + return + + # Write scaling parameter once per corr_name + if corr_name == cls.corr_name: + return + + cls.corr_name = corr_name + + scaling_surface_name = str(cls.scaling_param_number) + filename = cls.corr_name + "_" + node_name + "_map.irap" + + # Surface object is created from the known reference surface to get + # surface attributes. The param_for_surface array is c-order, + # but need f-order for surface to write with the libecl object + qc_surface = Surface(reference_surface_file) + nx = qc_surface.getNX() + ny = qc_surface.getNY() + for i in range(nx): + for j in range(ny): + c_indx = j + i * ny + f_indx = i + j * nx + qc_surface[f_indx] = param_for_surface[c_indx] + qc_surface.write(filename) + print( + "Write calculated scaling factor for SURFACE parameter " + f"in {cls.corr_name} to file: {filename}" + ) + debug_print( + f"Write calculated scaling factor with name: " + f"{scaling_surface_name} to file: {filename}", + LogLevel.LEVEL3, + log_level, + ) + + # Increase parameter number to define unique parameter name + cls.scaling_param_number = cls.scaling_param_number + 1 + + def get_param_from_ert(ens_config): new_params = Parameters() for key in ens_config.parameters: @@ -152,8 +364,8 @@ def build_decay_object( azimuth: float, grid: object, use_cutoff: bool, - tapering_range: float = None, -): + tapering_range: float = 1.5, +) -> Decay: # pylint: disable=too-many-arguments if method == "gaussian_decay": decay_obj = GaussianDecay( @@ -906,215 +1118,3 @@ def check_if_ref_point_in_grid(ref_point, grid, log_level): LogLevel.LEVEL3, log_level, ) - - -@dataclass -class Decay: - obs_pos: list - main_range: float - perp_range: float - azimuth: float - grid: object - - def __post_init__(self): - angle = (90.0 - self.azimuth) * math.pi / 180.0 - self.cosangle = math.cos(angle) - self.sinangle = math.sin(angle) - - def get_dx_dy(self, data_index): - try: - # Assume the grid is 3D Grid - x, y, _ = self.grid.get_xyz(active_index=data_index) - except AttributeError: - # Assume the grid is a 2D Surface grid - # pylint: disable=no-member - x, y = self.grid.getXY(data_index) - x_unrotated = x - self.obs_pos[0] - y_unrotated = y - self.obs_pos[1] - - dx = ( - x_unrotated * self.cosangle + y_unrotated * self.sinangle - ) / self.main_range - dy = ( - -x_unrotated * self.sinangle + y_unrotated * self.cosangle - ) / self.perp_range - return dx, dy - - def norm_dist_square(self, data_index): - dx, dy = self.get_dx_dy(data_index) - d2 = dx**2 + dy**2 - return d2 - - -@dataclass -class GaussianDecay(Decay): - cutoff: bool - - def __call__(self, data_index): - d2 = super().norm_dist_square(data_index) - if self.cutoff and d2 > 1.0: - return 0.0 - exp_arg = -3.0 * d2 - return math.exp(exp_arg) - - -@dataclass -class ConstGaussianDecay(Decay): - normalised_tapering_range: float - cutoff: bool - - def __call__(self, data_index): - d2 = super().norm_dist_square(data_index) - d = math.sqrt(d2) - if d <= 1.0: - return 1.0 - if self.cutoff and d > self.normalised_tapering_range: - return 0.0 - - distance_from_inner_ellipse = (d - 1) / (self.normalised_tapering_range - 1) - exp_arg = -3 * distance_from_inner_ellipse**2 - return math.exp(exp_arg) - - -@dataclass -class ExponentialDecay(Decay): - cutoff: bool - - def __call__(self, data_index): - d2 = super().norm_dist_square(data_index) - d = math.sqrt(d2) - if self.cutoff and d > 1.0: - return 0.0 - exp_arg = -3.0 * d - return math.exp(exp_arg) - - -@dataclass -class ConstExponentialDecay(Decay): - normalised_tapering_range: float - cutoff: bool - - def __call__(self, data_index): - d2 = super().norm_dist_square(data_index) - d = math.sqrt(d2) - if d <= 1.0: - return 1.0 - if self.cutoff and d > self.normalised_tapering_range: - return 0.0 - - distance_from_inner_ellipse = (d - 1) / (self.normalised_tapering_range - 1) - exp_arg = -3 * distance_from_inner_ellipse - return math.exp(exp_arg) - - -@dataclass -class ConstantScalingFactor: - value: float - - def __call__(self, _): - return self.value - - -class ScalingValues: - scaling_param_number = 1 - corr_name = None - - @classmethod - def initialize(cls): - cls.scaling_param_number = 1 - cls.corr_name = None - - @classmethod - def write_qc_parameter_field( - cls, - node_name, - corr_name, - field_scale, - grid, - param_for_field, - log_level=LogLevel.OFF, - ): - # pylint: disable=too-many-arguments - - if param_for_field is None or field_scale is None: - return - - # Write scaling parameter once per corr_name - if corr_name == cls.corr_name: - return - - cls.corr_name = corr_name - - # Need a parameter name <= 8 character long for GRDECL format - scaling_kw_name = "S_" + str(cls.scaling_param_number) - filename = cls.corr_name + "_" + node_name + "_" + scaling_kw_name + ".roff" - - scaling_values = np.reshape( - param_for_field, (grid.getNX(), grid.getNY(), grid.getNZ()), "C" - ) - save_field(scaling_values, scaling_kw_name, filename, FieldFileFormat.ROFF) - - print( - "Write calculated scaling factor with name: " - f"{scaling_kw_name} to file: {filename}" - ) - debug_print( - f"Write calculated scaling factor with name: " - f"{scaling_kw_name} to file: {filename}", - LogLevel.LEVEL3, - log_level, - ) - - # Increase parameter number to define unique parameter name - cls.scaling_param_number = cls.scaling_param_number + 1 - - @classmethod - def write_qc_parameter_surface( - # pylint: disable=too-many-locals - cls, - node_name: str, - corr_name: str, - surface_scale: bool, - reference_surface_file: str, - param_for_surface: object, - log_level=LogLevel.OFF, - ): - # pylint: disable=too-many-arguments - - if param_for_surface is None or surface_scale is None: - return - - # Write scaling parameter once per corr_name - if corr_name == cls.corr_name: - return - - cls.corr_name = corr_name - - scaling_surface_name = str(cls.scaling_param_number) - filename = cls.corr_name + "_" + node_name + "_map.irap" - - # Surface object is created from the known reference surface to get - # surface attributes. The param_for_surface array is c-order, - # but need f-order for surface to write with the libecl object - qc_surface = Surface(reference_surface_file) - nx = qc_surface.getNX() - ny = qc_surface.getNY() - for i in range(nx): - for j in range(ny): - c_indx = j + i * ny - f_indx = i + j * nx - qc_surface[f_indx] = param_for_surface[c_indx] - qc_surface.write(filename) - print( - "Write calculated scaling factor for SURFACE parameter " - f"in {cls.corr_name} to file: {filename}" - ) - debug_print( - f"Write calculated scaling factor with name: " - f"{scaling_surface_name} to file: {filename}", - LogLevel.LEVEL3, - log_level, - ) - - # Increase parameter number to define unique parameter name - cls.scaling_param_number = cls.scaling_param_number + 1