diff --git a/LoopStructural/modelling/core/geological_model.py b/LoopStructural/modelling/core/geological_model.py index e8a10a8f..ce501df3 100644 --- a/LoopStructural/modelling/core/geological_model.py +++ b/LoopStructural/modelling/core/geological_model.py @@ -125,9 +125,7 @@ def __init__( self.features = [] self.feature_name_index = {} self._data = pd.DataFrame() # None - - self.stratigraphic_column = None self.tol = 1e-10 * np.max(self.bounding_box.maximum - self.bounding_box.origin) @@ -179,8 +177,40 @@ def __str__(self): def _ipython_key_completions_(self): return self.feature_name_index.keys() - + def prepare_data(self, data: pd.DataFrame) -> pd.DataFrame: + data = data.copy() + data[['X', 'Y', 'Z']] = self.bounding_box.project(data[['X', 'Y', 'Z']].to_numpy()) + if "type" in data: + logger.warning("'type' is deprecated replace with 'feature_name' \n") + data.rename(columns={"type": "feature_name"}, inplace=True) + if "feature_name" not in data: + logger.error("Data does not contain 'feature_name' column") + raise BaseException("Cannot load data") + for h in all_heading(): + if h not in data: + data[h] = np.nan + if h == "w": + data[h] = 1.0 + if h == "coord": + data[h] = 0 + if h == "polarity": + data[h] = 1.0 + # LS wants polarity as -1 or 1, change 0 to -1 + data.loc[data["polarity"] == 0, "polarity"] = -1.0 + data.loc[np.isnan(data["w"]), "w"] = 1.0 + if "strike" in data and "dip" in data: + logger.info("Converting strike and dip to vectors") + mask = np.all(~np.isnan(data.loc[:, ["strike", "dip"]]), axis=1) + data.loc[mask, gradient_vec_names()] = ( + strikedip2vector(data.loc[mask, "strike"], data.loc[mask, "dip"]) + * data.loc[mask, "polarity"].to_numpy()[:, None] + ) + data.drop(["strike", "dip"], axis=1, inplace=True) + data[['X', 'Y', 'Z', 'val', 'nx', 'ny', 'nz', 'gx', 'gy', 'gz', 'tx', 'ty', 'tz']] = data[ + ['X', 'Y', 'Z', 'val', 'nx', 'ny', 'nz', 'gx', 'gy', 'gz', 'tx', 'ty', 'tz'] + ].astype(float) + return data @classmethod def from_processor(cls, processor): """Builds a model from a :class:`LoopStructural.modelling.input.ProcessInputData` object @@ -372,12 +402,7 @@ def fault_names(self): """ return [f.name for f in self.faults] - def check_inialisation(self): - if self.data is None: - logger.error("Data not associated with GeologicalModel. Run set_data") - return False - if self.data.shape[0] > 0: - return True + def to_file(self, file): """Save a model to a pickle file requires dill @@ -473,40 +498,9 @@ def data(self, data: pd.DataFrame): raise BaseException("Cannot load data") logger.info(f"Adding data to GeologicalModel with {len(data)} data points") self._data = data.copy() - self._data[['X','Y','Z']] = self.bounding_box.project(self._data[['X','Y','Z']].to_numpy()) - - - if "type" in self._data: - logger.warning("'type' is deprecated replace with 'feature_name' \n") - self._data.rename(columns={"type": "feature_name"}, inplace=True) - if "feature_name" not in self._data: - logger.error("Data does not contain 'feature_name' column") - raise BaseException("Cannot load data") - for h in all_heading(): - if h not in self._data: - self._data[h] = np.nan - if h == "w": - self._data[h] = 1.0 - if h == "coord": - self._data[h] = 0 - if h == "polarity": - self._data[h] = 1.0 - # LS wants polarity as -1 or 1, change 0 to -1 - self._data.loc[self._data["polarity"] == 0, "polarity"] = -1.0 - self._data.loc[np.isnan(self._data["w"]), "w"] = 1.0 - if "strike" in self._data and "dip" in self._data: - logger.info("Converting strike and dip to vectors") - mask = np.all(~np.isnan(self._data.loc[:, ["strike", "dip"]]), axis=1) - self._data.loc[mask, gradient_vec_names()] = ( - strikedip2vector(self._data.loc[mask, "strike"], self._data.loc[mask, "dip"]) - * self._data.loc[mask, "polarity"].to_numpy()[:, None] - ) - self._data.drop(["strike", "dip"], axis=1, inplace=True) - self._data[['X', 'Y', 'Z', 'val', 'nx', 'ny', 'nz', 'gx', 'gy', 'gz', 'tx', 'ty', 'tz']] = ( - self._data[ - ['X', 'Y', 'Z', 'val', 'nx', 'ny', 'nz', 'gx', 'gy', 'gz', 'tx', 'ty', 'tz'] - ].astype(float) - ) + # self._data[['X','Y','Z']] = self.bounding_box.project(self._data[['X','Y','Z']].to_numpy()) + + def set_model_data(self, data): logger.warning("deprecated method. Model data can now be set using the data attribute") @@ -601,9 +595,7 @@ def create_and_add_foliation( An interpolator will be chosen by calling :meth:`LoopStructural.GeologicalModel.get_interpolator` """ - if not self.check_inialisation(): - logger.warning(f"{series_surface_data} not added, model not initialised") - return + # if tol is not specified use the model default if tol is None: tol = self.tol @@ -623,7 +615,7 @@ def create_and_add_foliation( if series_surface_data.shape[0] == 0: logger.warning("No data for {series_surface_data}, skipping") return - series_builder.add_data_from_data_frame(series_surface_data) + series_builder.add_data_from_data_frame(self.prepare_data(series_surface_data)) self._add_faults(series_builder, features=faults) # build feature @@ -675,8 +667,7 @@ def create_and_add_fold_frame( fold_frame : FoldFrame the created fold frame """ - if not self.check_inialisation(): - return False + if tol is None: tol = self.tol @@ -697,7 +688,7 @@ def create_and_add_fold_frame( if fold_frame_data.shape[0] == 0: logger.warning(f"No data for {fold_frame_name}, skipping") return - fold_frame_builder.add_data_from_data_frame(fold_frame_data) + fold_frame_builder.add_data_from_data_frame(self.prepare_data(fold_frame_data)) self._add_faults(fold_frame_builder[0]) self._add_faults(fold_frame_builder[1]) self._add_faults(fold_frame_builder[2]) @@ -752,8 +743,7 @@ def create_and_add_folded_foliation( :class:`LoopStructural.modelling.features.builders.FoldedFeatureBuilder` """ - if not self.check_inialisation(): - return False + if tol is None: tol = self.tol @@ -783,7 +773,10 @@ def create_and_add_folded_foliation( logger.warning(f"No data for {foliation_name}, skipping") return series_builder.add_data_from_data_frame( -foliation_data ) + self.prepare_data( + foliation_data + ) + ) self._add_faults(series_builder) # series_builder.add_data_to_interpolator(True) # build feature @@ -850,8 +843,7 @@ def create_and_add_folded_fold_frame( see :class:`LoopStructural.modelling.features.fold.FoldEvent`, :class:`LoopStructural.modelling.features.builders.FoldedFeatureBuilder` """ - if not self.check_inialisation(): - return False + if tol is None: tol = self.tol @@ -878,7 +870,7 @@ def create_and_add_folded_fold_frame( ) if fold_frame_data is None: fold_frame_data = self.data[self.data["feature_name"] == fold_frame_name] - fold_frame_builder.add_data_from_data_frame(fold_frame_data) + fold_frame_builder.add_data_from_data_frame(self.prepare_data(fold_frame_data)) for i in range(3): self._add_faults(fold_frame_builder[i]) @@ -1331,7 +1323,7 @@ def create_and_add_fault( if fault_data.shape[0] == 0: logger.warning(f"No data for {fault_name}, skipping") return - + self._add_faults(fault_frame_builder, features=faults) # add data @@ -1344,7 +1336,7 @@ def create_and_add_fault( if intermediate_axis: intermediate_axis = intermediate_axis fault_frame_builder.create_data_from_geometry( - fault_frame_data=fault_data, + fault_frame_data=self.prepare_data(fault_data), fault_center=fault_center, fault_normal_vector=fault_normal_vector, fault_slip_vector=fault_slip_vector, @@ -1397,9 +1389,8 @@ def rescale(self, points: np.ndarray, *, inplace: bool = False) -> np.ndarray: points : np.array((N,3),dtype=double) """ - + return self.bounding_box.reproject(points,inplace=inplace) - # TODO move scale to bounding box/transformer def scale(self, points: np.ndarray, *, inplace: bool = False) -> np.ndarray: @@ -1418,7 +1409,6 @@ def scale(self, points: np.ndarray, *, inplace: bool = False) -> np.ndarray: """ return self.bounding_box.project(np.array(points).astype(float),inplace=inplace) - def regular_grid(self, *, nsteps=None, shuffle=True, rescale=False, order="C"): """ diff --git a/LoopStructural/modelling/features/builders/_fault_builder.py b/LoopStructural/modelling/features/builders/_fault_builder.py index b9974586..aa87dc58 100644 --- a/LoopStructural/modelling/features/builders/_fault_builder.py +++ b/LoopStructural/modelling/features/builders/_fault_builder.py @@ -150,6 +150,7 @@ def create_data_from_geometry( np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0), ["X", "Y"], ].to_numpy() + self.fault_dip = fault_dip if fault_normal_vector is None: if fault_frame_data.loc[ np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["nx"].notna())].shape[0]>0: diff --git a/LoopStructural/modelling/features/fault/_fault_segment.py b/LoopStructural/modelling/features/fault/_fault_segment.py index e98f945a..eec98598 100644 --- a/LoopStructural/modelling/features/fault/_fault_segment.py +++ b/LoopStructural/modelling/features/fault/_fault_segment.py @@ -104,7 +104,7 @@ def fault_major_axis(self): if self.builder is None: raise ValueError("Fault builder not set") return self.builder.fault_major_axis - + @property def fault_intermediate_axis(self): if self.builder is None: diff --git a/tests/unit/modelling/intrusions/test_intrusions.py b/tests/unit/modelling/intrusions/test_intrusions.py index 377ccc2e..0f075125 100644 --- a/tests/unit/modelling/intrusions/test_intrusions.py +++ b/tests/unit/modelling/intrusions/test_intrusions.py @@ -65,9 +65,9 @@ def test_intrusion_builder(): model.data = data model.nsteps = [10, 10, 10] - intrusion_data = data[data["feature_name"] == "tabular_intrusion"] - intrusion_frame_data = model.data[model.data["feature_name"] == "tabular_intrusion_frame"] - + intrusion_data = model.prepare_data(data[data["feature_name"] == "tabular_intrusion"]) + intrusion_frame_data = model.prepare_data(model.data[model.data["feature_name"] == "tabular_intrusion_frame"]) + conformable_feature = model.create_and_add_foliation("stratigraphy") intrusion_frame_parameters = { diff --git a/tests/unit/modelling/test_geological_model.py b/tests/unit/modelling/test_geological_model.py index c09405e9..09e4cf97 100644 --- a/tests/unit/modelling/test_geological_model.py +++ b/tests/unit/modelling/test_geological_model.py @@ -17,7 +17,7 @@ def test_rescale_model_data(): model.set_model_data(data) # Check that the model data is rescaled to local coordinates expected = data[['X', 'Y', 'Z']].values - bb[None, 0, :] - actual = model.data[['X', 'Y', 'Z']].values + actual = model.prepare_data(model.data)[['X', 'Y', 'Z']].values assert np.allclose(actual, expected, atol=1e-6) def test_access_feature_model(): data, bb = load_claudius()