Fix: refactor sorter (#233)

* refactor sorter

* fix sorter abstract class sort method to include optional parameters

* refactor calculate_stratigraphic_order in project class

* style: style fixes by ruff and autoformatting by black

* add libjpeg-turbo installation to yml file

* fix dtm in SorterObservationProjections

Author: noellehmcheng

.github/workflows/linting_and_testing.yml

@@ -53,23 +53,23 @@ 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
 
       - name: Install dependencies for windows python 3.11
         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
 
       - name: Install dependencies for other environments
         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: |

map2loop/project.py

@@ -569,7 +569,9 @@ def calculate_stratigraphic_order(self, take_best=False):
                     self.stratigraphic_column.stratigraphicUnits,
                     self.topology.get_unit_unit_relationships(),
                     self.contact_extractor.contacts,
-                    self.map_data,
+                    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
             ]
@@ -601,7 +603,9 @@ def calculate_stratigraphic_order(self, take_best=False):
                 self.stratigraphic_column.stratigraphicUnits,
                 self.topology.get_unit_unit_relationships(),
                 self.contact_extractor.contacts,
-                self.map_data,
+                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
@@ -1165,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

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,9 @@ def sort(
         units: pandas.DataFrame,
         unit_relationships: pandas.DataFrame,
         contacts: pandas.DataFrame,
-        map_data: MapData,
+        geology_data: geopandas.GeoDataFrame = None,
+        structure_data: geopandas.GeoDataFrame = None,
+        dtm_data: gdal.Dataset = None,
     ) -> list:
         """
         Execute sorter method (abstract method)
@@ -50,7 +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
-            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: sorted list of unit names
@@ -75,7 +81,9 @@ def sort(
         units: pandas.DataFrame,
         unit_relationships: pandas.DataFrame,
         contacts: pandas.DataFrame,
-        map_data: MapData,
+        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.
@@ -84,7 +92,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 +147,9 @@ def sort(
         units: pandas.DataFrame,
         unit_relationships: pandas.DataFrame,
         contacts: pandas.DataFrame,
-        map_data: MapData,
+        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.
@@ -150,7 +159,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 +200,9 @@ def sort(
         units: pandas.DataFrame,
         unit_relationships: pandas.DataFrame,
         contacts: pandas.DataFrame,
-        map_data: MapData,
+        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.
@@ -202,7 +212,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 +289,9 @@ def sort(
         units: pandas.DataFrame,
         unit_relationships: pandas.DataFrame,
         contacts: pandas.DataFrame,
-        map_data: MapData,
+        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.
@@ -290,7 +301,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 +372,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,19 +384,25 @@ 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
         """
         import networkx as nx
         import networkx.algorithms.approximation as nx_app
         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()
+        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 = []
@@ -434,11 +452,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