StructurePoint thickness calculator fails when the a unit from the stratigraphic column is not included on the map.
Example from the example notebooks
Minimal reproducible example
import time
import os
from map2loop.project import Project
from map2loop.m2l_enums import VerboseLevel, Datatype
from map2loop.sorter import SorterAlpha, SorterAgeBased, SorterUseHint, SorterUseNetworkX, SorterMaximiseContacts, SorterObservationProjections
from map2loop.thickness_calculator import ThicknessCalculatorAlpha, InterpolatedStructure, StructuralPoint
from map2loop.sampler import SamplerSpacing
from datetime import datetime
M2l should skip any units not included on the map. When using m2l for remote data or varying bounding boxes the column may be a larger area than the region.
Actual Behavior
File ~/miniconda3/envs/map2loop-new/lib/python3.11/site-packages/map2loop/thickness_calculator.py:486, in StructuralPoint.compute(self, units, stratigraphic_order, basal_contacts, structure_data, map_data)
483 bbox_poly = geology[geology['UNITNAME'] == litho_in][['minx', 'miny', 'maxx', 'maxy']]
485 # make a subset of the geology polygon & find neighbour units
--> 486 GEO_SUB = geology[geology['UNITNAME'] == litho_in]['geometry'].values[0]
487 neighbor_list = list(
488 basal_contacts[GEO_SUB.intersects(basal_contacts.geometry)]['basal_unit']
489 )
490 # draw orthogonal line to the strike (default value 10Km), and clip it by the bounding box of the lithology
File ~/miniconda3/envs/map2loop-new/lib/python3.11/site-packages/geopandas/array.py:422, in GeometryArray.getitem(self, idx)
420 def getitem(self, idx):
421 if isinstance(idx, numbers.Integral):
--> 422 return self._data[idx]
423 # array-like, slice
424 # validate and convert IntegerArray/BooleanArray
425 # to numpy array, pass-through non-array-like indexers
426 idx = pd.api.indexers.check_array_indexer(self, idx)
IndexError: index 0 is out of bounds for axis 0 with size 0
Version
m2l 3.1.8+, LPF 0.1.3
Bug Description
StructurePoint thickness calculator fails when the a unit from the stratigraphic column is not included on the map.
Example from the example notebooks
Minimal reproducible example
import time import os from map2loop.project import Project from map2loop.m2l_enums import VerboseLevel, Datatype from map2loop.sorter import SorterAlpha, SorterAgeBased, SorterUseHint, SorterUseNetworkX, SorterMaximiseContacts, SorterObservationProjections from map2loop.thickness_calculator import ThicknessCalculatorAlpha, InterpolatedStructure, StructuralPoint from map2loop.sampler import SamplerSpacing from datetime import datetime
os.environ["DOCUMENTATION_TEST"] = "True"
t0 = time.time()
Set the region of interest for the project
bbox_3d = { "minx": 515687.31005864, "miny": 7493446.76593407, "maxx": 562666.860106543, "maxy": 7521273.57407786, "base": -3200, "top": 3000, }
Specify minimum details (Australian state, projection, bounding box and output file)
Optional: the config information manually created in the config_dict can be used by
setting the parameter config_dictionary
loop_project_filename="wa_output.loop3d" proj = Project( use_australian_state_data = "WA",
config_dictionary = config_dict,
)
Specify that the text to look for (in the unit description) to indicate a sill is not just "sill"
proj.map_data.config.update_from_dictionary({"geology": {"sill_text":"is a sill"}})
Set the distance between sample points for arial and linestring geometry
proj.set_sampler(Datatype.GEOLOGY, SamplerSpacing(200.0)) proj.set_sampler(Datatype.FAULT, SamplerSpacing(200.0)) proj.set_minimum_fault_length(5000.0)
Choose which thickness calculator to use
proj.set_thickness_calculator(ThicknessCalculatorAlpha())
proj.set_thickness_calculator(InterpolatedStructure())
proj.set_thickness_calculator(StructuralPoint()) proj.map_data.config.geology_config["intrusive_text"] = "mafic intrusive"
Set specific layers from the geology map to be ignored (commonly "cover" or "water")
proj.set_ignore_codes(["cover", "Fortescue_Group", "A_FO_od"])
Choose which stratigraphic sorter to use or run_all with "take_best" flag to run them all
proj.set_sorter(SorterAlpha())
proj.set_sorter(SorterAgeBased())
proj.set_sorter(SorterUseHint())
proj.set_sorter(SorterUseNetworkX())
proj.set_sorter(SorterMaximiseContacts())
proj.set_sorter(SorterObservationProjections())
column = [
youngest
] proj.run_all(user_defined_stratigraphic_column=column)
or
proj.run_all(take_best=True)
t1 = time.time()
Expected Behavior
M2l should skip any units not included on the map. When using m2l for remote data or varying bounding boxes the column may be a larger area than the region.
Actual Behavior
File ~/miniconda3/envs/map2loop-new/lib/python3.11/site-packages/map2loop/thickness_calculator.py:486, in StructuralPoint.compute(self, units, stratigraphic_order, basal_contacts, structure_data, map_data) 483 bbox_poly = geology[geology['UNITNAME'] == litho_in][['minx', 'miny', 'maxx', 'maxy']] 485 # make a subset of the geology polygon & find neighbour units --> 486 GEO_SUB = geology[geology['UNITNAME'] == litho_in]['geometry'].values[0] 487 neighbor_list = list( 488 basal_contacts[GEO_SUB.intersects(basal_contacts.geometry)]['basal_unit'] 489 ) 490 # draw orthogonal line to the strike (default value 10Km), and clip it by the bounding box of the lithology
File ~/miniconda3/envs/map2loop-new/lib/python3.11/site-packages/geopandas/array.py:422, in GeometryArray.getitem(self, idx) 420 def getitem(self, idx): 421 if isinstance(idx, numbers.Integral): --> 422 return self._data[idx] 423 # array-like, slice 424 # validate and convert IntegerArray/BooleanArray 425 # to numpy array, pass-through non-array-like indexers 426 idx = pd.api.indexers.check_array_indexer(self, idx)
IndexError: index 0 is out of bounds for axis 0 with size 0
Additional Context
No response
Environment
ubuntu
Severity