Example 2 - Draw Your Own Area.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Loop Workflow Example 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> * High level approach to making a 3D model from just a bounding box you can draw\n",
    "> * To run this notebook for the first time, there are some dependencies needed to use the interactive map. Execute the cell immediately below and restart jupyter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#if not already installed:\n",
    "#!conda install -c loop3d map2loop loopstructural pyamg meshio\n",
    "#!pip install loopstructuralvisualisation[jupyter]\n",
    "\n",
    "# note on versions:\n",
    "\n",
    "## map2loop version = v3.2.0\n",
    "## LoopStructural version = 1.6.5\n",
    "## LoopProjectFile = 0.2.2\n",
    "## loopstructural-visualisation = 0.1.10\n",
    "\n",
    "\n",
    "#!conda install -c conda-forge leafmap ipyleaflet ipywidgets -y\n",
    "# and\n",
    "#!pip install geoh5py\n",
    "\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Leaflet Maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "25bb2148ce664c22b0167ce3c12985c1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Dropdown(description='Map area:', options=('Draw Your Own', 'Last Area Drawn'), value='Draw Your Own')"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import ipywidgets as widgets\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# load last saved map area and model engine (if they exist)\n",
    "if(not os.path.isdir('../scratch/')):\n",
    "    os.mkdir('../scratch/')\n",
    "if(os.path.isfile('../scratch/last_choices.txt')):\n",
    "    f=open('../scratch/last_choices.txt','r')\n",
    "    contents =f.readlines()\n",
    "    f.close()\n",
    "    default_map=contents[0].replace(\"\\n\",\"\")\n",
    "    default_engine=contents[1].replace(\"\\n\",\"\")\n",
    "else:\n",
    "    default_map='Turner_Syncline'\n",
    "    default_engine='loopstructural'\n",
    "\n",
    "options=['Draw Your Own','Last Area Drawn']\n",
    "\n",
    "if(not default_map in options):\n",
    "    default_map= options[0]\n",
    "\n",
    "map_choice=widgets.Dropdown(\n",
    "    options=options,\n",
    "    value=default_map,\n",
    "    description='Map area:',\n",
    "    disabled=False,\n",
    ")\n",
    "display(map_choice)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Draw Your Own\n"
     ]
    }
   ],
   "source": [
    "test_data_name=map_choice.value\n",
    "print(test_data_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import leafmap\n",
    "import pandas as pd\n",
    "from shapely.geometry import Polygon\n",
    "import geopandas as gpd\n",
    "src_crs = \"epsg:4326\"  # coordinate reference system for imported dtms (geodetic lat/long WGS84)\n",
    "dst_crs = \"epsg:28350\"  # coordinate reference system for imported dtms (geodetic lat/long WGS84)\n",
    "\n",
    "if(not test_data_name =='Draw Your Own'):\n",
    "    if(test_data_name=='Last Area Drawn'):\n",
    "        last_coords=pd.read_csv('../scratch/last_area.csv')\n",
    "        display(last_coords)\n",
    "        minx=last_coords.iloc[0]['minx']\n",
    "        miny=last_coords.iloc[0]['miny']\n",
    "        maxx=last_coords.iloc[0]['maxx']\n",
    "        maxy=last_coords.iloc[0]['maxy']\n",
    "    elif(not test_data_name =='Draw Your Own'):\n",
    "        y_point_list = [miny, miny, maxy, maxy, maxy]\n",
    "        x_point_list = [minx, maxx, maxx, minx, minx]\n",
    "        bbox_geom = Polygon(zip(x_point_list, y_point_list))\n",
    "        polygon = gpd.GeoDataFrame(index=[0], crs=dst_crs, geometry=[bbox_geom])\n",
    "        polygon_ll=polygon.to_crs(src_crs)\n",
    "\n",
    "        minx=polygon_ll.total_bounds[0]\n",
    "        maxx=polygon_ll.total_bounds[2]\n",
    "        miny=polygon_ll.total_bounds[1]\n",
    "        maxy=polygon_ll.total_bounds[3]\n",
    "\n",
    "        minlong=minx\n",
    "        maxlong=maxx\n",
    "        minlat=miny\n",
    "        maxlat=maxy\n",
    "        #print(\"x\",polygon_ll.total_bounds[0])\n",
    "        st_bbox=[minlong,minlat,maxlong,maxlat]\n",
    "        lat_point_list = [minlat, minlat, maxlat, maxlat,maxlat]\n",
    "        lon_point_list = [minlong, maxlong, maxlong, minlong, minlong]\n",
    "        bbox_geom = Polygon(zip(lon_point_list, lat_point_list))\n",
    "        rect = gpd.GeoDataFrame(index=[0], crs=src_crs, geometry=[bbox_geom]) \n",
    "\n",
    "    bbox2=str(minx)+\",\"+str(miny)+\",\"+str(maxx)+\",\"+str(maxy)\n",
    "    y_point_list = [miny, miny, maxy, maxy, maxy]\n",
    "    x_point_list = [minx, maxx, maxx, minx, minx]\n",
    "    bbox_geom = Polygon(zip(x_point_list, y_point_list))\n",
    "    polygon = gpd.GeoDataFrame(index=[0], crs=dst_crs, geometry=[bbox_geom])\n",
    "    polygon_ll=polygon.to_crs(src_crs)\n",
    "\n",
    "    minlong=polygon_ll.total_bounds[0]\n",
    "    maxlong=polygon_ll.total_bounds[2]\n",
    "    minlat=polygon_ll.total_bounds[1]\n",
    "    maxlat=polygon_ll.total_bounds[3]\n",
    "    \n",
    "    minlong=minx\n",
    "    maxlong=maxx\n",
    "    minlat=miny\n",
    "    maxlat=maxy\n",
    "\n",
    "    lat_point_list = [minlat, minlat, maxlat, maxlat,maxlat]\n",
    "    lon_point_list = [minlong, maxlong, maxlong, minlong, minlong]\n",
    "    bbox_geom = Polygon(zip(lon_point_list, lat_point_list))\n",
    "    rect = gpd.GeoDataFrame(index=[0], crs=src_crs, geometry=[bbox_geom]) \n",
    "\n",
    "    center=(minlat+((maxlat-minlat)/2),minlong+((maxlong-minlong)/2))\n",
    "else:\n",
    "    center=(-22.6,117.3)\n",
    "\n",
    "m = leafmap.Map( center=center, zoom=8,scroll_wheel_zoom=True)\n",
    "m.add_basemap( basemap='OpenTopoMap')\n",
    "m.add_wms_layer(url='https://www.loopwms.xyz/geoserver/loop/wms?',\n",
    "    layers='2_5m_interpgeop15_4326',\n",
    "    format='image/png',\n",
    "    transparent=True,\n",
    "    opacity=0.4,\n",
    "    attribution='Geology data from GSWA',\n",
    "    name='geology',\n",
    "    shown=True)\n",
    "# url = 'https://www.loopwms.xyz/geoserver/loop/wms?'\n",
    "m.add_wms_layer(url='https://www.loopwms.xyz/geoserver/loop/wms?',\n",
    "    layers='2_5m_interpstrucl15_4326',\n",
    "    format='image/png',\n",
    "    transparent=True,\n",
    "    opacity=0.4,\n",
    "    attribution='Linear features data from GSWA',\n",
    "    name='faults/folds')\n",
    "m.add_wms_layer(url='https://www.loopwms.xyz/geoserver/loop/wms?',\n",
    "    layers='waroxi_wa_4326_bed',\n",
    "    format='image/png',\n",
    "    transparent=True,\n",
    "    attribution='Outcrop data from GSWA',\n",
    "    name='outcrops')\n",
    "\n",
    "if( not test_data_name =='Draw Your Own'):\n",
    "    m.add_gdf(rect, layer_name=\"Last Rectangle\", \n",
    "                   fill_colors=['orange'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epsg:4326 [117.13698  -22.690712 117.564464 -22.396454]\n",
      "epsg:28350 [ 514069.70809875 7490607.48994945  558101.7853635  7523283.65581696]\n"
     ]
    }
   ],
   "source": [
    "polyclip=False # implemented in future version\n",
    "\n",
    "if(test_data_name=='Draw Your Own' or test_data_name=='Last Area Drawn'):\n",
    "    if(test_data_name=='Draw Your Own'):\n",
    "        src_crs = \"epsg:4326\"  # coordinate reference system for imported dtms (geodetic lat/long WGS84)\n",
    "        dst_crs = \"epsg:28350\" # coordinate system for example data\n",
    "        \n",
    "        rect=m.draw_features\n",
    "\n",
    "        if(len(rect)==0):\n",
    "            minlong=117.13698\n",
    "            maxlong=117.564464\n",
    "            minlat= -22.690712\n",
    "            maxlat=-22.396454\n",
    "\n",
    "        else:\n",
    "            minlong=rect[0]['geometry']['coordinates'][0][0][0]\n",
    "            maxlong=rect[0]['geometry']['coordinates'][0][2][0]\n",
    "            minlat= rect[0]['geometry']['coordinates'][0][0][1]\n",
    "            maxlat= rect[0]['geometry']['coordinates'][0][1][1]    \n",
    "    else:\n",
    "        use_roi_clip=False\n",
    "        roi_clip_path=''\n",
    "\n",
    "    bounds=(minlong,maxlong,minlat,maxlat)\n",
    "\n",
    "    lat_point_list = [minlat, minlat, maxlat, maxlat,maxlat]\n",
    "    lon_point_list = [minlong, maxlong, maxlong, minlong, minlong]\n",
    "    bbox_geom = Polygon(zip(lon_point_list, lat_point_list))\n",
    "    mbbox = gpd.GeoDataFrame(index=[0], crs=src_crs, geometry=[bbox_geom]) \n",
    "    bbox=mbbox.total_bounds\n",
    "    st_bbox=[bbox[0],bbox[1],bbox[2],bbox[3]]\n",
    "    print(src_crs,mbbox.total_bounds)\n",
    "    mbbox=mbbox.to_crs(dst_crs)\n",
    "    print(dst_crs,mbbox.total_bounds)\n",
    "    \n",
    "    f=open('../scratch/last_area.csv','w') \n",
    "    ostr='minx,miny,maxx,maxy\\n'\n",
    "    f.write(ostr)\n",
    "    ostr=str(minlong)+','+str(minlat)+','+str(maxlong)+','+str(maxlat)+'\\n'\n",
    "    f.write(ostr)\n",
    "    f.close()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Map2Loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: 2024-12-17 16:29:37,009: mapdata.py:466 -- Datatype FAULT_ORIENTATION is not set and so cannot be loaded\n",
      "\n",
      "WARNING: 2024-12-17 16:29:51,581: interpolators.py:352 -- Detected 12 collocated point clusters. Aggregating these points.\n",
      " \n",
      "WARNING: 2024-12-17 16:30:27,723: thickness_calculator.py:407 -- Thickness Calculator InterpolatedStructure: Cannot calculate thickness between Bunjinah_Formation and Pyradie_Formation\n",
      "\n",
      "WARNING: 2024-12-17 16:30:27,730: thickness_calculator.py:407 -- Thickness Calculator InterpolatedStructure: Cannot calculate thickness between Pyradie_Formation and Fortescue_Group\n",
      "\n",
      "WARNING: 2024-12-17 16:30:30,006: project.py:659 -- No fault orientation data found, skipping fault orientation calculation\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "map2loop 58.82671928405762 seconds\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "from map2loop.project import Project\n",
    "from map2loop.m2l_enums import VerboseLevel\n",
    "from datetime import datetime\n",
    "from map2loop.thickness_calculator import InterpolatedStructure\n",
    "import time\n",
    "import warnings\n",
    "\n",
    "\n",
    "t0 = time.time()\n",
    "warnings.filterwarnings('ignore')\n",
    "nowtime=datetime.now().isoformat(timespec='minutes')\n",
    "bbox_3d={\n",
    "    \"minx\": mbbox.total_bounds[0], #500000,\n",
    "    \"miny\": mbbox.total_bounds[1], #7490000,\n",
    "    \"maxx\": mbbox.total_bounds[2], #545000,\n",
    "    \"maxy\": mbbox.total_bounds[3], #7520000,\n",
    "    \"base\": -4800,\n",
    "    \"top\": 1200,\n",
    "}\n",
    "model_name=os.path.join(\".\",nowtime.replace(\":\",\"-\").replace(\"T\",\"-\"))\n",
    "loop_project_filename=os.path.join(\"wa_output.loop3d\")\n",
    "proj = Project(\n",
    "    use_australian_state_data = \"WA\",\n",
    "    verbose_level = VerboseLevel.NONE,\n",
    "    tmp_path  = model_name,\n",
    "    working_projection = \"EPSG:28350\",\n",
    "    bounding_box = bbox_3d,\n",
    "    loop_project_filename = loop_project_filename, \n",
    "    overwrite_loopprojectfile = True\n",
    ")\n",
    "proj.set_thickness_calculator(InterpolatedStructure())\n",
    "proj.run_all(take_best=True)\n",
    "\n",
    "t1 = time.time()\n",
    "m2l_time = t1-t0\n",
    "m2l_string = f\"{m2l_time} seconds\" if m2l_time < 60 else f\"{m2l_time/60.0} minutes\"\n",
    "print(f\"map2loop {m2l_string}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loop Structural"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ERROR: 2024-12-17 16:30:46,834: process_data.py:460 -- Thickness for Pyradie_Formation is less than or equal to 0\n",
      " Update the thickness value for Pyradie_Formation before continuing\n"
     ]
    }
   ],
   "source": [
    "import LoopProjectFile as LPF\n",
    "import LoopStructural\n",
    "from LoopStructural.modelling.input.project_file import LoopProjectfileProcessor as LPFProcessor\n",
    "from LoopStructural.visualisation import Loop3DView\n",
    "import numpy as np\n",
    "\n",
    "# LoopStructural Modelling\n",
    "fault_params = {\n",
    "    'interpolatortype':'FDI',\n",
    "    'nelements':1e4,\n",
    "}\n",
    "foliation_params = {\n",
    "    'interpolatortype':'FDI' , # 'interpolatortype':'PLI',\n",
    "    'nelements':1e5,  # how many tetras/voxels\n",
    "    'regularisation':5,\n",
    "}\n",
    "projFile = LPF.ProjectFile(loop_project_filename)\n",
    "processedData = LPFProcessor(projFile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "processedData.thicknesses['Pyradie_Formation']=600"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: 2024-12-17 16:31:37,117: _fault_builder.py:153 -- No orientation data for fault\n",
      "                    Defaulting to a dip of 90vertical fault\n",
      "WARNING: 2024-12-17 16:31:37,132: _fault_builder.py:208 -- There is no slip vector data for the fault, using vertical slip vector\n",
      "                          projected onto fault surface estimating from fault normal\n",
      "WARNING: 2024-12-17 16:31:37,245: _fault_builder.py:153 -- No orientation data for fault\n",
      "                    Defaulting to a dip of 90vertical fault\n",
      "WARNING: 2024-12-17 16:31:37,253: _fault_builder.py:208 -- There is no slip vector data for the fault, using vertical slip vector\n",
      "                          projected onto fault surface estimating from fault normal\n",
      "WARNING: 2024-12-17 16:31:37,335: _fault_builder.py:153 -- No orientation data for fault\n",
      "                    Defaulting to a dip of 90vertical fault\n",
      "WARNING: 2024-12-17 16:31:37,349: _fault_builder.py:208 -- There is no slip vector data for the fault, using vertical slip vector\n",
      "                          projected onto fault surface estimating from fault normal\n",
      "WARNING: 2024-12-17 16:31:37,419: _fault_builder.py:153 -- No orientation data for fault\n",
      "                    Defaulting to a dip of 90vertical fault\n",
      "WARNING: 2024-12-17 16:31:37,445: _fault_builder.py:208 -- There is no slip vector data for the fault, using vertical slip vector\n",
      "                          projected onto fault surface estimating from fault normal\n",
      "WARNING: 2024-12-17 16:31:37,534: _fault_builder.py:153 -- No orientation data for fault\n",
      "                    Defaulting to a dip of 90vertical fault\n",
      "WARNING: 2024-12-17 16:31:37,539: _fault_builder.py:208 -- There is no slip vector data for the fault, using vertical slip vector\n",
      "                          projected onto fault surface estimating from fault normal\n",
      "WARNING: 2024-12-17 16:31:37,612: _fault_builder.py:153 -- No orientation data for fault\n",
      "                    Defaulting to a dip of 90vertical fault\n",
      "WARNING: 2024-12-17 16:31:37,622: _fault_builder.py:208 -- There is no slip vector data for the fault, using vertical slip vector\n",
      "                          projected onto fault surface estimating from fault normal\n",
      "WARNING: 2024-12-17 16:31:37,703: _fault_builder.py:153 -- No orientation data for fault\n",
      "                    Defaulting to a dip of 90vertical fault\n",
      "WARNING: 2024-12-17 16:31:37,712: _fault_builder.py:208 -- There is no slip vector data for the fault, using vertical slip vector\n",
      "                          projected onto fault surface estimating from fault normal\n",
      "WARNING: 2024-12-17 16:31:37,786: _fault_builder.py:153 -- No orientation data for fault\n",
      "                    Defaulting to a dip of 90vertical fault\n",
      "WARNING: 2024-12-17 16:31:37,795: _fault_builder.py:208 -- There is no slip vector data for the fault, using vertical slip vector\n",
      "                          projected onto fault surface estimating from fault normal\n",
      "WARNING: 2024-12-17 16:31:37,893: _fault_builder.py:153 -- No orientation data for fault\n",
      "                    Defaulting to a dip of 90vertical fault\n",
      "WARNING: 2024-12-17 16:31:37,893: _fault_builder.py:208 -- There is no slip vector data for the fault, using vertical slip vector\n",
      "                          projected onto fault surface estimating from fault normal\n",
      "WARNING: 2024-12-17 16:31:37,985: _fault_builder.py:153 -- No orientation data for fault\n",
      "                    Defaulting to a dip of 90vertical fault\n",
      "WARNING: 2024-12-17 16:31:38,000: _fault_builder.py:208 -- There is no slip vector data for the fault, using vertical slip vector\n",
      "                          projected onto fault surface estimating from fault normal\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e8f037beb7924747bfade410647b3d6d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/31 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "processedData.foliation_properties['sg'] = foliation_params\n",
    "processedData.fault_properties['interpolatortype'] = fault_params['interpolatortype']\n",
    "processedData.fault_properties['nelements'] = fault_params['nelements']\n",
    "\n",
    "model = LoopStructural.GeologicalModel.from_processor(processedData)\n",
    "model.nsteps=np.array([100,100,50])\n",
    "model.update()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "61c50ffe191040debf7cf6622b208391",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Widget(value='<iframe src=\"http://localhost:52981/index.html?ui=P_0x1b99d7dd690_0&reconnect=auto\" class=\"pyvis…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view = Loop3DView(model)\n",
    "view.plot_model_surfaces()\n",
    "view.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}