Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions .github/workflows/linting_and_testing.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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: |
Expand Down
10 changes: 7 additions & 3 deletions map2loop/project.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
]
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
58 changes: 37 additions & 21 deletions map2loop/sorter.py
Original file line number Diff line number Diff line change
Expand Up @@ -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__)
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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.
Expand All @@ -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
Expand Down Expand Up @@ -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.
Expand All @@ -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
Expand Down Expand Up @@ -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.
Expand All @@ -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
Expand Down Expand Up @@ -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.
Expand All @@ -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
Expand Down Expand Up @@ -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.
Expand All @@ -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 = []
Expand Down Expand Up @@ -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
Expand Down
Loading