<jupyter_start><jupyter_text>map2loop: From geology layers to outputs to various mineral system outputs- Hamersley example This notebook reads in three layers from  local or remote sources:  geology polygons, orientation data and fault polylines; and calculates the topological relationships between the different features. Requires compiled cpp code from <NAME>This all gets fed into successive tolopogical and geometric transfroms that end up feeding into geomodeller to make a 3D model<jupyter_code>import geopandas as gpd
import pandas as pd
import numpy as np
import os
import sys
import stat
import functools 
import operator  
import matplotlib
import networkx as nx
import rasterio
from shapely.geometry import Polygon
from map2loop import m2l_utils
from map2loop import m2l_topology
from map2loop import m2l_geometry
from map2loop import m2l_interpolation
from map2loop import m2l_export
import time
%matplotlib inline
t0 = time.time()

#newwd="C:\\Users\\00073294\Dropbox\\loop_minex\\map2model\\"
#os.chdir(newwd)
print("Current Working Directory " )

gdal_data = os.environ['GDAL_DATA']
print("***",gdal_data)
print('is dir: ' + str(os.path.isdir(gdal_data)))
gcs_csv = os.path.join(gdal_data, 'gcs.csv')
print('is file: ' + str(os.path.isfile(gcs_csv)))
st = os.stat(gcs_csv)
print('is readable: ' + str(bool(st.st_mode & stat.S_IRGRP)))
os.environ['PROJ_LIB']=r"C:\\Users\\00073294\\AppData\\Local\\Continuum\\anaconda3\\Lib\\site-packages\\pyproj\\proj_dir\\share\\proj"
print(os.getenv('PROJ_LIB'))
print(os.getcwd())<jupyter_output><empty_output><jupyter_text>Create bounding box based on inputs<jupyter_code>test_data_name='test_minsys'

test_data_path='../'+test_data_name+'/'

os.chdir(test_data_path)
%run -i "m2l_config.py"
#%run -i "m2l_config_remote.py"
print(os.getcwd())

bbox2=str(minx)+","+str(miny)+","+str(maxx)+","+str(maxy)
lat_point_list = [miny, miny, maxy, maxy, maxy]
lon_point_list = [minx, maxx, maxx, minx, minx]
bbox_geom = Polygon(zip(lon_point_list, lat_point_list))
polygon = gpd.GeoDataFrame(index=[0], crs=dst_crs, geometry=[bbox_geom]) 
bbox=(minx,miny,maxx,maxy)


fold_decimate=2         

contact_decimate=5

use_interpolations=True       #use interpolated dips/contacts as additional constraints

use_fat=True                   #use fold axial trace orientation hints

pluton_form='domes'

fault_dip=90

min_fault_length=5000

compute_etc=False

#local_paths=False<jupyter_output><empty_output><jupyter_text>First we test to see if we have access to the online data we need<jupyter_code>loopwfs=m2l_utils.have_access("geo.loop-gis.org")
hawaii=m2l_utils.have_access("pae-paha.pacioos.hawaii.edu")

if(not (loopwfs & hawaii)):
    local_paths=True
    net=False
    print('using local paths')
else:
    net=True<jupyter_output><empty_output><jupyter_text>Optional WFS source  WFS brings in field names as lower case, so need to redefine codes too<jupyter_code>if((not local_paths) and net):
    structure_file='http://geo.loop-gis.org/geoserver/loop/wfs?service=WFS&version=1.1.0&request=GetFeature&typeName=warox_points_f5011&bbox='+bbox2+'&srs=EPSG:28350'
    fault_file='http://geo.loop-gis.org/geoserver/loop/wfs?service=WFS&version=1.1.0&request=GetFeature&typeName=linear_500k&bbox='+bbox2+'&srs=EPSG:28350'
    geology_file='http://geo.loop-gis.org/geoserver/loop/wfs?service=WFS&version=1.0.0&request=GetFeature&typeName=loop:geol_500k&bbox='+bbox2+'&srs=EPSG:28350'

    c_l= {
    #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
      "bedding": 'Bed',            #text to search for in field defined by sfcode to show that this is a bedding measurement
    #Stratigraphy
      "g": 'group_',               #field that contains coarser stratigraphic coding
      "c": 'code',                 #field that contains finer stratigraphic coding
      "ds": 'descriptn',           #field that contains information about lithology
      "u": 'unitname',             #field that contains alternate stratigraphic coding (not used??)
      "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 dscode to show that this is a sill
      "intrusive": 'intrusive',    #text to search for in field defined by dscode to show that this is an intrusion
      "volcanic": 'volcanic',      #text to search for in field defined by dscode to show that this is an intrusion
    #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
      "minf": 'Infrastructure',    #text to search for in field defined by mst code that shows site to ignore
    #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
      "fault": 'Fault',            #text to search for in field defined by fcode to show that this is a fault
      "fold": 'Fold axial trace',  #text to search for in field defined by fcode to show that this is a fold axial trace
      "n": 'name',                 #field that contains information on name of fault (not used??)
      "t": 'type',                 #field that contains information on type of fold
      "syn": 'syncline',           #text to search for in field defined by t to show that this is a syncline
    #ids
      "o": 'objectid',             #field that contains unique id of geometry object
      "gi": 'geopnt_id'            #field that contains unique id of structure point
    }<jupyter_output><empty_output><jupyter_text>Plot geology polygons and bounding box<jupyter_code>print(geology_file)
geology_ll = gpd.read_file(geology_file,bbox=bbox)

base=geology_ll.plot(column=c_l['c'],figsize=(10,10),edgecolor='#000000',linewidth=0.2)
polygon.plot(ax=base, color='none',edgecolor='black')<jupyter_output><empty_output><jupyter_text>Save geology to file as WKT<jupyter_code>sub_geol = geology_ll[['geometry', c_l['o'],c_l['c'],c_l['g'],c_l['u'],c_l['min'],c_l['max'],c_l['ds'],c_l['r1'],c_l['r2']]]
m2l_topology.save_geol_wkt(sub_geol,geology_file_csv, c_l)<jupyter_output><empty_output><jupyter_text>Save mineral deposits to file as WKTThis is not needed by map2loop to build 3D models, but is used by map2model to calculate mineral deposit/topology analyses.<jupyter_code>mindep = gpd.read_file(mindep_file,bbox=bbox)

sub_mindep = mindep[['geometry', c_l['msc'],c_l['msn'],c_l['mst'],c_l['mtc'],c_l['mscm'],c_l['mcom']]]
m2l_topology.save_mindep_wkt(sub_mindep,mindep_file_csv, c_l)

base=sub_mindep.plot()
polygon.plot(ax=base, color='none',edgecolor='black')<jupyter_output><empty_output><jupyter_text>Read and save WAROX point data as WKT<jupyter_code>warox = gpd.read_file(structure_file,bbox=bbox)

sub_pts = warox[['geometry', c_l['gi'],c_l['d'],c_l['dd']]]

m2l_topology.save_structure_wkt(sub_pts,structure_file_csv,c_l)

base=sub_pts.plot()
polygon.plot(ax=base, color='none',edgecolor='black')<jupyter_output><empty_output><jupyter_text>Plot faults and bounding box<jupyter_code>lines_ll=gpd.read_file(fault_file,bbox=bbox)

base2=lines_ll.plot(cmap='rainbow',column=c_l['f'],figsize=(10,10),linewidth=0.4)
polygon.plot(ax=base2, color='none',edgecolor='black')<jupyter_output><empty_output><jupyter_text>Save faults to file as WKT<jupyter_code>sub_lines = lines_ll[['geometry', c_l['o'],c_l['f']]]
m2l_topology.save_faults_wkt(sub_lines,fault_file_csv,c_l)<jupyter_output><empty_output><jupyter_text>Create map2model input file<jupyter_code>m2l_topology.save_Parfile(m2m_cpp_path,c_l,graph_path,geology_file_csv,fault_file_csv,structure_file_csv,minx,maxx,miny,maxy,500.0,'Fe,Cu,Au,NONE')<jupyter_output><empty_output><jupyter_text>Calculate topology<jupyter_code>os.chdir(m2m_cpp_path)
print(os.getcwd())
%system map2model.exe Parfile<jupyter_output><empty_output><jupyter_text>Simple network graph of the geology with legend<jupyter_code>G=nx.read_gml(strat_graph_file,label='id')
selected_nodes = [n for n,v in G.nodes(data=True) if n >=0]
nx.draw_networkx(G, pos=nx.kamada_kawai_layout(G), arrows=True, nodelist=selected_nodes)

nlist=list(G.nodes.data('LabelGraphics'))
nlist.sort()
for no in nlist:
    if(no[0]>=0):
        elem=str(no[1]).replace("{'text':","").replace(", 'fontSize': 14}","")
        #second=elem.split(":").replace("'","")
        print(no[0]," ",elem)<jupyter_output><empty_output><jupyter_text>Detailed Topology via dot or gml format filesFor a more detailed look we can open up the online version of yEd https://www.yworks.com/yed-live/  and load in the *.dot file that we will now create or dowload the yEd program and look at the *.gml files in the output directory from networkx.drawing.nx_pydot import write_dotimport filer,gfilerfilepath=filer.gui_fname().decode('UTF-8')write_dot(G,filepath) http://localhost:8888/notebooks/Dropbox/1_Jupyter_notebooks/map2loop/notebooks/2.%20map2loop_after_topology.ipynb  Process topography, stratigraphy, fold axial traces and faults Takes GML file produced by topology code, combines with geology polygons, structure points and dtm to create 3D model in gempy.Limitations:  no dykes, no sills. Sills require us to assign a unique surface to each instance of a sill (sill between units A and B needs to be different from sill of same age and strat codes as one found between E and F). Dykes via cokriging are really hard without just cookie cutting them in (but that is not our problem!). We are not checking for onlap relationships, which can perhaps been seen by having lots of units from one series adjacent to the youngest surface of the older series. Could also think about interpreting these as faults to introduce conceptual uncertainty. All mistakes belong to <NAME>, topology code that feeds this system by <NAME>.Geology layer needs to have some unique strat code or text, some group code or text to functionStructure layer needs dip/dip directionCurrently mostly hardwired to GSWA 500K map so needs work...<jupyter_code>os.chdir('../map2loop')

print(os.getcwd())
#from IPython.core.display import display, HTML
#display(HTML("<style>.container { width:80% !important; }</style>"))

sys.path.insert(0,"../..")
print(os.getcwd())


#print(os.getcwd())
#os.environ["PROJ_LIB"] = r"C:\Users\00073294\AppData\Local\Continuum\anaconda3\Lib\site-packages\pyproj\proj_dir\share\proj"<jupyter_output><empty_output><jupyter_text>Next we define an area of interest and some other basic stuff<jupyter_code>print(os.getcwd())
#test_data_path='../test_data3/'


#%run -i "../test_data3/m2l_config.py"

bbox2=str(minx)+","+str(miny)+","+str(maxx)+","+str(maxy)
lat_point_list = [miny, miny, maxy, maxy, maxy]
lon_point_list = [minx, maxx, maxx, minx, minx]
bbox_geom = Polygon(zip(lon_point_list, lat_point_list))
polygon = gpd.GeoDataFrame(index=[0], crs=dst_crs, geometry=[bbox_geom]) 
bbox=(minx,miny,maxx,maxy)

step_out=0.045 #add (in degrees) so edge pixel from dtm reprojection are not found


#contact_decimate=5<jupyter_output><empty_output><jupyter_text>Download and reproject the appropriate SRTM datamj: Getting this from GA, but could also get from Hawaii<jupyter_code>polygon_ll=polygon.to_crs(src_crs)

minlong=polygon_ll.total_bounds[0]-step_out
maxlong=polygon_ll.total_bounds[2]+step_out
minlat=polygon_ll.total_bounds[1]-step_out
maxlat=polygon_ll.total_bounds[3]+step_out
local_paths=False
print(minlong,maxlong,minlat,maxlat)
if((not local_paths) and net):     
    m2l_utils.get_dtm(dtm_file, minlong,maxlong,minlat,maxlat)
    geom_rp=m2l_utils.reproject_dtm(dtm_file,dtm_reproj_file,src_crs,dst_crs)<jupyter_output><empty_output><jupyter_text>Load stratigraphy graph and create list of series (aka groups)mj: The choice of what constitutes basic unit and what a group of units is hard-wired at the moment, but could be altered to any pair. Not even sure we need two levels but it seemed like a good idea at the time. Note that this needs the arcgis plugin version of the topology code (for now) as it seperates the different sub graphs. Text outputs list alternate topologies for series and surfaces, which if confirmed by comapring max-min ages will be a nice source of uncertainty.<jupyter_code>groups,glabels,G = m2l_topology.get_series(strat_graph_file,'id')
m2l_topology.save_units(G,tmp_path,glabels)<jupyter_output><empty_output><jupyter_text>Load geology & structure dataCurrently loading from local files, but could load geology from WFS server at GSWA EXCEPT that the WFS online map has less fields that the zipped shapefiles. Go figure. We don't use fault layer at the moment (except for Vitaliy's topology code) but same logic applies in terms of where to get it from. Already have fault/strat relationships and once we have fault/fault relationships will start to include faults in models.<jupyter_code># Extract point data from structure & geology layers for modelling
##First we readin the structure and map from shapefiles, or wherever...

bbox=(minx,miny,maxx,maxy)
geology = gpd.read_file(geology_file,bbox=bbox)


structure = gpd.read_file(structure_file,bbox=bbox)
structure.crs=dst_crs
print(fault_file)
faults = gpd.read_file(fault_file,bbox=bbox)
faults.crs=dst_crs
#display(faults)
sub_pts = structure[['geometry',c_l['d'],c_l['dd'],c_l['sf']]] 

base=geology.plot(column=c_l['c'],figsize=(10,10),edgecolor='#000000',linewidth=0.2)
sub_pts.plot(ax=base,edgecolor='black')
faults.plot(ax=base, column=c_l['f'],edgecolor='black')<jupyter_output><empty_output><jupyter_text>Clip geology, faults, structures and map geology to structures<jupyter_code>geology = m2l_utils.explode(geology)
geology.crs = dst_crs
pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows',None)
structure_code = gpd.sjoin(sub_pts, geology, how="left", op="within")

y_point_list = [miny, miny, maxy, maxy, miny]
x_point_list = [minx, maxx, maxx, minx, minx]

bbox_geom = Polygon(zip(x_point_list, y_point_list))

polygo = gpd.GeoDataFrame(index=[0], crs=dst_crs, geometry=[bbox_geom]) 
#display(polygo.geometry)
is_bed=structure_code[c_l['sf']].str.contains(c_l['bedding'], regex=False) 
    
all_beds = structure_code[is_bed]
#display(sfcode)
all_folds=faults[faults[c_l['f']].str.contains(c_l['fold'])]
folds_clip=m2l_utils.clip_shp(all_folds,polygo)
folds_clip.to_file(tmp_path+'folds_clip.shp')


all_faults=faults[faults[c_l['f']].str.contains(c_l['fault'])]



#display(structure_code)
geol_clip=m2l_utils.clip_shp(geology, polygo)
faults_clip=m2l_utils.clip_shp(all_faults,polygo)
#display(faults_clip)
structure_clip = m2l_utils.clip_shp(all_beds, polygo)

base = geol_clip.plot(column=c_l['c'],figsize=(7,7),edgecolor='#000000',linewidth=0.2)
faults_clip.plot(ax=base, column=c_l['t'],edgecolor='black')
#structure_clip.plot(ax=base, column=c_l['c'],edgecolor='black')
structure_clip.plot(ax=base, edgecolor='black')


if(c_l['dd']=='strike'):
    structure_clip['azimuth2'] = structure_clip.apply(lambda row: row[c_l['dd']]+90.0, axis = 1)
    c_l['dd']='azimuth2'
    
    
geol_clip.to_file(tmp_path+'geol_clip.shp')
faults_clip.to_file(tmp_path+'faults_clip.shp')
structure_clip.to_file(tmp_path+'structure_clip.shp')<jupyter_output><empty_output><jupyter_text>Create possible stratigraphy sets per groupmj: Uses first of each possible set of toplogies per unit and per group, which is arbitrary. On the other hand we are not checking relative ages again to see if this helps reduce ambiguity, which I think it would.<jupyter_code>m2l_topology.save_group(G,tmp_path,glabels,geol_clip,c_l)<jupyter_output><empty_output><jupyter_text>Export orientation data in csv  formatmj: Orientation data needs calculated height as file does not provide it, taken from SRTM data already downloaded. To calculate polarity (WHICH WE DON'T DO YET) we can calculate the dot product of the dip direction of a bedding plane and the vector to that points nearest basal contact node, if  abs(acos(dot product))>90  then right way up :Added code to not save intrusion orientation data as they won't have associated surfaces if sill...<jupyter_code>dtm = rasterio.open(dtm_reproj_file)
m2l_geometry.save_orientations(structure_clip,output_path,c_l,orientation_decimate,dtm)<jupyter_output><empty_output><jupyter_text>Find those series that don't have any orientation or contact point data  then create arbitrary point for series with no orientation dataNot sure if gempy needs this but geomodeller does. Currently just gives a point dipping 45 degrees to North, but could use dip direction normal to basal surface)<jupyter_code>m2l_geometry.create_orientations( tmp_path, output_path, dtm,geol_clip,structure_clip,c_l)<jupyter_output><empty_output><jupyter_text>Export contact information subset of each polygon to gempy formatmj: Orientation data needs calculated height as file does not provide it, taken from SRTM data already downloaded. Need to reduce number of points whilst retaining useful info (Ranee's job!)'To calculate which are the basal units contact for a polygon find the polygons which are older than the selected polygon, in the example below the central polygon has relative age 23 so its basal contact is with the polygons whose ages are 26 & 28. If there are no older units for a polygon it has no basal content. We keep every nth node based on the decimate term (simple count along polyline). gempy seems to need at least two points per surface, so we always take the first two points.<jupyter_code>ls_dict,ls_dict_decimate=m2l_geometry.save_basal_contacts(tmp_path,dtm,geol_clip,contact_decimate,c_l,intrusion_mode)<jupyter_output><empty_output><jupyter_text>Remove all basal contacts that are defined by faults and save to shapefile (no decimation)<jupyter_code>display(ls_dict)
m2l_geometry.save_basal_no_faults(tmp_path+'basal_contacts.shp',tmp_path+'faults_clip.shp',ls_dict,10,c_l,dst_crs)<jupyter_output><empty_output><jupyter_text>Remove faults from decimated basal contacts as save as csv file<jupyter_code>contacts=gpd.read_file(tmp_path+'basal_contacts.shp')

m2l_geometry.save_basal_contacts_csv(contacts,output_path,dtm,contact_decimate,c_l)

#m2l_geometry.save_contacts_with_faults_removed(tmp_path+'faults_clip.shp',output_path,10,ls_dict,ls_dict_decimate,c_l,dst_crs,dtm)<jupyter_output><empty_output><jupyter_text>Process fault geometrySave Faults as decimated points and representative orientation  Then, for each  fault string:- incementally advance along polyline every at each inter-node (no point in doing more?)- find local stratigraphy 10m to left and right of fault  Once full fault has been traversed:- Find list of contacts left - Find equivalent contacts on right- use interpolated orientations to estimate minimum true offset assuming vertical displacement and store - if no equivalent found, flag as domain fault and find min strat offset for contact, use cumulative minimum thickness estimate and store with flag (not implemented)- estimate median & sd of minimum fault offset and store with flag (not implemented)Local OrientationsSince much of the code is the same, we benefit by calculating local orientation data either side of fault so that geomodeller/gempy have satisfied fault compartment orientation data Save fault as contact info and and orientation info make vertical (for the moment)<jupyter_code>m2l_geometry.save_faults(tmp_path+'faults_clip.shp',output_path,dtm,c_l,fault_decimate,min_fault_length,fault_dip)

use_gcode=('Hamersley Group','Fortescue Group','Wyloo Group','Shingle Creek Group','Turee Creek Group' ) 
use_gcode2=('Hamersley_Group','Fortescue_Group','Wyloo_Group','Shingle_Creek_Group','Turee_Creek_Group' ) 

m2l_interpolation.process_fault_throw_and_near_orientations(tmp_path,output_path,dtm_reproj_file,c_l,use_gcode,use_gcode2,dst_crs,bbox,scheme)<jupyter_output><empty_output><jupyter_text>Process plutonsFor each instruve but not sill polygon, find older neighbours and store decimated contact points. Also store dipping contact orientations (user defined, just because) with four possible sub-surface configurations:saucers: \\_+++_/ batholiths: +++/__ __ _\\+++   domes: /‾+++‾\\ pendants: +++\\_  _/+++   Saves out orientations and contact points, as well as updated group level stratigraphic column.<jupyter_code>bbox=(minx,miny,maxx,maxy)

pluton_dip=str(pluton_dip)

dist_buffer=10
dtm = rasterio.open(dtm_reproj_file)
m2l_geometry.process_plutons(tmp_path,output_path,geol_clip,local_paths,dtm,pluton_form,pluton_dip,contact_decimate,c_l)<jupyter_output><empty_output><jupyter_text>Extract faults and basal contacts of groups from seismic section<jupyter_code>seismic_line_file=data_path+'seismic_line_10GA-CP1_rev.shp'   #input geology file (if local)
seismic_line = gpd.read_file(seismic_line_file) #import map
seismic_line.plot(figsize=(10,10),edgecolor='#000000',linewidth=0.2) #display map
display(seismic_line)


seismic_bbox_file=data_path+'seismic_bbox.shp'   #input geology file (if local)
seismic_bbox = gpd.read_file(seismic_bbox_file) #import map
seismic_bbox.set_index('POSITION',inplace=True)

seismic_interp_file=data_path+'seismic_interp.shp'   #input geology file (if local)
seismic_interp = gpd.read_file(seismic_interp_file) #import map
seismic_interp.plot(column='FEATURE',figsize=(10,10),edgecolor='#000000',linewidth=0.5) #display map
display(seismic_interp)

dtm_reproj_file=dtm_path+'dtm_rp.tif'
dtm = rasterio.open(dtm_reproj_file)

surface_cut=2000

m2l_geometry.extract_section(tmp_path,output_path,seismic_line,seismic_bbox,seismic_interp,dtm,surface_cut)

contacts=pd.read_csv(output_path+'contacts4.csv',",")
seismic_contacts=pd.read_csv(output_path+'seismic_base.csv',",")
all_contacts=pd.concat([contacts,seismic_contacts],sort=False)
all_contacts.to_csv (output_path+'contacts4.csv', index = None, header=True)

faults=pd.read_csv(output_path+'faults.csv',",")
seismic_faults=pd.read_csv(output_path+'seismic_faults.csv',",")
all_faults=pd.concat([faults,seismic_faults],sort=False)
all_faults.to_csv (output_path+'faults.csv', index = None, header=True)<jupyter_output><empty_output><jupyter_text>Interpolates a regular grid of orientations from an  shapefile of arbitrarily-located points and saves out four csv files of l,m & n direction cosines and dip dip direction dataCan choose between various RBF and IDW options     The purpose of these interpolations and associated code is to help in three cases:- Providing estimated dips and contacts in fault-bounded domains where no structural data are available- Needed to estimate true thickness of formations- Possibly useful for poulating parts of maps where little structural data is available<jupyter_code>use_gcode=('Hamersley Group','Fortescue Group','Wyloo Group','Shingle Creek Group','Turee Creek Group' ) ################# MOVE UP   #########################
structure_file=tmp_path+'structure_clip.shp'
bbox=(minx+inset,miny+inset,maxx-inset,maxy-inset)
print(c_l['g'])
print(use_gcode)
m2l_interpolation.interpolate_orientations(structure_file,tmp_path,bbox,c_l,use_gcode,scheme,gridx,gridy,False)<jupyter_output><empty_output><jupyter_text>Interpolates a regular grid from a shapefile of from contacts and saves out three csv files of l & m direction cosines and strike dataCan choose between various RBF and IDW options  Best with basal contacts only, using Notebook 1 as this avoids intrusive  and fault contacts being considered<jupyter_code>geology_file=tmp_path+'basal_contacts.shp'

dtm = rasterio.open(dtm_reproj_file)
use_gcode=('Hamersley_Group','Fortescue_Group','Wyloo_Group','Shingle_Creek_Group','Turee_Creek_Group' ) ################# MOVE UP   #########################

contact_decimate=2
null_scheme='null'

m2l_interpolation.save_contact_vectors(geology_file,tmp_path,dtm,bbox,c_l,null_scheme,contact_decimate)

m2l_interpolation.interpolate_contacts(geology_file,tmp_path,dtm,bbox,c_l,use_gcode,scheme,gridx,gridy,False)<jupyter_output><empty_output><jupyter_text>Combines rbf of contact orientations with rbf of dip/dipdir measurementsUses normalised direction cosines **(l,m,n)**:  - **l,m** from RBF of basal contact orientations  - **signs of l & m** from misorientation with RBF of orientation data and  - **n** from RBF of orientation data    May be useful for adding data where no orientations are available (e.g. in fault bounded domains) and for calculating true thickness of layers. Assumes a 2D plane of data, but if 3D RBF was calulated and projected contact info was used it should apply with topography too.<jupyter_code>geology_file=tmp_path+'geol_clip.shp'
combo_file=tmp_path+'combo.csv'

lc=np.loadtxt(tmp_path+'interpolation_contacts_l.csv',skiprows =1,delimiter =',',dtype=float)
mc=np.loadtxt(tmp_path+'interpolation_contacts_m.csv',skiprows =1,delimiter =',',dtype=float)
lo=np.loadtxt(tmp_path+'interpolation_l.csv',skiprows =1,delimiter =',',dtype=float)
mo=np.loadtxt(tmp_path+'interpolation_m.csv',skiprows =1,delimiter =',',dtype=float)
no=np.loadtxt(tmp_path+'interpolation_n.csv',skiprows =1,delimiter =',',dtype=float)
xy=np.loadtxt(tmp_path+'interpolation_'+scheme+'.csv',skiprows =1,delimiter =',',dtype=float)
m2l_interpolation.join_contacts_and_orientations(combo_file,geology_file,tmp_path,dtm_reproj_file,c_l,lo,mo,no,lc,mc,xy,dst_crs,bbox,False)<jupyter_output><empty_output><jupyter_text>Estimate formation thickness and normalised formation thickness<jupyter_code>buffer =5000
max_thickness_allowed=10000

m2l_geometry.calc_thickness(tmp_path,output_path,buffer,max_thickness_allowed,c_l)

m2l_geometry.normalise_thickness(output_path)<jupyter_output><empty_output><jupyter_text>Creates fold axial trace points<jupyter_code>m2l_geometry.save_fold_axial_traces(tmp_path+'folds_clip.shp',output_path,dtm,c_l,fold_decimate)

#Save fold axial trace near-hinge orientations
fat_step=750         # how much to step out normal to fold axial trace
close_dip=20.0       #dip to assign to all new orientations
   
m2l_geometry.save_fold_axial_traces_orientations(tmp_path+'folds_clip.shp',output_path,tmp_path,dtm,c_l,dst_crs,fold_decimate,fat_step,close_dip)<jupyter_output><empty_output><jupyter_text>Preprocess data to ensure it meets modelling requirements<jupyter_code>if(local_paths): ###############FUDGE#############
    use_gcode=('Hamersley_Group','Fortescue_Group','Turee_Creek_Group','A_mgn_PRK',  'A_mgn_PMI' ) ################# MOVE UP   #########################
    #use_gcode=('Hamersley_Group','Fortescue_Group','Wyloo_Group','Shingle_Creek_Group','Turee_Creek_Group','A_mgn_PRK',  'A_mgn_PMI' ) ################# MOVE UP   #########################
else:
    use_gcode=('Hamersley_Group','Fortescue_Group','Wyloo_Group','Shingle_Creek_Group','Turee_Creek_Group','A_mgn_PMI',  'A_mgn_PRK' ) ################# MOVE UP   #########################


print('only processing',use_gcode)

#inputs=('invented_orientations','interpolated_orientations','intrusive_orientations','fat_orientations','near_fault_orientations')
#inputs=('invented_orientations','interpolated_orientations','intrusive_orientations','fat_orientations','near_fault_orientations')
inputs=('invented_orientations','intrusive_orientations','fat_orientations','near_fault_orientations')

m2l_geometry.tidy_data(output_path,tmp_path,use_gcode,use_interpolations,use_fat,pluton_form,inputs)<jupyter_output><empty_output><jupyter_text>Analyse fault-fault topology<jupyter_code>m2l_topology.parse_fault_relationships(graph_path,tmp_path,output_path)<jupyter_output><empty_output><jupyter_text>loop2gemodeller test<jupyter_code>save_faults=True
compute_etc=True
m2l_export.loop2geomodeller(test_data_path,tmp_path,output_path,'./dtm/dtm_rp.tif',bbox,save_faults,compute_etc)
t1 = time.time()

os.chdir(test_data_path)
%system geomodellerbatch.exe -batch m2l.taskfile
t2 = time.time()
#%system geomodellerbatch.exe -batch m2l_compute.taskfile
t3 = time.time()
import winsound
duration = 500  # milliseconds
freq = 1100  # Hz
winsound.Beep(freq, duration)
winsound.Beep(freq, duration)
duration=999
winsound.Beep(freq, duration)

print("m2l",(t1-t0)/60.0,"batch process",(t2-t1)/60.0,"batch calculate",(t3-t2)/60.0,"minutes")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11456707930233365, 'fraction_numerical': 0.011928291549637643, 'mean_word_length': 6.849755035383778, 'pattern_counts': {'":': 31, '<': 150, '<?xml version=': 0, '>': 153, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '36', 'dir': 'jupyter-structured-clean-dedup', 'id': '23622', 'license': 'MIT', 'n_tokens_mistral': 10371, 'n_tokens_neox': 9903, 'n_words': 2314, 'path': 'notebooks/103 Mineral System Analysis.ipynb', 'repo_name': 'bluetyson/map2loop', 'size': '28861'}