From 8394557f7668bc3446294d430a8d319b982e978b Mon Sep 17 00:00:00 2001 From: Noelle Cheng Date: Wed, 10 Sep 2025 19:07:12 +0800 Subject: [PATCH 1/6] refactor sorter --- map2loop/project.py | 47 ++++++++++++++++++++++++++++++--------------- map2loop/sorter.py | 40 +++++++++++++++++++------------------- 2 files changed, 52 insertions(+), 35 deletions(-) diff --git a/map2loop/project.py b/map2loop/project.py index 2e1b46e6..37919482 100644 --- a/map2loop/project.py +++ b/map2loop/project.py @@ -564,15 +564,23 @@ def calculate_stratigraphic_order(self, take_best=False): f"Calculating best stratigraphic column from {[sorter.sorter_label for sorter in sorters]}" ) - columns = [ - sorter.sort( - self.stratigraphic_column.stratigraphicUnits, - self.topology.get_unit_unit_relationships(), - self.contact_extractor.contacts, - self.map_data, - ) - for sorter in sorters - ] + columns = [] + for sorter in sorters: + if sorter.sorter_label == "SorterObservationProjections": + columns.append(sorter.sort( + self.stratigraphic_column.stratigraphicUnits, + self.topology.get_unit_unit_relationships(), + self.contact_extractor.contacts, + self.map_data.get_map_data(Datatype.GEOLOGY), + self.map_data.get_map_data(Datatype.STRUCTURE), + self.map_data.get_map_data(Datatype.DTM), + )) + else: + columns.append(sorter.sort( + self.stratigraphic_column.stratigraphicUnits, + self.topology.get_unit_unit_relationships(), + self.contact_extractor.contacts, + )) basal_contacts = [ self.contact_extractor.extract_basal_contacts( column, save_contacts=False @@ -597,12 +605,21 @@ def calculate_stratigraphic_order(self, take_best=False): self.stratigraphic_column.column = column else: logger.info(f'Calculating stratigraphic column using sorter {self.sorter.sorter_label}') - self.stratigraphic_column.column = self.sorter.sort( - self.stratigraphic_column.stratigraphicUnits, - self.topology.get_unit_unit_relationships(), - self.contact_extractor.contacts, - self.map_data, - ) + if self.sorter.sorter_label == "SorterObservationProjections": + self.stratigraphic_column.column = self.sorter.sort( + self.stratigraphic_column.stratigraphicUnits, + self.topology.get_unit_unit_relationships(), + self.contact_extractor.contacts, + self.map_data.get_map_data(Datatype.GEOLOGY), + self.map_data.get_map_data(Datatype.STRUCTURE), + self.map_data.get_map_data(Datatype.DTM), + ) + else: + self.stratigraphic_column.column = self.sorter.sort( + self.stratigraphic_column.stratigraphicUnits, + self.topology.get_unit_unit_relationships(), + self.contact_extractor.contacts, + ) @beartype.beartype def set_thickness_calculator( diff --git a/map2loop/sorter.py b/map2loop/sorter.py index da4dab76..3b0d6ab9 100644 --- a/map2loop/sorter.py +++ b/map2loop/sorter.py @@ -3,9 +3,11 @@ import pandas import numpy as np import math -from .mapdata import MapData from typing import Union +from osgeo import gdal +import geopandas +from map2loop.utils import value_from_raster from .logging import getLogger logger = getLogger(__name__) @@ -41,7 +43,6 @@ def sort( units: pandas.DataFrame, unit_relationships: pandas.DataFrame, contacts: pandas.DataFrame, - map_data: MapData, ) -> list: """ Execute sorter method (abstract method) @@ -50,7 +51,6 @@ def sort( units (pandas.DataFrame): the data frame to sort (columns must contain ["layerId", "name", "minAge", "maxAge", "group"]) units_relationships (pandas.DataFrame): the relationships between units (columns must contain ["Index1", "Unitname1", "Index2", "Unitname2"]) contacts (pandas.DataFrame): unit contacts with length of the contacts in metres - map_data (map2loop.MapData): a catchall so that access to all map data is available Returns: list: sorted list of unit names @@ -75,7 +75,6 @@ def sort( units: pandas.DataFrame, unit_relationships: pandas.DataFrame, contacts: pandas.DataFrame, - map_data: MapData, ) -> list: """ Execute sorter method takes unit data, relationships and a hint and returns the sorted unit names based on this algorithm. @@ -84,7 +83,6 @@ def sort( units (pandas.DataFrame): the data frame to sort units_relationships (pandas.DataFrame): the relationships between units contacts (pandas.DataFrame): unit contacts with length of the contacts in metres - map_data (map2loop.MapData): a catchall so that access to all map data is available Returns: list: the sorted unit names @@ -140,7 +138,6 @@ def sort( units: pandas.DataFrame, unit_relationships: pandas.DataFrame, contacts: pandas.DataFrame, - map_data: MapData, ) -> list: """ Execute sorter method takes unit data, relationships and a hint and returns the sorted unit names based on this algorithm. @@ -150,7 +147,6 @@ def sort( units_relationships (pandas.DataFrame): the relationships between units stratigraphic_order_hint (list): a list of unit names to use as a hint to sorting the units contacts (pandas.DataFrame): unit contacts with length of the contacts in metres - map_data (map2loop.MapData): a catchall so that access to all map data is available Returns: list: the sorted unit names @@ -192,7 +188,6 @@ def sort( units: pandas.DataFrame, unit_relationships: pandas.DataFrame, contacts: pandas.DataFrame, - map_data: MapData, ) -> list: """ Execute sorter method takes unit data, relationships and a hint and returns the sorted unit names based on this algorithm. @@ -202,7 +197,6 @@ def sort( units_relationships (pandas.DataFrame): the relationships between units stratigraphic_order_hint (list): a list of unit names to use as a hint to sorting the units contacts (pandas.DataFrame): unit contacts with length of the contacts in metres - map_data (map2loop.MapData): a catchall so that access to all map data is available Returns: list: the sorted unit names @@ -280,7 +274,6 @@ def sort( units: pandas.DataFrame, unit_relationships: pandas.DataFrame, contacts: pandas.DataFrame, - map_data: MapData, ) -> list: """ Execute sorter method takes unit data, relationships and a hint and returns the sorted unit names based on this algorithm. @@ -290,7 +283,6 @@ def sort( units_relationships (pandas.DataFrame): the relationships between units stratigraphic_order_hint (list): a list of unit names to use as a hint to sorting the units contacts (pandas.DataFrame): unit contacts with length of the contacts in metres - map_data (map2loop.MapData): a catchall so that access to all map data is available Returns: list: the sorted unit names @@ -362,7 +354,9 @@ def sort( units: pandas.DataFrame, unit_relationships: pandas.DataFrame, contacts: pandas.DataFrame, - map_data: MapData, + geology_data: geopandas.GeoDataFrame, + structure_data: geopandas.GeoDataFrame, + dtm_data: gdal.Dataset ) -> list: """ Execute sorter method takes unit data, relationships and a hint and returns the sorted unit names based on this algorithm. @@ -372,7 +366,9 @@ def sort( units_relationships (pandas.DataFrame): the relationships between units stratigraphic_order_hint (list): a list of unit names to use as a hint to sorting the units contacts (pandas.DataFrame): unit contacts with length of the contacts in metres - map_data (map2loop.MapData): a catchall so that access to all map data is available + geology_data (geopandas.GeoDataFrame): the geology data + structure_data (geopandas.GeoDataFrame): the structure data + dtm_data (ggdal.Dataset): the dtm data Returns: list: the sorted unit names @@ -382,9 +378,15 @@ def sort( from shapely.geometry import LineString, Point from map2loop.m2l_enums import Datatype - geol = map_data.get_map_data(Datatype.GEOLOGY).copy() - geol = geol.drop(geol.index[np.logical_or(geol["INTRUSIVE"], geol["SILL"])]) - orientations = map_data.get_map_data(Datatype.STRUCTURE).copy() + geol = geology_data.copy() + if "INTRUSIVE" in geol.columns: + geol = geol.drop(geol.index[geol["INTRUSIVE"]]) + if "SILL" in geol.columns: + geol = geol.drop(geol.index[geol["SILL"]]) + orientations = structure_data.copy() + dtm = dtm_data.copy() + inv_geotransform = gdal.InvGeoTransform(dtm.GetGeoTransform()) + dtm_array = np.array(dtm.GetRasterBand(1).ReadAsArray().T) # Create a map of maps to store younger/older observations ordered_unit_observations = [] @@ -434,11 +436,9 @@ def sort( continue # Get heights for intersection point and start of ray - height = map_data.get_value_from_raster(Datatype.DTM, start.x, start.y) + height = value_from_raster(inv_geotransform, dtm_array, start.x, start.y) first_intersect_point = Point(start.x, start.y, height) - height = map_data.get_value_from_raster( - Datatype.DTM, second_intersect_point.x, second_intersect_point.y - ) + height = value_from_raster(inv_geotransform, dtm_array, second_intersect_point.x, second_intersect_point.y) second_intersect_point = Point(second_intersect_point.x, start.y, height) # Check vertical difference between points and compare to projected dip angle From 5d04aae9843c61dc11c721abf779f0e1de94a035 Mon Sep 17 00:00:00 2001 From: Noelle Cheng Date: Fri, 12 Sep 2025 14:25:13 +0800 Subject: [PATCH 2/6] fix sorter abstract class sort method to include optional parameters --- map2loop/sorter.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/map2loop/sorter.py b/map2loop/sorter.py index 3b0d6ab9..c822e3ef 100644 --- a/map2loop/sorter.py +++ b/map2loop/sorter.py @@ -43,6 +43,9 @@ def sort( units: pandas.DataFrame, unit_relationships: pandas.DataFrame, contacts: pandas.DataFrame, + geology_data: geopandas.GeoDataFrame = None, + structure_data: geopandas.GeoDataFrame = None, + dtm_data: gdal.Dataset = None, ) -> list: """ Execute sorter method (abstract method) @@ -51,6 +54,9 @@ def sort( units (pandas.DataFrame): the data frame to sort (columns must contain ["layerId", "name", "minAge", "maxAge", "group"]) units_relationships (pandas.DataFrame): the relationships between units (columns must contain ["Index1", "Unitname1", "Index2", "Unitname2"]) contacts (pandas.DataFrame): unit contacts with length of the contacts in metres + geology_data (geopandas.GeoDataFrame): the geology data + structure_data (geopandas.GeoDataFrame): the structure data + dtm_data (ggdal.Dataset): the dtm data Returns: list: sorted list of unit names @@ -75,6 +81,9 @@ def sort( units: pandas.DataFrame, unit_relationships: pandas.DataFrame, contacts: pandas.DataFrame, + geology_data: geopandas.GeoDataFrame = None, + structure_data: geopandas.GeoDataFrame = None, + dtm_data: gdal.Dataset = None, ) -> list: """ Execute sorter method takes unit data, relationships and a hint and returns the sorted unit names based on this algorithm. @@ -138,6 +147,9 @@ def sort( units: pandas.DataFrame, unit_relationships: pandas.DataFrame, contacts: pandas.DataFrame, + geology_data: geopandas.GeoDataFrame = None, + structure_data: geopandas.GeoDataFrame = None, + dtm_data: gdal.Dataset = None, ) -> list: """ Execute sorter method takes unit data, relationships and a hint and returns the sorted unit names based on this algorithm. @@ -188,6 +200,9 @@ def sort( units: pandas.DataFrame, unit_relationships: pandas.DataFrame, contacts: pandas.DataFrame, + geology_data: geopandas.GeoDataFrame = None, + structure_data: geopandas.GeoDataFrame = None, + dtm_data: gdal.Dataset = None, ) -> list: """ Execute sorter method takes unit data, relationships and a hint and returns the sorted unit names based on this algorithm. @@ -274,6 +289,9 @@ def sort( units: pandas.DataFrame, unit_relationships: pandas.DataFrame, contacts: pandas.DataFrame, + geology_data: geopandas.GeoDataFrame = None, + structure_data: geopandas.GeoDataFrame = None, + dtm_data: gdal.Dataset = None, ) -> list: """ Execute sorter method takes unit data, relationships and a hint and returns the sorted unit names based on this algorithm. From b83252a8fcb2abc99fb76859d9e2361da714b777 Mon Sep 17 00:00:00 2001 From: Noelle Cheng Date: Fri, 12 Sep 2025 14:52:24 +0800 Subject: [PATCH 3/6] refactor calculate_stratigraphic_order in project class --- map2loop/project.py | 53 +++++++++++++++++---------------------------- 1 file changed, 20 insertions(+), 33 deletions(-) diff --git a/map2loop/project.py b/map2loop/project.py index 37919482..73b741ab 100644 --- a/map2loop/project.py +++ b/map2loop/project.py @@ -564,23 +564,17 @@ def calculate_stratigraphic_order(self, take_best=False): f"Calculating best stratigraphic column from {[sorter.sorter_label for sorter in sorters]}" ) - columns = [] - for sorter in sorters: - if sorter.sorter_label == "SorterObservationProjections": - columns.append(sorter.sort( - self.stratigraphic_column.stratigraphicUnits, - self.topology.get_unit_unit_relationships(), - self.contact_extractor.contacts, - self.map_data.get_map_data(Datatype.GEOLOGY), - self.map_data.get_map_data(Datatype.STRUCTURE), - self.map_data.get_map_data(Datatype.DTM), - )) - else: - columns.append(sorter.sort( - self.stratigraphic_column.stratigraphicUnits, - self.topology.get_unit_unit_relationships(), - self.contact_extractor.contacts, - )) + columns = [ + sorter.sort( + self.stratigraphic_column.stratigraphicUnits, + self.topology.get_unit_unit_relationships(), + self.contact_extractor.contacts, + self.map_data.get_map_data(Datatype.GEOLOGY), + self.map_data.get_map_data(Datatype.STRUCTURE), + self.map_data.get_map_data(Datatype.DTM), + ) + for sorter in sorters + ] basal_contacts = [ self.contact_extractor.extract_basal_contacts( column, save_contacts=False @@ -605,21 +599,14 @@ def calculate_stratigraphic_order(self, take_best=False): self.stratigraphic_column.column = column else: logger.info(f'Calculating stratigraphic column using sorter {self.sorter.sorter_label}') - if self.sorter.sorter_label == "SorterObservationProjections": - self.stratigraphic_column.column = self.sorter.sort( - self.stratigraphic_column.stratigraphicUnits, - self.topology.get_unit_unit_relationships(), - self.contact_extractor.contacts, - self.map_data.get_map_data(Datatype.GEOLOGY), - self.map_data.get_map_data(Datatype.STRUCTURE), - self.map_data.get_map_data(Datatype.DTM), - ) - else: - self.stratigraphic_column.column = self.sorter.sort( - self.stratigraphic_column.stratigraphicUnits, - self.topology.get_unit_unit_relationships(), - self.contact_extractor.contacts, - ) + self.stratigraphic_column.column = self.sorter.sort( + self.stratigraphic_column.stratigraphicUnits, + self.topology.get_unit_unit_relationships(), + self.contact_extractor.contacts, + self.map_data.get_map_data(Datatype.GEOLOGY), + self.map_data.get_map_data(Datatype.STRUCTURE), + self.map_data.get_map_data(Datatype.DTM), + ) @beartype.beartype def set_thickness_calculator( @@ -1182,4 +1169,4 @@ def save_geotiff_raster( target_ds.SetProjection(geol.crs.to_string()) target_ds = None source_layer = None - source_ds = None + source_ds = None \ No newline at end of file From b5767ef5660f14365c044583a50fc6110619e58c Mon Sep 17 00:00:00 2001 From: noellehmcheng <143368485+noellehmcheng@users.noreply.github.com> Date: Fri, 12 Sep 2025 07:20:00 +0000 Subject: [PATCH 4/6] style: style fixes by ruff and autoformatting by black --- map2loop/sorter.py | 1 - 1 file changed, 1 deletion(-) diff --git a/map2loop/sorter.py b/map2loop/sorter.py index c822e3ef..bcb8b56f 100644 --- a/map2loop/sorter.py +++ b/map2loop/sorter.py @@ -394,7 +394,6 @@ def sort( import networkx as nx import networkx.algorithms.approximation as nx_app from shapely.geometry import LineString, Point - from map2loop.m2l_enums import Datatype geol = geology_data.copy() if "INTRUSIVE" in geol.columns: From a5c272869dbe8005a1584555deb18a381f073630 Mon Sep 17 00:00:00 2001 From: Noelle Cheng Date: Fri, 12 Sep 2025 15:44:10 +0800 Subject: [PATCH 5/6] add libjpeg-turbo installation to yml file --- .github/workflows/linting_and_testing.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/linting_and_testing.yml b/.github/workflows/linting_and_testing.yml index c432d6ce..0259ccb5 100644 --- a/.github/workflows/linting_and_testing.yml +++ b/.github/workflows/linting_and_testing.yml @@ -53,7 +53,7 @@ jobs: if: ${{ matrix.os == 'windows-latest' && (matrix.python-version == '3.9' || matrix.python-version == '3.10') }} run: | conda run -n test conda info - conda run -n test conda install -c loop3d -c conda-forge "gdal=3.4.3" python=${{ matrix.python-version }} -y + conda run -n test conda install -c loop3d -c conda-forge "gdal=3.4.3" python=${{ matrix.python-version }} libjpeg-turbo -y conda run -n test conda install -c loop3d -c conda-forge --file dependencies.txt python=${{ matrix.python-version }} -y conda run -n test conda install pytest python=${{ matrix.python-version }} -y @@ -61,7 +61,7 @@ jobs: if: ${{ matrix.os == 'windows-latest' && matrix.python-version == '3.11' }} run: | conda run -n test conda info - conda run -n test conda install -c loop3d -c conda-forge gdal "netcdf4=*=nompi_py311*" python=${{ matrix.python-version }} -y + conda run -n test conda install -c loop3d -c conda-forge gdal "netcdf4=*=nompi_py311*" python=${{ matrix.python-version }} libjpeg-turbo -y conda run -n test conda install -c loop3d -c conda-forge --file dependencies.txt python=${{ matrix.python-version }} -y conda run -n test conda install pytest python=${{ matrix.python-version }} -y @@ -69,7 +69,7 @@ jobs: if: ${{ !(matrix.os == 'windows-latest' && (matrix.python-version == '3.9' || matrix.python-version == '3.10' || matrix.python-version == '3.11')) }} run: | conda run -n test conda info - conda run -n test conda install -c loop3d -c conda-forge python=${{ matrix.python-version }} gdal pytest --file dependencies.txt -y + conda run -n test conda install -c loop3d -c conda-forge python=${{ matrix.python-version }} gdal pytest libjpeg-turbo --file dependencies.txt -y - name: Install map2loop run: | From c1984e1ed43f63607cd8fb55ec06888c39e46913 Mon Sep 17 00:00:00 2001 From: Noelle Cheng Date: Mon, 15 Sep 2025 20:22:38 +0800 Subject: [PATCH 6/6] fix dtm in SorterObservationProjections --- map2loop/sorter.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/map2loop/sorter.py b/map2loop/sorter.py index bcb8b56f..0428ace6 100644 --- a/map2loop/sorter.py +++ b/map2loop/sorter.py @@ -401,9 +401,8 @@ def sort( if "SILL" in geol.columns: geol = geol.drop(geol.index[geol["SILL"]]) orientations = structure_data.copy() - dtm = dtm_data.copy() - inv_geotransform = gdal.InvGeoTransform(dtm.GetGeoTransform()) - dtm_array = np.array(dtm.GetRasterBand(1).ReadAsArray().T) + inv_geotransform = gdal.InvGeoTransform(dtm_data.GetGeoTransform()) + dtm_array = np.array(dtm_data.GetRasterBand(1).ReadAsArray().T) # Create a map of maps to store younger/older observations ordered_unit_observations = []