diff --git a/.DS_Store b/.DS_Store index 72b7758..2455104 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/plato/__pycache__/globe.cpython-310.pyc b/plato/__pycache__/globe.cpython-310.pyc index 46271f2..6c5f878 100644 Binary files a/plato/__pycache__/globe.cpython-310.pyc and b/plato/__pycache__/globe.cpython-310.pyc differ diff --git a/plato/__pycache__/grids.cpython-310.pyc b/plato/__pycache__/grids.cpython-310.pyc index 35ae369..5a93000 100644 Binary files a/plato/__pycache__/grids.cpython-310.pyc and b/plato/__pycache__/grids.cpython-310.pyc differ diff --git a/plato/__pycache__/optimisation.cpython-310.pyc b/plato/__pycache__/optimisation.cpython-310.pyc index 314b3dc..e96325d 100644 Binary files a/plato/__pycache__/optimisation.cpython-310.pyc and b/plato/__pycache__/optimisation.cpython-310.pyc differ diff --git a/plato/__pycache__/plate_torques.cpython-310.pyc b/plato/__pycache__/plate_torques.cpython-310.pyc index cea6266..cb3eb34 100644 Binary files a/plato/__pycache__/plate_torques.cpython-310.pyc and b/plato/__pycache__/plate_torques.cpython-310.pyc differ diff --git a/plato/__pycache__/plates.cpython-310.pyc b/plato/__pycache__/plates.cpython-310.pyc index 479e835..92ed6ea 100644 Binary files a/plato/__pycache__/plates.cpython-310.pyc and b/plato/__pycache__/plates.cpython-310.pyc differ diff --git a/plato/__pycache__/points.cpython-310.pyc b/plato/__pycache__/points.cpython-310.pyc index d2b7702..52ccd31 100644 Binary files a/plato/__pycache__/points.cpython-310.pyc and b/plato/__pycache__/points.cpython-310.pyc differ diff --git a/plato/__pycache__/settings.cpython-310.pyc b/plato/__pycache__/settings.cpython-310.pyc index cfd2112..0e55904 100644 Binary files a/plato/__pycache__/settings.cpython-310.pyc and b/plato/__pycache__/settings.cpython-310.pyc differ diff --git a/plato/__pycache__/slabs.cpython-310.pyc b/plato/__pycache__/slabs.cpython-310.pyc index 2f06cc9..2a43d9e 100644 Binary files a/plato/__pycache__/slabs.cpython-310.pyc and b/plato/__pycache__/slabs.cpython-310.pyc differ diff --git a/plato/__pycache__/utils_calc.cpython-310.pyc b/plato/__pycache__/utils_calc.cpython-310.pyc index 02f95b4..b416d53 100644 Binary files a/plato/__pycache__/utils_calc.cpython-310.pyc and b/plato/__pycache__/utils_calc.cpython-310.pyc differ diff --git a/plato/__pycache__/utils_data.cpython-310.pyc b/plato/__pycache__/utils_data.cpython-310.pyc index ff44cbb..22ea5a2 100644 Binary files a/plato/__pycache__/utils_data.cpython-310.pyc and b/plato/__pycache__/utils_data.cpython-310.pyc differ diff --git a/plato/globe.py b/plato/globe.py index 5dd7623..9c68920 100644 --- a/plato/globe.py +++ b/plato/globe.py @@ -187,6 +187,7 @@ def calculate_net_rotation( ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Calculate the net rotation of the Earth's lithosphere. @@ -198,7 +199,11 @@ def calculate_net_rotation( _cases = utils_data.select_cases(cases, self.settings.cases) # Calculate the net rotation of the Earth's lithosphere - for i, _age in enumerate(_ages): + for i, _age in enumerate(_tqdm( + _ages, + desc="Calculating net rotation", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR), + )): for _case in _cases: # Check if plate data is provided if utils_init.check_object_data(plates, Plates, _age, _case): @@ -217,10 +222,10 @@ def calculate_net_rotation( ) # Check if plate data is provided - if utils_init.check_object_data(plates, Plates, _age, _case): + if utils_init.check_object_data(points, Points, _age, _case): logging.info(f"Calculating net rotation for case {_case} at age {_age} using provided Points data") _points = points - elif utils_init.check_object_data(self.plates, Plates, _age, _case): + elif utils_init.check_object_data(self.points, Points, _age, _case): logging.info(f"Calculating net rotation for case {_case} at age {_age} using stored Points data") _points = self.points else: @@ -239,10 +244,17 @@ def calculate_net_rotation( # Select plates and points data selected_plates = _plates.data[_age][_case] selected_points = _points.data[_age][_case] + + # Filter by plateIDs if plateIDs is not None: selected_plates = _plates.data[_age][_case][_plates.data[_age][_case].plateID.isin(_plateIDs)] selected_points = _points.data[_age][_case][_points.data[_age][_case].plateID.isin(_plateIDs)] + # Filter by minimum plate area + if self.settings.options[_case]["Minimum plate area"] > 0.: + selected_plates = selected_plates[selected_plates.area >= self.settings.options[_case]["Minimum plate area"]] + selected_points = selected_points[selected_points.plateID.isin(selected_plates.plateID)] + # Calculate net rotation net_rotation_pole = utils_calc.compute_net_rotation( selected_plates, @@ -311,6 +323,7 @@ def save( self, cases: Optional[Union[str, List[str]]] = None, file_dir: Optional[str] = None, + PROGRESS_BAR: Optional[bool] = True, ): """ Function to save 'Globe' object. @@ -323,7 +336,11 @@ def save( _file_dir = self.settings.dir_path if file_dir is None else file_dir # Loop through ages - for _case in _tqdm(_cases, desc="Saving Globe", disable=self.settings.logger.level==logging.INFO): + for _case in _tqdm( + _cases, + desc="Saving Globe", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): utils_data.DataFrame_to_parquet( self.data[_case], "Globe", @@ -337,6 +354,7 @@ def export( self, cases: Optional[Union[str, List[str]]] = None, file_dir: Optional[str] = None, + PROGRESS_BAR: Optional[bool] = True, ): """ Function to export 'Globe' object. @@ -349,7 +367,11 @@ def export( _file_dir = self.settings.dir_path if file_dir is None else file_dir # Loop through ages - for _case in _tqdm(_cases, desc="Saving Globe", disable=self.settings.logger.level==logging.INFO): + for _case in _tqdm( + _cases, + desc="Exporting Globe", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): utils_data.DataFrame_to_csv( self.data[_case], "Globe", diff --git a/plato/grids.py b/plato/grids.py index 664d91d..13b75f2 100644 --- a/plato/grids.py +++ b/plato/grids.py @@ -353,73 +353,78 @@ def interpolate_data_to_grid( logging.info(f"{grid_type} updated:", getattr(self, grid_type)) def save_all( - self, - ages: Union[None, List[int], List[float], _numpy.ndarray] = None, - cases: Union[None, str, List[str]] = None, - file_dir: Optional[str] = None, + self, + ages: Union[None, List[int], List[float], _numpy.ndarray] = None, + cases: Union[None, str, List[str]] = None, + file_dir: Optional[str] = None, + PROGRESS_BAR: bool = True, ): """ Function to save all the grids """ # Save seafloor grid - self.save_seafloor_age(ages, file_dir) + self.save_seafloor_age(ages, file_dir, PROGRESS_BAR) # Save sediment grid - self.save_sediment(ages, cases, file_dir) + self.save_sediment(ages, cases, file_dir, PROGRESS_BAR) # Save continental grid - self.save_continent(ages, cases, file_dir) + self.save_continent(ages, cases, file_dir, PROGRESS_BAR) # Save velocity grid - self.save_velocity(ages, cases, file_dir) + self.save_velocity(ages, cases, file_dir, PROGRESS_BAR) def save_seafloor_age( self, ages: Union[None, List[int], List[float], _numpy.ndarray] = None, file_dir: Optional[str] = None, + PROGRESS_BAR: bool = True, ): """ Function to save the the seafloor age grid. """ - self.save_grid(self.seafloor_age, "Seafloor_age", ages, None, file_dir) + self.save_grid(self.seafloor_age, "Seafloor_age", ages, None, file_dir, PROGRESS_BAR) def save_sediment( self, ages: Union[None, List[int], List[float], _numpy.ndarray] = None, cases: Union[None, str, List[str]] = None, file_dir: Optional[str] = None, + PROGRESS_BAR: bool = True, ): """ Function to save the the sediment grid. """ if self.sediment is not None: - self.save_grid(self.sediment, "Sediment", ages, cases, file_dir) + self.save_grid(self.sediment, "Sediment", ages, cases, file_dir, PROGRESS_BAR) def save_continent( self, ages: Union[None, List[int], List[float], _numpy.ndarray] = None, cases: Union[None, str, List[str]] = None, file_dir: Optional[str] = None, + PROGRESS_BAR: bool = True, ): """ Function to save the the continental grid. """ # Check if grids exists if self.continent is not None: - self.save_grid(self.continent, "Continent", ages, cases, file_dir) + self.save_grid(self.continent, "Continent", ages, cases, file_dir, PROGRESS_BAR) def save_velocity( self, ages: Union[None, List[int], List[float], _numpy.ndarray] = None, cases: Union[None, str, List[str]] = None, file_dir: Optional[str] = None, + PROGRESS_BAR: bool = True, ): """ Function to save the the velocity grid. """ # Check if grids exists if self.velocity is not None: - self.save_grid(self.velocity, "Velocity", ages, cases, file_dir) + self.save_grid(self.velocity, "Velocity", ages, cases, file_dir, PROGRESS_BAR) def save_grid( self, @@ -428,6 +433,7 @@ def save_grid( ages: Union[None, List[int], List[float], _numpy.ndarray] = None, cases: Union[None, str, List[str]] = None, file_dir: Optional[str] = None, + PROGRESS_BAR: bool = True, ): """ Function to save a grid @@ -442,7 +448,11 @@ def save_grid( _file_dir = self.settings.dir_path if file_dir is None else file_dir # Loop through ages - for _age in _tqdm(_ages, desc=f"Saving {type} grids", disable=self.settings.logger.level==logging.INFO): + for _age in _tqdm( + _ages, + desc=f"Saving {type} grids", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): if data[_age] is Dict: # Loop through cases for _case in _cases: diff --git a/plato/optimisation.py b/plato/optimisation.py index 67838b4..167bc12 100644 --- a/plato/optimisation.py +++ b/plato/optimisation.py @@ -197,7 +197,7 @@ def minimise_residual_torque( def optimise_slab_pull_coefficient( self, age = None, - case = None, + cases = None, plateIDs = None, grid_size = 500, viscosity = 1.23e20, @@ -224,7 +224,7 @@ def optimise_slab_pull_coefficient( _ages = utils_data.select_ages(age, self.settings.ages) # Select cases, if not provided - _cases = utils_data.select_cases(case, self.settings.reconstructed_cases) + _cases = utils_data.select_cases(cases, self.settings.reconstructed_cases) for _age in _tqdm(_ages, desc="Optimising slab pull coefficient"): for _case in _cases: @@ -297,15 +297,21 @@ def optimise_slab_pull_coefficient( mask = self.slabs.data[_age][_case]["lower_plateID"] == _plateID self.slabs.data[_age][_case].loc[mask, "slab_pull_constant"] = opt_sp_const # self.opt_sp_const[_age][_case][_data.index[k]] = opt_sp_const + + # Recalculate all the relevant torques + self.plate_torques.calculate_slab_pull_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs, PROGRESS_BAR=False) + self.plate_torques.calculate_driving_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs, PROGRESS_BAR=False) + self.plate_torques.calculate_residual_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs, PROGRESS_BAR=False) def invert_residual_torque( self, age = None, - case = None, + cases = None, plateIDs = None, parameter = "Slab pull constant", - grid_size = 500, - viscosity = 1.23e20, + vmin = 8.5, + vmax = 13.5, + step = 0.1, plot = False, ): """ @@ -333,45 +339,51 @@ def invert_residual_torque( _ages = utils_data.select_ages(age, self.settings.ages) # Select cases, if not provided - _cases = utils_data.select_cases(case, self.settings.reconstructed_cases) + _cases = utils_data.select_cases(cases, self.settings.reconstructed_cases) # Define plateIDs, if not provided _plateIDs = utils_data.select_plateIDs( plateIDs, - [101, # North America - 201, # South America - 301, # Eurasia - 501, # India - 801, # Australia - 802, # Antarctica - 901, # Pacific - 902, # Farallon - 911, # Nazca - 919, # Phoenix - 926, # Izanagi + [ + 101, # North America + 201, # South America + 301, # Eurasia + 501, # India + 801, # Australia + 802, # Antarctica + 901, # Pacific + 902, # Farallon + 911, # Nazca + 919, # Phoenix + 926, # Izanagi ] ) # Define constants - constants = _numpy.arange(8.5, 13.5, .1) + constants = _numpy.arange(vmin, vmax, step) driving_torque_opt_stack = {_case: {_plateID: _numpy.zeros((len(constants), len(_ages))) for _plateID in _plateIDs} for _case in _cases} residual_torque_opt_stack = {_case: {_plateID: _numpy.zeros((len(constants), len(_ages))) for _plateID in _plateIDs} for _case in _cases} - for i, constant in _tqdm(enumerate(constants), desc="Optimising slab pull coefficient"): + for i, constant in enumerate(_tqdm(constants, desc="Inverting residual torque to determine optimal slab pull coefficient")): logging.info(f"Optimising for {constant}") _data = {} # Calculate the torques the normal way - self.plate_torques.calculate_slab_pull_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs) - self.plate_torques.calculate_driving_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs) - self.plate_torques.calculate_residual_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs) + self.plate_torques.calculate_slab_pull_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs, PROGRESS_BAR=False) + self.plate_torques.calculate_driving_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs, PROGRESS_BAR=False) + self.plate_torques.calculate_residual_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs, PROGRESS_BAR=False) # Loop through ages for _age in _ages: - for _case in _cases: + if _age not in _data: _data[_age] = {} - _data[_age][_case] = self.plate_torques.slabs.data[_age][_case].copy() + + for _case in _cases: + _data[_age][_case] = self.slabs.data[_age][_case].copy() + + # Get the slab pull force magnitude + max_slab_pull_force_mag = _data[_age][_case]["slab_pull_force_mag"] / _data[_age][_case]["slab_pull_constant"] # Modify the magnitude of the slab pull force using the 2D dot product of the residual force and the slab pull force and the constant _data[_age][_case]["slab_pull_force_mag"] -= ( @@ -379,14 +391,18 @@ def invert_residual_torque( _data[_age][_case]["residual_force_lon"] * _data[_age][_case]["slab_pull_force_lon"] ) * 10**-constant + # Ensure the slab pull force magnitude is positive and not larger than the original slab pull force magnitude + _data[_age][_case].loc[_data[_age][_case]["slab_pull_force_mag"] < 0, "slab_pull_force_mag"] = 0 + _data[_age][_case].loc[_data[_age][_case]["slab_pull_force_mag"] > max_slab_pull_force_mag, "slab_pull_force_mag"] = max_slab_pull_force_mag[_data[_age][_case]["slab_pull_force_mag"] > max_slab_pull_force_mag] + # Decompose the slab pull force into latitudinal and longitudinal components using the trench normal azimuth _data[_age][_case]["slab_pull_force_lat"] = _numpy.cos(_numpy.deg2rad(_data[_age][_case]["trench_normal_azimuth"])) * _data[_age][_case]["slab_pull_force_mag"] _data[_age][_case]["slab_pull_force_lon"] = _numpy.sin(_numpy.deg2rad(_data[_age][_case]["trench_normal_azimuth"])) * _data[_age][_case]["slab_pull_force_mag"] # Calculate the torques with the modified slab pull forces - self.plate_torques.plates.calculate_torque_on_plates(_data, ages=_ages, cases=_cases, plateIDs=_plateIDs, torque_var="slab_pull", PROGRESS_BAR = False) - self.plate_torques.calculate_driving_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs, PROGRESS_BAR = False) - self.plate_torques.calculate_residual_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs, PROGRESS_BAR = False) + self.plates.calculate_torque_on_plates(_data, ages=_ages, cases=_cases, plateIDs=_plateIDs, torque_var="slab_pull", PROGRESS_BAR=False) + self.plate_torques.calculate_driving_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs, PROGRESS_BAR=False) + self.plate_torques.calculate_residual_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs, PROGRESS_BAR=False) # Extract the driving and residual torques _iter_driving_torque = self.plate_torques.extract_data_through_time(ages=_ages, cases=_cases, plateIDs=_plateIDs, var="driving_torque_mag") @@ -394,33 +410,69 @@ def invert_residual_torque( for _case in _cases: for _plateID in _plateIDs: - driving_torque_opt_stack[_case][_plateID][i, :] = _iter_driving_torque[_case][_plateID].values - residual_torque_opt_stack[_case][_plateID][i, :] = _iter_residual_torque[_case][_plateID].values + if len(_cases) == 1: + driving_torque_opt_stack[_case][_plateID][i, :] = _iter_driving_torque[_plateID].values + residual_torque_opt_stack[_case][_plateID][i, :] = _iter_residual_torque[_plateID].values + else: + driving_torque_opt_stack[_case][_plateID][i, :] = _iter_driving_torque[_case][_plateID].values + residual_torque_opt_stack[_case][_plateID][i, :] = _iter_residual_torque[_case][_plateID].values # Get minimum residual torque for each unique combination of case and plate ID. minimum_residual_torque = {_case: {_plateID: None for _plateID in _plateIDs} for _case in _cases} - opt_constant = {_case: {_plateID: None for _plateID in _plateIDs} for _case in _cases} + opt_constants = {_case: {_plateID: None for _plateID in _plateIDs} for _case in _cases} for _case in _cases: for _plateID in _plateIDs: - minimum_residual_torque[_case][_plateID] = _numpy.min(residual_torque_opt_stack[_plateID]/driving_torque_opt_stack[_plateID], axis = 0) - opt_index = _numpy.argmin(residual_torque_opt_stack[_plateID]/driving_torque_opt_stack[_plateID], axis = 0) - opt_constant[_case][_plateID] = constants[opt_index] + minimum_residual_torque[_case][_plateID] = _numpy.min(residual_torque_opt_stack[_case][_plateID]/driving_torque_opt_stack[_case][_plateID], axis = 0) + opt_index = _numpy.argmin(residual_torque_opt_stack[_case][_plateID]/driving_torque_opt_stack[_case][_plateID], axis = 0) + opt_constants[_case][_plateID] = constants[opt_index] + + # Recalculate all the relevant torques + self.plate_torques.calculate_slab_pull_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs, PROGRESS_BAR=False) + self.plate_torques.calculate_residual_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs, PROGRESS_BAR=False) # Calculate optimal slab pull constant and store in slab data for k, _age in enumerate(_ages): for _case in _cases: for _plateID in _plateIDs: # Select data - _data = self.slabs.data[_age][_case] + _data = self.slabs.data[_age][_case].copy() _data = _data[_data["lower_plateID"] == _plateID] - # Calculate slab pull constant - _data["slab_pull_constant"] = ( - _data["slab_pull_force_mag"] - ( - _data[_age][_case]["residual_force_lat"] * _data[_age][_case]["slab_pull_force_lat"] + \ - _data[_age][_case]["residual_force_lon"] * _data[_age][_case]["slab_pull_force_lon"] - ) * 10**-constant - ) / _data["slab_pull_force_mag"] + if not _data.empty: + # Get the old slab pull force magnitude + max_slab_pull_force_mag = _data["slab_pull_force_mag"].values / _data["slab_pull_constant"].values + + # Calculate the slab pull force + _data.loc[_data.index, "slab_pull_force_mag"] -= ( + _data["residual_force_lat"] * _data["slab_pull_force_lat"] + \ + _data["residual_force_lon"] * _data["slab_pull_force_lon"] + ) * 10**-opt_constants[_case][_plateID][k] + + # Make sure the slab pull force magnitude is positive and not larger than the original slab pull force magnitude + _data.loc[_data["slab_pull_force_mag"] < 0, "slab_pull_force_mag"] = 0 + _data.loc[_data["slab_pull_force_mag"] > max_slab_pull_force_mag, "slab_pull_force_mag"] = max_slab_pull_force_mag[_data["slab_pull_force_mag"] > max_slab_pull_force_mag] + + # Decompose the slab pull force into latitudinal and longitudinal components using the trench normal azimuth + _data.loc[_data.index, "slab_pull_force_lat"] = _numpy.cos(_numpy.deg2rad(_data["trench_normal_azimuth"])) * _data["slab_pull_force_mag"] + _data.loc[_data.index, "slab_pull_force_lon"] = _numpy.sin(_numpy.deg2rad(_data["trench_normal_azimuth"])) * _data["slab_pull_force_mag"] + + # Calculate the slab pull constant + _data.loc[_data.index, "slab_pull_constant"] = _data.loc[_data.index, "slab_pull_force_mag"] / max_slab_pull_force_mag + + # Make sure the slab pull constant is between 0 and 1 + _data.loc[_data["slab_pull_constant"] < 0, "slab_pull_constant"] = 0 + _data.loc[_data["slab_pull_constant"] > 1, "slab_pull_constant"] = 1 + + # Feed optimal values back into slab data + self.slabs.data[_age][_case].loc[_data.index, "slab_pull_force_mag"] = _data["slab_pull_force_mag"].values + self.slabs.data[_age][_case].loc[_data.index, "slab_pull_force_lat"] = _data["slab_pull_force_lat"].values + self.slabs.data[_age][_case].loc[_data.index, "slab_pull_force_lon"] = _data["slab_pull_force_lon"].values + self.slabs.data[_age][_case].loc[_data.index, "slab_pull_constant"] = _data["slab_pull_constant"].values + + # Recalculate all the relevant torques + self.plates.calculate_torque_on_plates(self.slabs.data, ages=_ages, cases=_cases, plateIDs=_plateIDs, torque_var="slab_pull", PROGRESS_BAR=False) + self.plate_torques.calculate_driving_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs, PROGRESS_BAR=False) + self.plate_torques.calculate_residual_torque(ages=_ages, cases=_cases, plateIDs=_plateIDs, PROGRESS_BAR=False) def optimise_torques( self, diff --git a/plato/plate_torques.py b/plato/plate_torques.py index d478e5a..58fcf69 100644 --- a/plato/plate_torques.py +++ b/plato/plate_torques.py @@ -79,7 +79,7 @@ def __init__( grids = None, globe = None, DEBUG_MODE = False, - PARALLEL_MODE = False + PARALLEL_MODE = False, ): """ Constructor for the PlateTorques class. @@ -181,7 +181,7 @@ def add_grid( grid_type = "seafloor_age", target_variable = "z", mask_continents = False, - prefactor = 1. + prefactor = 1., ): """ Function to add a grid to the grids object. This calls the add_grid method in the Grids class. @@ -214,6 +214,7 @@ def calculate_rms_velocity( ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to calculate root mean square velocities for plates. @@ -231,6 +232,7 @@ def calculate_rms_velocity( ages, cases, plateIDs, + PROGRESS_BAR, ) def calculate_net_rotation( @@ -238,6 +240,7 @@ def calculate_net_rotation( ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to calculate net rotation of the entire lithosphere. @@ -257,6 +260,7 @@ def calculate_net_rotation( ages, cases, plateIDs, + PROGRESS_BAR, ) def sample_all( @@ -264,6 +268,7 @@ def sample_all( ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to sample all variables relevant to the plate torques calculation. @@ -277,22 +282,23 @@ def sample_all( :type plateIDs: int, float, list, numpy.ndarray """ # Sample point seafloor ages - self.sample_point_seafloor_ages(ages, cases, plateIDs) + self.sample_point_seafloor_ages(ages, cases, plateIDs, PROGRESS_BAR) # Sample slab seafloor ages - self.sample_slab_seafloor_ages(ages, cases, plateIDs) + self.sample_slab_seafloor_ages(ages, cases, plateIDs, PROGRESS_BAR) # Sample arc seafloor ages - self.sample_arc_seafloor_ages(ages, cases, plateIDs) + self.sample_arc_seafloor_ages(ages, cases, plateIDs, PROGRESS_BAR) # Sample slab sediment thicknesses - self.sample_slab_sediment_thicknesses(ages, cases, plateIDs) + self.sample_slab_sediment_thicknesses(ages, cases, plateIDs, PROGRESS_BAR) def sample_seafloor_ages( self, ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to sample the seafloor ages and other variables (if available). @@ -316,6 +322,7 @@ def sample_point_seafloor_ages( ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to sample the seafloor ages and other variables (if available). @@ -327,13 +334,14 @@ def sample_point_seafloor_ages( :param plateIDs: plateIDs of interest (default: None) :type plateIDs: int, float, list, numpy.ndarray """ - self.points.sample_seafloor_ages(ages, cases, plateIDs, self.grids.seafloor_age) + self.points.sample_seafloor_ages(ages, cases, plateIDs, self.grids.seafloor_age, PROGRESS_BAR) def sample_slab_seafloor_ages( self, ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to sample the seafloor ages and other variables (if available). @@ -345,13 +353,14 @@ def sample_slab_seafloor_ages( :param plateIDs: plateIDs of interest (default: None) :type plateIDs: int, float, list, numpy.ndarray """ - self.slabs.sample_slab_seafloor_ages(ages, cases, plateIDs, self.grids.seafloor_age) + self.slabs.sample_slab_seafloor_ages(ages, cases, plateIDs, self.grids.seafloor_age, PROGRESS_BAR) def sample_arc_seafloor_ages( self, ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to sample the seafloor ages and other variables (if available). @@ -363,13 +372,16 @@ def sample_arc_seafloor_ages( :param plateIDs: plateIDs of interest (default: None) :type plateIDs: int, float, list, numpy.ndarray """ - self.slabs.sample_arc_seafloor_ages(ages, cases, plateIDs, self.grids.seafloor_age) + self.slabs.sample_arc_seafloor_ages(ages, cases, plateIDs, self.grids.seafloor_age, PROGRESS_BAR) + + self.slabs.set_continental_arc(ages, cases, plateIDs, PROGRESS_BAR) def sample_slab_sediment_thicknesses( self, ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to sample the seafloor ages and other variables (if available). @@ -382,16 +394,17 @@ def sample_slab_sediment_thicknesses( :type plateIDs: int, float, list, numpy.ndarray """ # Sample arc seafloor ages - self.sample_arc_seafloor_ages(ages, cases, plateIDs) + self.sample_arc_seafloor_ages(ages, cases, plateIDs, PROGRESS_BAR) # Sample slab sediment thickness - self.slabs.sample_slab_sediment_thickness(ages, cases, plateIDs, self.grids.sediment) + self.slabs.sample_slab_sediment_thickness(ages, cases, plateIDs, self.grids.sediment, PROGRESS_BAR) def calculate_all_torques( self, ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to calculate all torques. @@ -404,25 +417,25 @@ def calculate_all_torques( :type plateIDs: int, float, list, numpy.ndarray """ # Calculate slab pull torque - self.calculate_slab_pull_torque(ages, cases, plateIDs) + self.calculate_slab_pull_torque(ages, cases, plateIDs, PROGRESS_BAR) # Calculate GPE torque - self.calculate_gpe_torque(ages, cases, plateIDs) + self.calculate_gpe_torque(ages, cases, plateIDs, PROGRESS_BAR) # Calculate mantle drag torque - self.calculate_mantle_drag_torque(ages, cases, plateIDs) + self.calculate_mantle_drag_torque(ages, cases, plateIDs, PROGRESS_BAR) # Calculate slab bend torque - self.calculate_slab_bend_torque(ages, cases, plateIDs) + self.calculate_slab_bend_torque(ages, cases, plateIDs, PROGRESS_BAR) # Calculate synthetic velocity - self.calculate_synthetic_velocity(ages, cases, plateIDs) + self.calculate_synthetic_velocity(ages, cases, plateIDs, PROGRESS_BAR) # Calculate driving torque - self.calculate_driving_torque(ages, cases, plateIDs) + self.calculate_driving_torque(ages, cases, plateIDs, PROGRESS_BAR) # Calculate residual torque - self.calculate_residual_torque(ages, cases, plateIDs) + self.calculate_residual_torque(ages, cases, plateIDs, PROGRESS_BAR) logging.info("Calculated all torques!") @@ -431,6 +444,7 @@ def calculate_gpe_torque( ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to calculate the GPE torque. @@ -448,6 +462,7 @@ def calculate_gpe_torque( cases, plateIDs, self.grids.seafloor_age, + PROGRESS_BAR, ) # Calculate GPE torque acting on plate @@ -457,6 +472,7 @@ def calculate_gpe_torque( cases, plateIDs, torque_var = "GPE", + PROGRESS_BAR = PROGRESS_BAR, ) logging.info("Calculated GPE torque!") @@ -466,6 +482,7 @@ def calculate_slab_pull_torque( ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to calculate the slab pull torque. @@ -482,6 +499,7 @@ def calculate_slab_pull_torque( ages, cases, plateIDs, + PROGRESS_BAR, ) # Calculate slab pull force on plates @@ -491,6 +509,7 @@ def calculate_slab_pull_torque( cases, plateIDs, torque_var = "slab_pull", + PROGRESS_BAR = PROGRESS_BAR, ) logging.info("Calculated slab pull torque!") @@ -500,6 +519,7 @@ def calculate_mantle_drag_torque( ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to calculate the mantle drag torque. @@ -516,6 +536,7 @@ def calculate_mantle_drag_torque( ages, cases, plateIDs, + PROGRESS_BAR, ) # Calculate mantle drag force on plates @@ -525,6 +546,7 @@ def calculate_mantle_drag_torque( cases, plateIDs, torque_var = "mantle_drag", + PROGRESS_BAR = PROGRESS_BAR, ) logging.info("Calculated mantle drag torque!") @@ -534,6 +556,7 @@ def calculate_slab_bend_torque( ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to calculate the slab bend torque. @@ -546,7 +569,7 @@ def calculate_slab_bend_torque( :type plateIDs: int, float, list, numpy.ndarray """ # Calculate the slab bend force along the trenches - self.slabs.calculate_slab_bend_force(ages, cases, plateIDs) + self.slabs.calculate_slab_bend_force(ages, cases, plateIDs, PROGRESS_BAR) # Calculate the torque on self.plates.calculate_torque_on_plates( @@ -555,6 +578,7 @@ def calculate_slab_bend_torque( cases, plateIDs, torque_var = "slab_bend", + PROGRESS_BAR = PROGRESS_BAR, ) logging.info("Calculated slab bend torque!") @@ -564,6 +588,7 @@ def calculate_driving_torque( ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to calculate the driving torque. @@ -580,6 +605,7 @@ def calculate_driving_torque( ages, cases, plateIDs, + PROGRESS_BAR, ) def calculate_residual_torque( @@ -587,6 +613,7 @@ def calculate_residual_torque( ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to calculate the driving torque. @@ -603,6 +630,7 @@ def calculate_residual_torque( ages, cases, plateIDs, + PROGRESS_BAR, ) # Calculate residual torque at slabs @@ -611,6 +639,7 @@ def calculate_residual_torque( cases, plateIDs, type = "slabs", + PROGRESS_BAR = PROGRESS_BAR, ) # Calculate residual torque at points @@ -619,6 +648,7 @@ def calculate_residual_torque( cases, plateIDs, type = "points", + PROGRESS_BAR = PROGRESS_BAR, ) def calculate_residual_force( @@ -627,6 +657,7 @@ def calculate_residual_force( cases = None, plateIDs = None, type = "slabs", + PROGRESS_BAR = True, ): """ Function to calculate the residual forces. @@ -645,6 +676,7 @@ def calculate_residual_force( cases, plateIDs, self.plates.data, + PROGRESS_BAR, ) elif type == "points": @@ -654,6 +686,7 @@ def calculate_residual_force( cases, plateIDs, self.plates.data, + PROGRESS_BAR, ) else: @@ -664,6 +697,7 @@ def calculate_synthetic_velocity( ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to compute synthetic velocities. @@ -686,16 +720,23 @@ def calculate_synthetic_velocity( _ages, _cases, plateIDs, + PROGRESS_BAR, ) # Calculate net rotation - self.calculate_net_rotation(_ages, _cases) + self.calculate_net_rotation( + _ages, + _cases, + plateIDs, + PROGRESS_BAR, + ) # Calculate velocities at points self.points.calculate_velocities( _ages, _cases, self.plates.data, + PROGRESS_BAR, ) # Calculate velocities at slabs @@ -703,6 +744,7 @@ def calculate_synthetic_velocity( _ages, _cases, self.plates.data, + PROGRESS_BAR, ) # Calculate RMS velocity of plates @@ -711,6 +753,7 @@ def calculate_synthetic_velocity( _ages, _cases, plateIDs, + PROGRESS_BAR, ) def rotate_torque( @@ -721,6 +764,7 @@ def rotate_torque( ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to rotate a torque vector stored in another the Plates object to the reference frame of this Plates object. @@ -746,6 +790,7 @@ def rotate_torque( ages, cases, plateIDs, + PROGRESS_BAR, ) def extract_data_through_time( @@ -784,6 +829,7 @@ def save_all( cases = None, plateIDs = None, file_dir = None, + PROGRESS_BAR = True, ): """ Function to save all classes within the PlateTorques object. @@ -798,19 +844,19 @@ def save_all( :type file_dir: str """ # Save plates - self.save_plates(ages, cases, plateIDs, file_dir) + self.save_plates(ages, cases, plateIDs, file_dir, PROGRESS_BAR) # Save points - self.save_points(ages, cases, plateIDs, file_dir) + self.save_points(ages, cases, plateIDs, file_dir, PROGRESS_BAR) # Save slabs - self.save_slabs(ages, cases, plateIDs, file_dir) + self.save_slabs(ages, cases, plateIDs, file_dir, PROGRESS_BAR) # Save grids - self.save_grids(ages, cases, file_dir) + self.save_grids(ages, cases, file_dir, PROGRESS_BAR) # Save globe - self.save_globe(cases, file_dir) + self.save_globe(cases, file_dir, PROGRESS_BAR) def save_plates( self, @@ -818,6 +864,7 @@ def save_plates( cases = None, plateIDs = None, file_dir = None, + PROGRESS_BAR = True, ): """ Function to save plates. @@ -832,7 +879,7 @@ def save_plates( :type file_dir: str """ # Save plates - self.plates.save(ages, cases, plateIDs, file_dir) + self.plates.save(ages, cases, plateIDs, file_dir, PROGRESS_BAR) def save_points( self, @@ -840,6 +887,7 @@ def save_points( cases = None, plateIDs = None, file_dir = None, + PROGRESS_BAR = True, ): """ Function to save points. @@ -854,7 +902,7 @@ def save_points( :type file_dir: str """ # Save points - self.points.save(ages, cases, plateIDs, file_dir) + self.points.save(ages, cases, plateIDs, file_dir, PROGRESS_BAR) def save_slabs( self, @@ -862,6 +910,7 @@ def save_slabs( cases = None, plateIDs = None, file_dir = None, + PROGRESS_BAR = True, ): """ Function to save slabs. @@ -876,13 +925,14 @@ def save_slabs( :type file_dir: str """ # Save slabs - self.slabs.save(ages, cases, plateIDs, file_dir) + self.slabs.save(ages, cases, plateIDs, file_dir, PROGRESS_BAR) def save_grids( self, ages = None, cases = None, file_dir = None, + PROGRESS_BAR = True, ): """ Function to save grids. @@ -897,12 +947,13 @@ def save_grids( :type file_dir: str """ # Save grids - self.grids.save_all(ages, cases, file_dir) + self.grids.save_all(ages, cases, file_dir, PROGRESS_BAR) def save_globe( self, cases = None, file_dir = None, + PROGRESS_BAR = True, ): """ Function to save globe. @@ -917,7 +968,7 @@ def save_globe( :type file_dir: str """ # Save globe - self.globe.save(cases, file_dir) + self.globe.save(cases, file_dir, PROGRESS_BAR) def export_all( self, @@ -925,6 +976,7 @@ def export_all( cases = None, plateIDs = None, file_dir = None, + PROGRESS_BAR = True, ): """ Function to export all classes within the PlateTorques object. @@ -939,19 +991,19 @@ def export_all( :type file_dir: str """ # Save plates - self.export_plates(ages, cases, plateIDs, file_dir) + self.export_plates(ages, cases, plateIDs, file_dir, PROGRESS_BAR) # Save points - self.export_points(ages, cases, plateIDs, file_dir) + self.export_points(ages, cases, plateIDs, file_dir, PROGRESS_BAR) # Save slabs - self.export_slabs(ages, cases, plateIDs, file_dir) + self.export_slabs(ages, cases, plateIDs, file_dir, PROGRESS_BAR) # Save grids - self.save_grids(ages, cases, file_dir) + self.save_grids(ages, cases, file_dir, PROGRESS_BAR) # Save globe - self.export_globe(cases, file_dir) + self.export_globe(cases, file_dir, PROGRESS_BAR) def export_plates( self, @@ -959,6 +1011,7 @@ def export_plates( cases = None, plateIDs = None, file_dir = None, + PROGRESS_BAR = True, ): """ Function to export plates. @@ -973,7 +1026,7 @@ def export_plates( :type file_dir: str """ # Save plates - self.plates.export(ages, cases, plateIDs, file_dir) + self.plates.export(ages, cases, plateIDs, file_dir, PROGRESS_BAR) def export_points( self, @@ -981,6 +1034,7 @@ def export_points( cases = None, plateIDs = None, file_dir = None, + PROGRESS_BAR = True, ): """ Function to export points. @@ -995,7 +1049,7 @@ def export_points( :type file_dir: str """ # Save points - self.points.export(ages, cases, plateIDs, file_dir) + self.points.export(ages, cases, plateIDs, file_dir, PROGRESS_BAR) def export_slabs( self, @@ -1003,6 +1057,7 @@ def export_slabs( cases = None, plateIDs = None, file_dir = None, + PROGRESS_BAR = True, ): """ Function to export slabs. @@ -1017,13 +1072,14 @@ def export_slabs( :type file_dir: str """ # Save slabs - self.slabs.export(ages, cases, plateIDs, file_dir) + self.slabs.export(ages, cases, plateIDs, file_dir, PROGRESS_BAR) def export_grids( self, ages = None, cases = None, file_dir = None, + PROGRESS_BAR = True, ): """ Function to export grids. @@ -1036,12 +1092,13 @@ def export_grids( :type file_dir: str """ # Save grids - self.grids.export(ages, cases, file_dir) + self.grids.export(ages, cases, file_dir, PROGRESS_BAR) def export_globe( self, cases = None, file_dir = None, + PROGRESS_BAR = True, ): """ Function to expprt globe. @@ -1052,4 +1109,4 @@ def export_globe( :type file_dir: str """ # Save globe - self.globe.export(cases, file_dir) \ No newline at end of file + self.globe.export(cases, file_dir, PROGRESS_BAR) \ No newline at end of file diff --git a/plato/plates.py b/plato/plates.py index 36014e6..564b4b1 100644 --- a/plato/plates.py +++ b/plato/plates.py @@ -130,7 +130,7 @@ def __init__( # Copy all DataFrames from the available case for entries in entry: if entry not in available_cases: - self.geometries[_age][entry] = self.geometries[_age][available_cases[0]].copy() + self.resolved_geometries[_age][entry] = self.resolved_geometries[_age][available_cases[0]].copy() else: # Initialise missing geometries self.resolved_geometries[_age][key] = utils_data.get_resolved_geometries( @@ -184,7 +184,7 @@ def __init__( # Check if any DataFrames were loaded if len(available_cases) > 0: # Copy all DataFrames from the available case - for entries in entry: + for entry in entries: if entry not in available_cases: self.data[_age][entry] = self.data[_age][available_cases[0]].copy() else: @@ -445,11 +445,11 @@ def calculate_driving_torque( logging.info("Computing driving torques...") # Loop through ages - for i, _age in _tqdm( - enumerate(_ages), + for i, _age in enumerate(_tqdm( + _ages, desc="Calculating driving torque", disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) - ): + )): # Loop through cases for _case in _cases: # Select plates diff --git a/plato/points.py b/plato/points.py index 6b90248..e6270d6 100644 --- a/plato/points.py +++ b/plato/points.py @@ -1,15 +1,18 @@ # Standard libraries import logging +from typing import Dict, List, Optional, Union # Third-party libraries import geopandas as _geopandas +import gplately as _gplately +import numpy as _numpy import pandas as _pandas import xarray as _xarray - from tqdm import tqdm as _tqdm # Local libraries from . import utils_data, utils_calc, utils_init +from .settings import Settings class Points: """ @@ -53,20 +56,21 @@ class Points: """ def __init__( self, - settings = None, - reconstruction = None, - rotation_file = None, - topology_file = None, - polygon_file = None, - reconstruction_name = None, - ages = None, - cases_file = None, - cases_sheet = "Sheet1", - files_dir = None, - resolved_geometries = None, - PARALLEL_MODE = False, - DEBUG_MODE = False, - CALCULATE_VELOCITIES = True, + settings: Optional[Settings] = None, + reconstruction: Optional[_gplately.PlateReconstruction] = None, + rotation_file: Optional[str] = None, + topology_file: Optional[str] = None, + polygon_file: Optional[str] = None, + reconstruction_name: Optional[str] = None, + ages: Optional[Union[int, float, List[Union[int, float]], _numpy.ndarray]] = None, + cases_file: Optional[str] = None, + cases_sheet: str = "Sheet1", + files_dir: Optional[str] = None, + resolved_geometries: Dict[float, Dict[str, _geopandas.GeoDataFrame]] = None, + PARALLEL_MODE: bool = False, + DEBUG_MODE: bool = False, + CALCULATE_VELOCITIES: bool = True, + PROGRESS_BAR: bool = True, ): """ Constructor for the 'Points' object. @@ -118,9 +122,10 @@ def __init__( # Check if any DataFrames were loaded if len(available_cases) > 0: # Copy all DataFrames from the available case - for entries in entry: + for entry in entries: if entry not in available_cases: self.data[_age][entry] = self.data[_age][available_cases[0]].copy() + else: # Initialise missing data if not isinstance(resolved_geometries, dict) or not isinstance(resolved_geometries.get(key), _geopandas.GeoDataFrame): @@ -163,6 +168,7 @@ def calculate_velocities( ages = None, cases = None, stage_rotation = None, + PROGRESS_BAR = True, ): """ Function to compute velocities at points. @@ -181,7 +187,11 @@ def calculate_velocities( _cases = utils_data.select_cases(cases, self.settings.point_cases) # Loop through ages and cases - for _age in _ages: + for _age in _tqdm( + _ages, + desc="Calculating velocities at points", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): for _case in _cases: for plateID in self.data[_age][_case].plateID.unique(): # Get stage rotation, if not provided @@ -234,6 +244,8 @@ def sample_seafloor_ages( cases = None, plateIDs = None, seafloor_grids = None, + vars = ["seafloor_age"], + PROGRESS_BAR = True, ): """ Samples seafloor age at points. @@ -254,6 +266,7 @@ def sample_seafloor_ages( plateIDs, seafloor_grids, vars, + PROGRESS_BAR = PROGRESS_BAR, ) # Set sampling flag to true @@ -268,6 +281,7 @@ def sample_grid( vars = ["seafloor_age"], sampling_coords = ["lat", "lon"], cols = ["seafloor_age"], + PROGRESS_BAR = True, ): """ Samples any grid at points. @@ -302,7 +316,11 @@ def sample_grid( _vars = [] # Loop through valid times - for _age in _tqdm(_ages, desc="Sampling points", disable=self.settings.logger.level == logging.INFO): + for _age in _tqdm( + _ages, + desc="Sampling points", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): for key, entries in _iterable.items(): # Define plateIDs if not provided _plateIDs = utils_data.select_plateIDs(plateIDs, self.data[_age][key].plateID.unique()) @@ -355,6 +373,7 @@ def calculate_gpe_force( cases = None, plateIDs = None, seafloor_grid = None, + PROGRESS_BAR = True, ): """ Function to compute gravitational potential energy (GPE) force acting at points. @@ -375,7 +394,11 @@ def calculate_gpe_force( _iterable = utils_data.select_iterable(cases, self.settings.gpe_cases) # Loop through reconstruction times - for _age in _tqdm(_ages, desc="Computing GPE forces", disable=(self.settings.logger.level==logging.INFO)): + for _age in _tqdm( + _ages, + desc="Calculating GPE forces", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): # Loop through gpe cases for key, entries in _iterable.items(): if self.settings.options[key]["GPE torque"]: @@ -423,6 +446,7 @@ def calculate_mantle_drag_force( ages = None, cases = None, plateIDs = None, + PROGRESS_BAR = True, ): """ Function to compute mantle drag force acting at points. @@ -441,7 +465,11 @@ def calculate_mantle_drag_force( _iterable = utils_data.select_iterable(cases, self.settings.gpe_cases) # Loop through reconstruction times - for _age in _tqdm(_ages, desc="Computing mantle drag forces", disable=(self.settings.logger.level==logging.INFO)): + for _age in _tqdm( + _ages, + desc="Calculating mantle drag forces", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): # Loop through gpe cases for key, entries in _iterable.items(): if self.settings.options[key]["Mantle drag torque"] and self.settings.options[key]["Reconstructed motions"]: @@ -483,6 +511,7 @@ def calculate_residual_force( cases = None, plateIDs = None, residual_torque = None, + PROGRESS_BAR = True, ): """ Function to calculate residual torque along trenches. @@ -503,37 +532,42 @@ def calculate_residual_force( _cases = utils_data.select_iterable(cases, self.settings.slab_pull_cases) # Loop through ages and cases - for _age in _ages: - for _case in _cases: - # Select plateIDs - _plateIDs = utils_data.select_plateIDs(plateIDs, self.data[_age][_case]["plateID"].unique()) - - for _plateID in _plateIDs: - if ( - isinstance(residual_torque, dict) - and _age in residual_torque.keys() - and _case in residual_torque[_age].keys() - and isinstance(residual_torque[_age][_case], _pandas.DataFrame) - ): - # Get stage rotation from the provided DataFrame in the dictionary - _residual_torque = residual_torque[_age][_case][residual_torque[_age][_case].plateID == _plateID] - - # Make mask for plate - mask = self.data[_age][_case]["plateID"] == _plateID - - # Compute velocities - forces = utils_calc.compute_residual_force( - self.data[_age][_case].loc[mask], - _residual_torque, - plateID_col = "plateID", - weight_col = "segment_area", - ) + for _case in _tqdm( + _cases, + desc="Calculating residual forces at points", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): + if self.settings.options[_case]["Reconstructed motions"]: + for _age in _ages: + # Select plateIDs + _plateIDs = utils_data.select_plateIDs(plateIDs, self.data[_age][_case]["plateID"].unique()) + + for _plateID in _plateIDs: + if ( + isinstance(residual_torque, dict) + and _age in residual_torque.keys() + and _case in residual_torque[_age].keys() + and isinstance(residual_torque[_age][_case], _pandas.DataFrame) + ): + # Get stage rotation from the provided DataFrame in the dictionary + _residual_torque = residual_torque[_age][_case][residual_torque[_age][_case].plateID == _plateID] + + # Make mask for plate + mask = self.data[_age][_case]["plateID"] == _plateID + + # Compute velocities + forces = utils_calc.compute_residual_force( + self.data[_age][_case].loc[mask], + _residual_torque, + plateID_col = "plateID", + weight_col = "segment_area", + ) - # Store velocities - self.data[_age][_case].loc[mask, "residual_force_lat"] = forces[0] - self.data[_age][_case].loc[mask, "residual_force_lon"] = forces[1] - self.data[_age][_case].loc[mask, "residual_force_mag"] = forces[2] - self.data[_age][_case].loc[mask, "residual_force_azi"] = forces[3] + # Store velocities + self.data[_age][_case].loc[mask, "residual_force_lat"] = forces[0] + self.data[_age][_case].loc[mask, "residual_force_lon"] = forces[1] + self.data[_age][_case].loc[mask, "residual_force_mag"] = forces[2] + self.data[_age][_case].loc[mask, "residual_force_azi"] = forces[3] def save( self, @@ -541,6 +575,7 @@ def save( cases = None, plateIDs = None, file_dir = None, + PROGRESS_BAR = True, ): """ Function to save the 'Points' object. @@ -565,7 +600,11 @@ def save( _file_dir = self.settings.dir_path if file_dir is None else file_dir # Loop through ages - for _age in _tqdm(_ages, desc="Saving Points", disable=self.settings.logger.level==logging.INFO): + for _age in _tqdm( + _ages, + desc="Saving Points", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): # Loop through cases for _case in _cases: # Define plateIDs if not provided @@ -592,6 +631,7 @@ def export( cases = None, plateIDs = None, file_dir = None, + PROGRESS_BAR = True, ): """ Function to export the 'Points' object. @@ -616,7 +656,11 @@ def export( _file_dir = self.settings.dir_path if file_dir is None else file_dir # Loop through ages - for _age in _tqdm(_ages, desc="Exporting Points", disable=self.settings.logger.level==logging.INFO): + for _age in _tqdm( + _ages, + desc="Exporting Points", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): # Loop through cases for _case in _cases: # Define plateIDs if not provided diff --git a/plato/settings.py b/plato/settings.py index 2ee0626..5b0b4ab 100644 --- a/plato/settings.py +++ b/plato/settings.py @@ -127,7 +127,7 @@ def __init__( # NOTE: this is hardcoded for lack of a better alternative # ^* Earthbyte reconstructions: Seton et al. (2012), Müller et al. (2016, 2019), Matthews et al. (2016), Clennett et al. (2020) self.oceanic_arc_plateIDs = [ - 529, # Kohistan-Ladakh + 518, 529, # Kohistan-Ladakh 608, 659, 699, # Izu-Bonin-Marianas 612, # Luzon 645, # East Sunda @@ -140,12 +140,13 @@ def __init__( 827, # New Hebrides 847, # Vityaz 853, # West Solomon Sea - 844, 841, 865, 943, # Junction + 609, 844, 841, 865, 869, 943, # Junction 1072, 1073, 1080, # Insular 2007, # Antilles 9052, 95104, # Central America 9022, # Cascadia root 9040, # Angayucham + 67350, # Woyla ] logging.info("Settings initialisation complete.") diff --git a/plato/slabs.py b/plato/slabs.py index 2a59dc5..44962bc 100644 --- a/plato/slabs.py +++ b/plato/slabs.py @@ -1,6 +1,8 @@ +# Standard libraries import logging from typing import Dict, List, Optional, Union +# Third party libraries import geopandas as _geopandas import gplately as _gplately import numpy as _numpy @@ -8,6 +10,7 @@ import xarray as _xarray from tqdm import tqdm as _tqdm +# Local libraries from . import utils_data, utils_calc, utils_init from .settings import Settings @@ -128,7 +131,7 @@ def __init__( # Check if any DataFrames were loaded if len(available_cases) > 0: # Copy all DataFrames from the available case - for entries in entry: + for entry in entries: if entry not in available_cases: self.data[_age][entry] = self.data[_age][available_cases[0]].copy() else: @@ -198,7 +201,11 @@ def calculate_velocities( _cases = utils_data.select_cases(cases, self.settings.cases) # Loop through ages and cases - for _age in _ages: + for _age in _tqdm( + _ages, + desc="Calculating velocities at slabs", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): for _case in _cases: for plate in ["upper_plate", "lower_plate", "trench"]: plateID_col = f"{plate}ID" if plate != "trench" else "trench_plateID" @@ -256,6 +263,7 @@ def sample_slab_seafloor_ages( cases: Optional[Union[str, List[str]]] = None, plateIDs: Optional[Union[int, float, List[Union[int, float]], _numpy.ndarray]] = None, grids = None, + PROGRESS_BAR: bool = True, ): """ Samples seafloor age at slabs. @@ -280,6 +288,7 @@ def sample_slab_seafloor_ages( plate = "lower", vars = ["seafloor_age"], cols = ["slab_seafloor_age"], + PROGRESS_BAR = PROGRESS_BAR, ) # Set sampling flag to true @@ -291,6 +300,7 @@ def sample_arc_seafloor_ages( cases = None, plateIDs = None, grids = None, + PROGRESS_BAR = True, ): """ Samples seafloor age at arcs. @@ -315,6 +325,7 @@ def sample_arc_seafloor_ages( plate = "upper", vars = ["seafloor_age"], cols = ["arc_seafloor_age"], + PROGRESS_BAR = PROGRESS_BAR, ) # Set sampling flag to true @@ -326,6 +337,7 @@ def sample_slab_sediment_thickness( cases: Optional[Union[str, List[str]]] = None, plateIDs: Optional[Union[int, float, List[Union[int, float]], _numpy.ndarray]] = None, grids = None, + PROGRESS_BAR: bool = True, ): """ Samples sediment thickness at slabs. @@ -357,6 +369,7 @@ def sample_slab_sediment_thickness( plate = "lower", vars = None, cols = ["sediment_thickness"], + PROGRESS_BAR = PROGRESS_BAR, ) # Get continental arcs @@ -364,6 +377,7 @@ def sample_slab_sediment_thickness( ages, cases, plateIDs, + PROGRESS_BAR, ) # Set active margin sediment thickness @@ -372,6 +386,7 @@ def sample_slab_sediment_thickness( cases, plateIDs, cols = "sediment_thickness", + PROGRESS_BAR = PROGRESS_BAR, ) # Set sampling flag to true @@ -382,6 +397,7 @@ def set_continental_arc( ages: Optional[Union[int, float, List[Union[int, float]], _numpy.ndarray]] = None, cases: Optional[Union[str, List[str]]] = None, plateIDs: Optional[Union[int, float, List[Union[int, float]], _numpy.ndarray]] = None, + PROGRESS_BAR: bool = True, ): """ Function to set whether a trench has a continental arc. @@ -397,7 +413,7 @@ def set_continental_arc( """ # Check whether arc seafloor ages have been sampled if not self.sampled_seafloor_at_arcs: - self.sample_arc_seafloor_ages(ages, cases, plateIDs) + self.sample_arc_seafloor_ages(ages, cases, plateIDs, PROGRESS_BAR) # Define ages if not provided _ages = utils_data.select_ages(ages, self.settings.ages) @@ -441,6 +457,7 @@ def sample_grid( plate: Optional[str] = "lower", vars: Optional[Union[str, List[str]]] = ["seafloor_age"], cols = ["slab_seafloor_age"], + PROGRESS_BAR: bool = True, ): """ Samples any grid at slabs. @@ -456,7 +473,11 @@ def sample_grid( # Loop through valid cases # Order of loops is flipped to skip cases where no grid needs to be sampled - for key, entries in _tqdm(_iterable.items(), desc="Sampling slabs", disable=self.settings.logger.level == logging.INFO): + for key, entries in _tqdm( + _iterable.items(), + desc="Sampling slabs", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): # Skip if sediment grid is not sampled if cols == ["sediment_thickness"] and not self.settings.options[key]["Sample sediment grid"]: logging.info(f"Skipping sampling of sediment thickness for case {key}.") @@ -559,6 +580,7 @@ def set_values( plate: Optional[str] = "lower", cols: Optional[Union[List[str], str]] = None, vals: Optional[Union[List[float], float]] = None, + PROGRESS_BAR: bool = True, ): """ Function to add values to the 'Slabs' object. @@ -573,7 +595,11 @@ def set_values( type = "arc" if plate == "upper" else "slab" # Loop through valid cases - for key, entries in _tqdm(_iterable.items(), desc="Adding values", disable=(self.settings.logger.level==logging.INFO)): + for key, entries in _tqdm( + _iterable.items(), + desc="Adding values", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): # Special case for active margin sediments if cols == "sediment_thickness" and vals == None: if not self.settings.options[key]["Active margin sediments"]: @@ -630,6 +656,7 @@ def calculate_slab_pull_force( ages: Optional[Union[int, float, List[Union[int, float]], _numpy.ndarray]] = None, cases: Optional[Union[str, List[str]]] = None, plateIDs: Optional[Union[int, float, List[Union[int, float]], _numpy.ndarray]] = None, + PROGRESS_BAR: bool = True, ): """ Function to compute slab pull force along trenches. @@ -642,7 +669,11 @@ def calculate_slab_pull_force( # Loop through valid cases # Order of loops is flipped to skip cases where no slab pull torque needs to be sampled - for key, entries in _tqdm(_iterable.items(), desc="Computing slab pull forces", disable=(self.settings.logger.level==logging.INFO)): + for key, entries in _tqdm( + _iterable.items(), + desc="Calculating slab pull forces", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): # Skip if slab pull torque is not sampled if self.settings.options[key]["Slab pull torque"]: # Loop through ages @@ -700,6 +731,7 @@ def calculate_slab_bend_force( ages: Optional[Union[int, float, List[Union[int, float]], _numpy.ndarray]] = None, cases: Optional[Union[str, List[str]]] = None, plateIDs: Optional[Union[int, float, List[Union[int, float]], _numpy.ndarray]] = None, + PROGRESS_BAR: bool = True, ): """ Function to compute slab bend force along trenches. @@ -712,7 +744,11 @@ def calculate_slab_bend_force( # Loop through valid cases # Order of loops is flipped to skip cases where no slab bend torque needs to be sampled - for key, entries in _tqdm(_iterable.items(), desc="Computing slab bend forces", disable=(self.settings.logger.level==logging.INFO)): + for key, entries in _tqdm( + _iterable.items(), + desc="Calculating slab bend forces", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): # Skip if slab pull torque is not sampled if self.settings.options[key]["Slab bend torque"]: # Loop through ages @@ -762,6 +798,7 @@ def calculate_residual_force( cases: Optional[Union[str, List[str]]] = None, plateIDs: Optional[Union[int, float, List[Union[int, float]], _numpy.ndarray]] = None, residual_torque: Optional[Dict] = None, + PROGRESS_BAR: bool = True, ): """ Function to calculate residual torque along trenches. @@ -773,37 +810,42 @@ def calculate_residual_force( _cases = utils_data.select_iterable(cases, self.settings.slab_pull_cases) # Loop through ages and cases - for _age in _ages: - for _case in _cases: - # Select plateIDs - _plateIDs = utils_data.select_plateIDs(plateIDs, self.data[_age][_case]["lower_plateID"].unique()) - - for _plateID in _plateIDs: - if ( - isinstance(residual_torque, Dict) - and _age in residual_torque.keys() - and _case in residual_torque[_age].keys() - and isinstance(residual_torque[_age][_case], _pandas.DataFrame) - ): - # Get stage rotation from the provided DataFrame in the dictionary - _residual_torque = residual_torque[_age][_case][residual_torque[_age][_case].plateID == _plateID] - - # Make mask for plate - mask = self.data[_age][_case]["lower_plateID"] == _plateID - - # Compute velocities - forces = utils_calc.compute_residual_force( - self.data[_age][_case].loc[mask], - _residual_torque, - plateID_col = "lower_plateID", - weight_col = "trench_normal_azimuth", - ) + for _case in _tqdm( + _cases, + desc="Calculating residual forces at slabs", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): + if self.settings.options[_case]["Reconstructed motions"]: + for _age in _ages: + # Select plateIDs + _plateIDs = utils_data.select_plateIDs(plateIDs, self.data[_age][_case]["lower_plateID"].unique()) - # Store velocities - self.data[_age][_case].loc[mask, "residual_force_lat"] = forces[0] - self.data[_age][_case].loc[mask, "residual_force_lon"] = forces[1] - self.data[_age][_case].loc[mask, "residual_force_mag"] = forces[2] - self.data[_age][_case].loc[mask, "residual_force_azi"] = forces[3] + for _plateID in _plateIDs: + if ( + isinstance(residual_torque, Dict) + and _age in residual_torque.keys() + and _case in residual_torque[_age].keys() + and isinstance(residual_torque[_age][_case], _pandas.DataFrame) + ): + # Get stage rotation from the provided DataFrame in the dictionary + _residual_torque = residual_torque[_age][_case][residual_torque[_age][_case].plateID == _plateID] + + # Make mask for plate + mask = self.data[_age][_case]["lower_plateID"] == _plateID + + # Compute velocities + forces = utils_calc.compute_residual_force( + self.data[_age][_case].loc[mask], + _residual_torque, + plateID_col = "lower_plateID", + weight_col = "trench_normal_azimuth", + ) + + # Store velocities + self.data[_age][_case].loc[mask, "residual_force_lat"] = forces[0] + self.data[_age][_case].loc[mask, "residual_force_lon"] = forces[1] + self.data[_age][_case].loc[mask, "residual_force_mag"] = forces[2] + self.data[_age][_case].loc[mask, "residual_force_azi"] = forces[3] def extract_data_through_time( self, @@ -884,6 +926,7 @@ def save( cases: Optional[Union[str, List[str]]] = None, plateIDs: Optional[Union[int, float, List[Union[int, float]], _numpy.ndarray]] = None, file_dir: Optional[str] = None, + PROGRESS_BAR: bool = True, ): """ Function to save the 'Slabs' object. @@ -901,7 +944,11 @@ def save( _file_dir = self.settings.dir_path if file_dir is None else file_dir # Loop through ages - for _age in _tqdm(_ages, desc="Saving Slabs", disable=self.settings.logger.level==logging.INFO): + for _age in _tqdm( + _ages, + desc="Saving Slabs", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): # Loop through cases for _case in _cases: utils_data.DataFrame_to_parquet( @@ -921,6 +968,7 @@ def export( cases: Optional[Union[str, List[str]]] = None, plateIDs: Optional[Union[int, float, List[Union[int, float]], _numpy.ndarray]] = None, file_dir: Optional[str] = None, + PROGRESS_BAR: bool = True, ): """ Function to export the 'Slabs' object. @@ -938,7 +986,11 @@ def export( _file_dir = self.settings.dir_path if file_dir is None else file_dir # Loop through ages - for _age in _tqdm(_ages, desc="Exporting Slabs", disable=self.settings.logger.level==logging.INFO): + for _age in _tqdm( + _ages, + desc="Exporting Slabs", + disable=(self.settings.logger.level in [logging.INFO, logging.DEBUG] or not PROGRESS_BAR) + ): # Loop through cases for _case in _cases: utils_data.DataFrame_to_csv( diff --git a/plato/utils_calc.py b/plato/utils_calc.py index af37c7b..ab7429c 100644 --- a/plato/utils_calc.py +++ b/plato/utils_calc.py @@ -561,7 +561,7 @@ def compute_synthetic_stage_rotation( stage_rotation_poles_mag = _numpy.nan_to_num(stage_rotation_poles_mag) # Normalise the rotation poles by the drag coefficient and the square of the Earth's radius - stage_rotation_poles_mag /= options["Mantle viscosity"] / mech.La * constants.mean_Earth_radius_m**2 + stage_rotation_poles_mag /= options["Mantle viscosity"] / mech.La * plates.area #constants.mean_Earth_radius_m**2 # Convert to degrees because the 'geocentric_cartesian2spherical' does not convert the magnitude to degrees stage_rotation_poles_mag = _numpy.rad2deg(stage_rotation_poles_mag) diff --git a/plato/utils_data.py b/plato/utils_data.py index b4ac634..e3a2ea9 100644 --- a/plato/utils_data.py +++ b/plato/utils_data.py @@ -594,7 +594,7 @@ def get_options( 1.23e20, 250, 1, - 7.5e12, + 0, 0, 1, 0.1,