|
| 1 | +import geopandas |
| 2 | +import pandas |
| 3 | +import shapely |
| 4 | +from .logging import getLogger |
| 5 | +from typing import Optional |
| 6 | + |
| 7 | +logger = getLogger(__name__) |
| 8 | + |
| 9 | +class ContactExtractor: |
| 10 | + def __init__(self, geology: geopandas.GeoDataFrame, faults: Optional[geopandas.GeoDataFrame] = None): |
| 11 | + self.geology = geology |
| 12 | + self.faults = faults |
| 13 | + self.contacts = None |
| 14 | + self.basal_contacts = None |
| 15 | + self.all_basal_contacts = None |
| 16 | + |
| 17 | + def extract_all_contacts(self, save_contacts: bool = True) -> geopandas.GeoDataFrame: |
| 18 | + logger.info("Extracting contacts") |
| 19 | + geology = self.geology.copy() |
| 20 | + geology = geology.dissolve(by="UNITNAME", as_index=False) |
| 21 | + geology = geology[~geology["INTRUSIVE"]] |
| 22 | + geology = geology[~geology["SILL"]] |
| 23 | + if self.faults is not None: |
| 24 | + faults = self.faults.copy() |
| 25 | + faults["geometry"] = faults.buffer(50) |
| 26 | + geology = geopandas.overlay(geology, faults, how="difference", keep_geom_type=False) |
| 27 | + units = geology["UNITNAME"].unique().tolist() |
| 28 | + column_names = ["UNITNAME_1", "UNITNAME_2", "geometry"] |
| 29 | + contacts = geopandas.GeoDataFrame(crs=geology.crs, columns=column_names, data=None) |
| 30 | + while len(units) > 1: |
| 31 | + unit1 = units[0] |
| 32 | + units = units[1:] |
| 33 | + for unit2 in units: |
| 34 | + if unit1 != unit2: |
| 35 | + join = geopandas.overlay( |
| 36 | + geology[geology["UNITNAME"] == unit1], |
| 37 | + geology[geology["UNITNAME"] == unit2], |
| 38 | + keep_geom_type=False, |
| 39 | + )[column_names] |
| 40 | + join["geometry"] = join.buffer(1) |
| 41 | + buffered = geology[geology["UNITNAME"] == unit2][["geometry"]].copy() |
| 42 | + buffered["geometry"] = buffered.boundary |
| 43 | + end = geopandas.overlay(buffered, join, keep_geom_type=False) |
| 44 | + if len(end): |
| 45 | + contacts = pandas.concat([contacts, end], ignore_index=True) |
| 46 | + contacts["length"] = [row.length for row in contacts["geometry"]] |
| 47 | + if save_contacts: |
| 48 | + self.contacts = contacts |
| 49 | + return contacts |
| 50 | + |
| 51 | + def extract_basal_contacts(self, |
| 52 | + stratigraphic_column: list, |
| 53 | + save_contacts: bool = True) -> geopandas.GeoDataFrame: |
| 54 | + |
| 55 | + logger.info("Extracting basal contacts") |
| 56 | + units = stratigraphic_column |
| 57 | + |
| 58 | + if self.contacts is None: |
| 59 | + self.extract_all_contacts(save_contacts=True) |
| 60 | + basal_contacts = self.contacts.copy() |
| 61 | + else: |
| 62 | + basal_contacts = self.contacts.copy() |
| 63 | + if any(unit not in units for unit in basal_contacts["UNITNAME_1"].unique()): |
| 64 | + missing_units = ( |
| 65 | + basal_contacts[~basal_contacts["UNITNAME_1"].isin(units)]["UNITNAME_1"] |
| 66 | + .unique() |
| 67 | + .tolist() |
| 68 | + ) |
| 69 | + logger.error( |
| 70 | + "There are units in the Geology dataset, but not in the stratigraphic column: " |
| 71 | + + ", ".join(missing_units) |
| 72 | + + ". Please readjust the stratigraphic column if this is a user defined column." |
| 73 | + ) |
| 74 | + raise ValueError( |
| 75 | + "There are units in stratigraphic column, but not in the Geology dataset: " |
| 76 | + + ", ".join(missing_units) |
| 77 | + + ". Please readjust the stratigraphic column if this is a user defined column." |
| 78 | + ) |
| 79 | + basal_contacts["ID"] = basal_contacts.apply( |
| 80 | + lambda row: min(units.index(row["UNITNAME_1"]), units.index(row["UNITNAME_2"])), axis=1 |
| 81 | + ) |
| 82 | + basal_contacts["basal_unit"] = basal_contacts.apply(lambda row: units[row["ID"]], axis=1) |
| 83 | + basal_contacts["stratigraphic_distance"] = basal_contacts.apply( |
| 84 | + lambda row: abs(units.index(row["UNITNAME_1"]) - units.index(row["UNITNAME_2"])), axis=1 |
| 85 | + ) |
| 86 | + basal_contacts["type"] = basal_contacts.apply( |
| 87 | + lambda row: "ABNORMAL" if abs(row["stratigraphic_distance"]) > 1 else "BASAL", axis=1 |
| 88 | + ) |
| 89 | + basal_contacts = basal_contacts[["ID", "basal_unit", "type", "geometry"]] |
| 90 | + basal_contacts["geometry"] = [ |
| 91 | + shapely.line_merge(shapely.snap(geo, geo, 1)) for geo in basal_contacts["geometry"] |
| 92 | + ] |
| 93 | + if save_contacts: |
| 94 | + self.all_basal_contacts = basal_contacts |
| 95 | + self.basal_contacts = basal_contacts[basal_contacts["type"] == "BASAL"] |
| 96 | + return basal_contacts |
0 commit comments