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: | diff --git a/map2loop/project.py b/map2loop/project.py index 2e1b46e6..73b741ab 100644 --- a/map2loop/project.py +++ b/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 \ No newline at end of file diff --git a/map2loop/sorter.py b/map2loop/sorter.py index da4dab76..0428ace6 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,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,7 +384,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 @@ -380,11 +394,15 @@ 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 = 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