Configuring map2loop parameters
- There are four types of parameters we need to define in order to use map2loop:
a) The paths or URLs that tell map2loop where the information layers (GIS files or online sources) are stored. These are passed to the map2loop update_config() method.
b) The names of fields and some text flags that tell map2loop where specific information can be retrieved from these layers. These are stored in an hjson format text file the path of which is passed to the map2loop Project() method.
c) A data output path, bounding box and Coordinate Reference System information to define the extent of the model and where to put it. This is passed to the map2loop update_config() method
d) Parameters that control the specific functioning of map2loop: what to calculate, what decimation factors to apply to the augmented data outputs, what resolution interpolations to use etc. These is passed to the map2looop project.run() method.
a) paths (via project())
These is passed to the map2loop Project() method.
Examples:
Remote WFS layers: See Example 1 notebook
model_name = "./proj_output"
proj = Project(
geology_filename="./source_data/geol.shp",
fault_filename="./source_data/faults_50.shp",
fold_filename="./source_data/fold_axis.shp",
structure_filename="./source_data/strike_dip.shp",
mindep_filename="./source_data/null_mindeps.shp",
metadata_filename='./source_data/m2l_config.hjson',
dtm_filename='./DTM/projected_dtm.tif',
working_projection="EPSG:27700",
verbose_level=VerboseLevel.NONE,
overwrite="true",
project_path=model_name
)
Remote WFS for GSWA with interactive selection of region of interest: See Example 2 notebook
proj=Project(loopdata_state="WA")
where remote=True signifies that WFS-served data will be accessed.
Standard Australia State Geological Surveys datasets, we have predefined the paths for all data and the following code is sufficient: See Example 2 Notebook, do not define if not using standard datasets
proj = Project(
loopdata_state="WA", # choice between 'WA','QLD','NT','NSW','VIC','SA', 'TAS'
)
Local GIS layers: See Example 3 notebook
model_name = "./proj_output"
proj = Project(
geology_filename="./source_data/geol.shp",
fault_filename="./source_data/faults_50.shp",
fold_filename="./source_data/fold_axis.shp",
structure_filename="./source_data/strike_dip.shp",
mindep_filename="./source_data/null_mindeps.shp",
metadata_filename='./source_data/m2l_config.hjson',
dtm_filename='./DTM/projected_dtm.tif',
working_projection="EPSG:27700",
verbose_level=VerboseLevel.NONE,
overwrite="true",
project_path=model_name
)
where remote=False signifies that local GIS files will be accessed. Paths can be relative or absolute, or even a URL, however for URLs, the components of the shapefile or TAB file have to be zipped up.
b) Layer field codes:
You will need to create or modify an hjson format file that provides the names of fields and some text flags that tell map2loop where and what specific information can be retrieved from these layers. These are stored in an hjson format text file the path of which is passed to the map2loop Project() method. The easiest way to get started is to use a jupyter notebook allows you to reduce errors by providing a primitive GUI for creating an hjson config file and associated python script, named: See Utility 1 - Config file generator.ipynb notebook. Alternatively if you are brave you can edit the values to the right of the colon in each row of an existing hjson file. For example to specify that the field in the geospatial layer that contains bedding dip information is called MYDIP, replace the appropriate code in the hjson file below with:
“d”:”MYDIP”,
Some verification is carried out by map2loop to ensure the required parameters have been defined. In the following section field refers to a field name in a geospatial layer; text refers to some text in the contents of a field for a specific geometric object. You shouldn’t use the same field for different codes as this may cause problems.
{
# Orientations-----------------------------
"d": "DIP", # field that contains dip information
"dd": "DIP_DIR", # field that contains dip direction information
"sf": "FEATURE", # field that contains information on type of structure
# text to search for in field defined by sf code to show that this is a bedding measurement
"bedding": "Bed",
# flag to determine measurement convention (currently "strike" or "dip direction")
"otype": "dip direction",
"bo": "TYPE", # field that contains type of foliation
# text to search for in field defined by bo code to show that this is an overturned bedding measurement
"btype": "overturned",
# Stratigraphy-----------------------------
"g": "GROUP_", # field that contains coarser stratigraphic coding
# field that contains alternate coarser stratigraphic coding if "g" is blank
"g2": "SUPERSUITE",
"c": "UNITNAME", # field that contains finer stratigraphic coding
"ds": "DESCRIPTN", # field that contains information about lithology
# field that contains alternate stratigraphic coding (not used??)
"u": "CODE",
"r1": "ROCKTYPE1", # field that contains extra lithology information
"r2": "ROCKTYPE2", # field that contains even more lithology information
"sill": "sill", # text to search for in field defined by ds code to show that this is a sill
# text to search for in field defined by r1 code to show that this is an intrusion
"intrusive": "intrusive", # text to search for in field defined by ds code to show that this is an volcanic (not intrusion) "volcanic": "volcanic",
# Mineral Deposits-----------------------------
"msc": "SITE_CODE", # field that contains site code of deposit
"msn": "SHORT_NAME", # field that contains short name of deposit
"mst": "SITE_TYPE_", # field that contains site type of deposit
"mtc": "TARGET_COM", # field that contains target commodity of deposit
"mscm": "SITE_COMMO", # field that contains site commodity of deposit
"mcom": "COMMODITY_", # field that contains commodity group of deposit
# text to search for in field defined by mst code that shows site to ignore
"minf": "Infrastructure",
# Timing-----------------------------
"min": "MIN_AGE_MA", # field that contains minimum age of unit defined by ccode
"max": "MAX_AGE_MA", # field that contains maximum age of unit defined by ccode
#faults and folds-----------------------------
"f": "FEATURE", # field that contains information on type of structure
# text to search for in field defined by f code to show that this is a fault
"fault": "Fault",
"ff": "FEATURE", # field that contains information on type of structure
# text to search for in field defined by f code to show that this is a fold axial trace
"fold": "Fold axial trace",
"fdip": "DIP", # field for numeric fault dip value
# text to search for in field defined by fdip to show that this has no known dip
"fdipnull": "0",
"fdipdir": "DIP_DIR", # field for text fault dip direction value
# flag for text fault dip direction type num e.g. 045 or alpha e.g. southeast
"fdipdir_flag": "alpha",
"fdipest": "DIP_EST", # field for text fault dip estimate value
# text to search for in field defined by fdipest to give fault dip estimate in increasing steepness
"fdipest_vals": "gentle,moderate,steep",
# field that contains information on name of fault (not used??)
"n": "NAME",
"t": "TYPE", # field that contains information on type of fold
# text to search for in field defined by t to show that this is a syncline
"syn": "syncline",
# ids-----------------------------
"o": "OBJECTID", # field that contains unique id of geometry object
"gi": "GEOPNT_ID", # field that contains unique id of structure point
"deposit_dist": 500
}
c) ROI, Projection, output paths (via update_config())
A data output path which points to a new or existing directory (a new directory will be created if needed), bounding box and Coordinate Reference System information to define the extent of the model. This is be passed to the map2loop update_config() method
proj.update_config(
clut_path='./source_data/Wales_Colour.csv',
bbox_3d={
"minx": 324241,
#"minx": 407000,
"miny": 209276,
"maxx": 328607,
"maxy": 213331,
#"maxy": 370000,
"base": 0,
"top": 600,
},
run_flags={
'aus': False,
'close_dip': -999,
'contact_decimate': 25,
'cover_dip':10,
'cover_spacing':10000,
'contact_dip': -999,
'contact_orientation_decimate': 5,
'deposits': "Fe,Cu,Au,NONE",
'dist_buffer': 10,
'dtb_nul':None,
'fat_step': 750,
'fault_decimate': 0, # only a few points per fault, so no need to decimate
'fault_dip': 90,
'fold_decimate': 5,
'interpolation_scheme': 'scipy_rbf',
'interpolation_spacing': 500,
'intrusion_mode': 0,
'max_thickness_allowed': 10000,
'min_fault_length': 500,
'misorientation': 30,
'null_scheme': 'null',
'orientation_decimate': 0, # only 75 points so no need to decimate
'pluton_dip': 45,
'pluton_form': 'saucers',
'thickness_buffer': 5000,
'use_fat': False,
'use_interpolations': False
},
drift_prefix=['MCQGB']
)
Full list of update_config flags:
bbox_3d 3D bounding box of coordinates and base/top values defining the area, defaults to { “minx”: 0, “maxx”: 0, “maxx”: 0, “maxy”: 0, “base”: -10000, “top”: 1200, } :type bbox_3d: dict, optional
clut_path Path to custom map colours file :type clut_path: string, optional
run_flags Global dictionary that defines custom params such as decimation and minimum fault length, see below :type run_flags: dict, optional
**kwargs
Calculation control parameters (via run_flags dictionary)
These control the specific functionality of map2loop: what to calculate, what decimation factors to apply to the augmented data outputs, what resolution interpolations to use etc. These are passed to the map2looop project run() method:
proj.run()
This method performs the data processing steps of the map2loop workflow, and can be modified by including the following parameters as run_flags [defaults](data type):
aus: Indicates if area is in Australia for using ASUD, the Australian Stratigraphic Units Database to redfine stratigraphic relationships. Should only be True in Australia, and when the finest stratigraphic level is the ASUD standard Formation name. [True] (bool)
close_dip: Dip to assign to all new fold axial trace orientations. If -999 then the nearest interpolated dip for that supergroup will be used instead. [-999] In degrees (int)
contact_decimate: Save every nth contact data point. 0 means save all data. [5] (int)
contact_dip: Dip to assign to all new basal contact orientations. If -999 then the nearest interpolated dip for that supergroup will be used instead. [-999] In degrees (int)
contact_orientation_decimate: Save every nth contact orientation point. 0 means save all data. [5] (int)
deposits: Mineral deposit names for focused topology extraction. [“Fe,Cu,Au,NONE”] Topological analysis of faults and strat will only be carried out relative to these deposit type. NONE must always be one of the types (str)
dist_buffer: Buffer for processing basal contacts. Basal contact vertices less than this distance from the fault will be ignored. [10] In metres. (int)
dtb: Path to depth to basement grid. Geotif of depths in the same projection system as everything else. [‘’] (str)
fat_step: How much to step out normal to the fold axial trace. Distance in metres. [750] In metres. (int)
fault_decimate: Save every nth fault data point along fault tace. 0 means save all data. [5] (int)
fault_dip: default fault dip [90] In degrees (int)
fold_decimate: Save every nth fold axial trace data point. 0 means save all data. [5] (int)
interpolation_scheme: What interpolation method to use of scipy_rbf (radial basis) or scipy_idw (inverse distance weighted). [‘scipy_rbf’] (str)
interpolation_spacing: Interpolation grid spacing in meters. Used to interpolation bedding orientations [500] In metres or if a negative value defines fixed number of grid points in x & y (int)
intrusion_mode: 1 to exclude all intrusions from basal contacts, [0] to only exclude sills. [0] (int)
max_thickness_allowed: when estimating local formation thickness [10000] in metres. (int)
min_fault_length: Min fault length to be considered. In metres. [5000] In meters. (int)
misorientation: [30] Maximum misorientation in pole to great circle of bedding between groups to be considered part of same supergroup (int)
null_scheme: How null values present in the depth to basement geotif. [‘null’] (str)
orientation_decimate: Save every nth orientation data point. 0 means save all data. [0] type int
pluton_dip: default pluton contact dip [45] In degrees (int)
pluton_form: Possible forms from domes, saucers or pendant. [‘domes’] (str)
thickness_buffer: How far away to look for next highest unit when calculating formation thickness [5000] In metres. (int)
use_fat: Use fold axial trace info to add near-axis bedding info [True] (bool)
use_interpolations: Use all interpolated dips for modelling [True] (bool)
fault_orientation_clusters:[2] number of clusters for kmeans clustering of faults by orientation (int)
fault_length_clusters: number of clusters for kmeans clustering of faults by length [2] (int)
use_roi_clip: use non-rectangular ROI polygon [False] (bool)
roi_clip_path: path to non-rectangular ROI polygon shapefile [‘’] (bool)
d) Calculation workflow parameters
seismic_section: Add data from a single seismic section (paths hardwired for the moment) [False] (bool)
cover_map: Add data from a depth to basement grid (paths hardwired for the moment) [False] (bool)
near_fault_interpolations: Add stratigraphic info near faults [False] (bool)
fold_axial_traces: Add dip info either side of fold axial trace to enhance second order folds [False] (bool)
stereonets: Calculate stereonets to define supergroups [True] (bool)
formation_thickness: Calculate formation thickness [True] (bool)
polarity: Calculate bedding polarity (doesn’t work!!) [False] (bool)
strat_offset: Calculate stratigraphic offset across faults [True] (bool)
contact_dips: Add fixed or interpolated dips to contacts [True] (bool)
Individual workflow parameters can be overwritten AFTER the call to proj.update_config() as follows:
proj.workflow['contact_dips'] = False
Example minimum code
An example minimum code to run map2loop with mostly default settings might look like this See Example 3 notebook ):
import os
from map2loop.project import Project
from map2loop.m2l_enums import VerboseLevel
import time
from datetime import datetime
nowtime=datetime.now().isoformat(timespec='minutes')
model_name=nowtime.replace(":","-").replace("T","-")
t0 = time.time()
proj = Project(
geology_filename="./source_data/geol_clip.shp",
fault_filename="./source_data/faults_clip.shp",
fold_filename="./source_data/folds_clip.shp",
structure_filename="./source_data/structure_clip.shp",
mindep_filename="./source_data/mindeps_clip.shp",
dtm_filename='./source_data/dtm_rp.tif',
metadata_filename='./source_data/example.hjson',
overwrite="true",
verbose_level=VerboseLevel.NONE,
project_path=model_name,
working_projection="EPSG:28350",
)
proj.update_config(
out_dir=model_name,
bbox_3d={
"minx": 520000,
"miny": 7490000,
"maxx": 550000,
"maxy": 7510000,
"base": -3200,
"top": 1200,
},
run_flags={ },
#proj_crs={'init': 'EPSG:28350'},
clut_path='./source_data/500kibg_colours.csv',
quiet='all' # change this to 'None' (with quotes) to see intermediate output
)
proj.run()
This code will take information from the GSWA 1:500K Interpeted Bedrock Geology map:
which generates a series of augmented outputs, including: a topological analysis of the spatial and temporal relationships of fetaures in the map:
extraction of basal contacts of units, fault locations, orientaion of bedding and foliations, and derived products such as local formation thickness estimates:
these can then be fed into 3D modelling packages such as LoopStructural (top Link to interactive 3D model) and gempy (lower).
