search

Loop3D / map2loop

MIT License
9 stars 9 forks source link

[Bug] - InterpolatedStructure thickness calculator fails for collocated points #123

Closed lachlangrose closed 1 month ago

lachlangrose commented 1 month ago

Version

3.1.8

Bug Description

InterpolatedStructure thickness calculator raises singular matrix error when trying to calculate thickness for the example map2loop notebooks. From a quick look it seems like there are collocated points and it is possible that the RBF doesn't like collocated points.

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,
"minx": 520000,
"miny": 7490000,
"maxx": 550000,
"maxy": 7510000,
"base": -3200,
"top": 1200,

}

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,

working_projection = "EPSG:28350",
bounding_box = bbox_3d,
verbose_level = VerboseLevel.NONE,
loop_project_filename = loop_project_filename,
overwrite_loopprojectfile = True

)

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

'Turee_Creek_Group',
'Boolgeeda_Iron_Formation',
'Woongarra_Rhyolite',
'Weeli_Wolli_Formation',
'Brockman_Iron_Formation',
'Mount_McRae_Shale_and_Mount_Sylvia_Formation',
'Wittenoom_Formation',
'Marra_Mamba_Iron_Formation',
'Jeerinah_Formation',
'Bunjinah_Formation',
'Pyradie_Formation',
'Fortescue_Group',
# oldest

] proj.run_all(user_defined_stratigraphic_column=column)

or

proj.run_all(take_best=True)

t1 = time.time()

Expected Behavior

It should ignore collocated points, or add a small value to the diagonal of the A matrix for the interpolation.

Actual Behavior

exception

Additional Context

No response

Environment

ubuntu

Severity