3b. South Australia Model with m2l.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Case study 3: Flinders ranges\n",
    "This example should be run using the Loop3d/Loop docker container. It requires \n",
    "* LoopStructural\n",
    "* map2loop\n",
    "* folium\n",
    "* ipyleaflet\n",
    "* ipywidgets\n",
    "* lavavu-osmesa\n",
    "\n",
    "It cannot be easily run on google colab, if you are using google colab run [this](South%20Australia%20Model%20No%20map2loop.ipynb) notebook instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import folium\n",
    "from ipyleaflet import Map, basemaps, GeoJSON, LayersControl, DrawControl,WMSLayer, GeoData\n",
    "from ipywidgets import Label\n",
    "from ipywidgets import Label\n",
    "import ipywidgets as widgets\n",
    "import numpy as np\n",
    "from shapely.geometry import Polygon\n",
    "import geopandas\n",
    "import pandas as pd\n",
    "import lavavu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_crs = 'EPSG:4326'\n",
    "proj_crs = 'EPSG:28354'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we define a bounding box for the model, it will be displayed on a map in the next cell. You can use the map to draw your own bounding box if you want to try different model areas. \n",
    "\n",
    "**Warning** if you try new areas it may not work as expected because the modelling process assumes chronostratigraphic stratigraphy and requires adequate structural data. The structural data are shown by the red dots. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "existing_extent = (138.362206, -32.463842, 139.262268, -31.999178)\n",
    "center = (existing_extent[1]+existing_extent[3])/2, (existing_extent[0]+existing_extent[2])/2\n",
    "minlat = existing_extent[3]\n",
    "maxlat = existing_extent[1]\n",
    "minlong = existing_extent[2]\n",
    "maxlong = existing_extent[0]\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 = geopandas.GeoDataFrame(index=[0], crs=map_crs, geometry=[bbox_geom])\n",
    "example_rect = GeoData(geo_dataframe = mbbox,\n",
    "                   style={'color': 'purple', 'opacity':3, 'weight':1.9, 'dashArray':'2', 'fillOpacity':0.6},                  \n",
    "                   name = 'Example')\n",
    "print(map_crs,mbbox.total_bounds)\n",
    "mbbox=mbbox.to_crs(proj_crs)\n",
    "print(proj_crs,mbbox.total_bounds)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# center=(-33.5,138.7)\n",
    "\n",
    "wms_warox = WMSLayer(\n",
    "    url='http://geo.loop-gis.org/geoserver/GSSA/wms?',\n",
    "    layers='GSSA:sth_flinders',\n",
    "    format='image/png',\n",
    "    transparent=True,\n",
    "    attribution='Outcrop data from GSSA',\n",
    "    name='outcrops'\n",
    "\n",
    ")\n",
    "wms_geol = WMSLayer(\n",
    "    url='http://geo.loop-gis.org/geoserver/GSSA/wms?',\n",
    "    layers='GSSA:2m surface geology',\n",
    "    format='image/png',\n",
    "    transparent=True,\n",
    "    opacity=0.4,\n",
    "    attribution='Geology data from GSSA',\n",
    "    name='geology'\n",
    "\n",
    ")\n",
    "m =Map(basemap=basemaps.OpenTopoMap, center=center, zoom=8,scroll_wheel_zoom=True)\n",
    "label = Label()\n",
    "display(label)\n",
    "\n",
    "def handle_interaction(**kwargs):\n",
    "    if kwargs.get('type') == 'mousemove':\n",
    "        label.value = str(kwargs.get('coordinates'))\n",
    "\n",
    "m.on_interaction(handle_interaction)\n",
    "m.add_layer(wms_geol)\n",
    "m.add_layer(wms_warox)\n",
    "# if(not test_data_name =='Draw Your Own'):\n",
    "#     m.add_layer(example_rect)\n",
    "\n",
    "\n",
    "m.add_control(LayersControl())\n",
    "dc = DrawControl(rectangle={'shapeOptions': {'color': '#0000FF'}})\n",
    "m.add_control(dc)\n",
    "m.add_layer(example_rect)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geopandas import GeoDataFrame\n",
    "from shapely.geometry import shape \n",
    "draw = dc.last_draw\n",
    "try:\n",
    "    poly = shape(draw['geometry'])\n",
    "    poly_gdf = GeoDataFrame([1], geometry=[poly], crs=map_crs)\n",
    "    bbox_gdf = poly_gdf.to_crs(proj_crs)\n",
    "    bbox_3d = {\"minx\": np.min(bbox_gdf.loc[:,'geometry'][0].exterior.xy[0]),\n",
    "         \"miny\": np.min(bbox_gdf.loc[:,'geometry'][0].exterior.xy[1]),\n",
    "         \"maxx\": np.max(bbox_gdf.loc[:,'geometry'][0].exterior.xy[0]),\n",
    "         \"maxy\": np.max(bbox_gdf.loc[:,'geometry'][0].exterior.xy[1]),\n",
    "         \"base\": -3200,\n",
    "         \"top\": 1200,}\n",
    "except:\n",
    "    minx = mbbox.total_bounds[0]\n",
    "    maxx = mbbox.total_bounds[2]\n",
    "    miny = mbbox.total_bounds[1]\n",
    "    maxy = mbbox.total_bounds[3]# [ 274934.13200956 6297758.41393543  323980.90024616 6329307.64682912]\n",
    "    bbox_3d = {'minx': minx,\n",
    "     'miny': miny,\n",
    "     'maxx': maxx,\n",
    "     'maxy': maxy,\n",
    "     'base': -3200,\n",
    "     'top': 1200}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Map2Loop\n",
    "This process uses map2loop to process the goelogical datasets. Here we use local files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import hjson\n",
    "import map2loop\n",
    "map2loop._clut_path='SA/data/GSSA_2M_colours.csv'\n",
    "from map2loop.project import Project\n",
    "proj = Project( \n",
    "                 structure_file='SA/data/sth_flinders_28354.shp',\n",
    "                fault_file='SA/data/2M Linear Structures_28354.shp',\n",
    "                fold_file='SA/data/2M Linear Structures_28354.shp',\n",
    "                geology_file='SA/data/2M_Surface_Geology_28354_relage.shp',\n",
    "                mindep_file=\"SA/data/null_mindeps.shp\",\n",
    "                 metadata='SA/data/meta.hjson',\n",
    "#                  state = \"SA\",\n",
    "#                  remote = False,\n",
    "#                 clut_path='SA/data/GSSA_2M_colours.csv'\n",
    "                 # path to hosted metadata describing the remote sources\n",
    "#                  metadata='https://gist.githubusercontent.com/yohanderose/8f843de0dde531f009a3973cbdadcc9f/raw/918f412ae488ce1a6bca188306f7730061ecf551/meta_remote.hjson')                \n",
    ")\n",
    "proj.update_config(\n",
    "                    out_dir='sa-test2',\n",
    "                    bbox_3d=bbox_3d,\n",
    "                    proj_crs={'init': 'EPSG:28354'},\n",
    "#                     quiet=True,\n",
    "#                     overwrite=\"true\",\n",
    "                    drift_prefix=['T','Q','water','void']\n",
    "                  )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "proj.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loop Structural"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define project pathing from m2l\n",
    "proj_path = proj.config.project_path\n",
    "graph_path = proj.config.graph_path\n",
    "tmp_path = proj.config.tmp_path\n",
    "data_path = proj.config.data_path\n",
    "dtm_path = proj.config.dtm_path\n",
    "output_path = proj.config.output_path\n",
    "vtk_path = proj.config.vtk_path\n",
    "\n",
    "# Define project bounds\n",
    "minx,miny,maxx,maxy = proj.config.bbox\n",
    "model_base = proj.config.bbox_3d['base']\n",
    "model_top = proj.config.bbox_3d['top']\n",
    "\n",
    "fault_file = proj.config.fault_file_csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import random\n",
    "import os\n",
    "import time\n",
    "from datetime import datetime\n",
    "import shutil\n",
    "import logging\n",
    "logging.getLogger().setLevel(logging.INFO)\n",
    "\n",
    "import numpy as np\n",
    "from LoopStructural import GeologicalModel\n",
    "import lavavu\n",
    "from LoopStructural.visualisation import LavaVuModelViewer\n",
    "from LoopStructural import GeologicalModel\n",
    "\n",
    "f=open(os.path.join(tmp_path, 'bbox.csv'),'w')\n",
    "f.write('minx,miny,maxx,maxy,lower,upper\\n')\n",
    "ostr='{},{},{},{},{},{}\\n'.format(minx,miny,maxx,maxy,model_base,model_top)\n",
    "f.write(ostr)\n",
    "f.close()\n",
    "\n",
    "fault_params = {'interpolatortype':'FDI',\n",
    "                'nelements':3e5,\n",
    "                'fault_buffer':.5,\n",
    "                'solver':'pyamg',\n",
    "#                 overprints:overprints,\n",
    "                'cpw':10,\n",
    "                'npw':10}\n",
    "foliation_params = {'interpolatortype':'FDI' , # 'interpolatortype':'PLI', 'FDI', 'surfe'\n",
    "                    'nelements':3e5,  # how many tetras/voxels\n",
    "                    'buffer':1.8,  # how much to extend interpolation around box\n",
    "                    'solver':'pyamg',\n",
    "                    'damp':True,\n",
    "                    'cpw':5,\n",
    "                    'npw':5}\n",
    "\n",
    "model, m2l_data = GeologicalModel.from_map2loop_directory(proj_path,\n",
    "                                                          skip_faults=False,\n",
    "                                                          rescale=False,\n",
    "                                                          fault_params=fault_params,\n",
    "                                                          foliation_params=foliation_params)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view = LavaVuModelViewer(model,vertical_exaggeration=1) \n",
    "view.nsteps = np.array([200,200,100])\n",
    "#view.set_zscale(2)\n",
    "# view.add_model()\n",
    "# view.nelements = 3e5#steps=np.array([,50,50])\n",
    "#view.add_model_surfaces(filename=filename)\n",
    "view.add_model_surfaces(opacity=0.8)\n",
    "view.add_data(model['supergroup_0'])\n",
    "view.interactive()  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# to export the figure to webgl\n",
    "# view.export_to_webgl('south_australia_model.html')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geol_clip = proj.config.geol_clip\n",
    "asc=pd.read_csv(tmp_path+'/all_sorts_clean.csv',\",\")\n",
    "#display(asc)\n",
    "colours=pd.read_csv(proj.config.clut_path,\",\")\n",
    "if( proj.config.c_l['c']=='CODE'):\n",
    "    code=proj.config.c_l['c'].lower()\n",
    "else:\n",
    "    code=proj.config.c_l['c']\n",
    "\n",
    "colours = [] #container for the discrete colours we are using\n",
    "i=0\n",
    "geol_clip['colour_index'] = np.nan #initialise a colour index attribute column\n",
    "for ind,strat in asc.iterrows():\n",
    "    geol_clip[proj.config.c_l['c']].str.replace(\" \",\"_\")\n",
    "    geol_clip.loc[geol_clip[proj.config.c_l['c']]==strat['code'].replace(\"_\",\" \"),'colour_index'] = i\n",
    "    colours.append(strat['colour'])\n",
    "    i=i+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from LoopStructural.visualisation import MapView\n",
    "from LoopStructural.visualisation.stratigraphic_column import StratigraphicColumnView\n",
    "from map2loop.map import MapUtil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maps =proj.config.geol_clip#geopandas.read_file(proj.config.tmp_path+'geol_clip.shp')\n",
    "masp = maps[maps['STRATNAM_1']!='Unnamed GIS Unit - see description']\n",
    "\n",
    "extent=[model.bounding_box[0,0], model.bounding_box[1,0], model.bounding_box[0,1], model.bounding_box[1,1]]\n",
    "\n",
    "\n",
    "mapview = MapView(model)\n",
    "proj.config.load_dtm(True)\n",
    "geomap = MapUtil(proj.config.bbox_3d,maps,proj.config.dtm)\n",
    "mapview.nsteps = (200,200)\n",
    "xy = np.array([mapview.xx\n",
    "               .flatten(),mapview.yy.flatten()]).T\n",
    "dtm = geomap.evaluate_dtm_at_points(xy)\n",
    "\n",
    "geo_map = geomap.evaluate_geology_at_points(xy)\n",
    "# dtm[:]=0\n",
    "model_map = model.evaluate_model(np.vstack([xy.T,dtm]).T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = []\n",
    "colours = []\n",
    "boundaries = []\n",
    "for u, v in m2l_data['stratigraphic_column']['supergroup_0'].items():\n",
    "    data.append((m2l_data['strat_va'][u],v['colour']))\n",
    "    colours.append(v['colour'])\n",
    "    boundaries.append(m2l_data['strat_va'][u])    \n",
    "from matplotlib import colors\n",
    "cmap = colors.ListedColormap(colours)\n",
    "cmap = colors.ListedColormap(cmap.colors[::-1])\n",
    "b = np.array(boundaries)\n",
    "print()\n",
    "b = b[::-1]\n",
    "boundaries = colors.BoundaryNorm(b,ncolors=len(colours))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot geology map and model comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'font.size': 32})\n",
    "fig3 = plt.figure(constrained_layout=True,figsize=(60,15))\n",
    "nx = 6\n",
    "gs = fig3.add_gridspec(5, 2*nx+1)\n",
    "legend = fig3.add_subplot(gs[:-1, :1])\n",
    "strike = fig3.add_subplot(gs[-1:,:1])\n",
    "\n",
    "#add stike symbol\n",
    "gradient_data = np.array([[0,1.,0]])\n",
    "t = gradient_data[:, [1, 0]] * np.array([1, -1]).T\n",
    "n = gradient_data[:, 0:2]\n",
    "# t *= symb_scale\n",
    "n *= 0.5\n",
    "p1 = gradient_data[:, [0, 1]] - t\n",
    "p2 = gradient_data[:, [0, 1]] + t\n",
    "# plt.scatter(val[:,0],val[:,1],c='black')\n",
    "strike.plot([p1[:, 0], p2[:, 0]], [p1[:, 1], p2[:, 1]], 'black')\n",
    "p1 = gradient_data[:, [0, 1]]\n",
    "p2 = gradient_data[:, [0, 1]] + n\n",
    "strike.plot([p1[:, 0], p2[:, 0]], [p1[:, 1], p2[:, 1]], 'black')\n",
    "strike.annotate(\"Bedding orientation\", xy=[0,0],xytext=[1.5,.75],annotation_clip=False,size=32)\n",
    "strike.set_xlim(-2,2)\n",
    "# strike.plot()\n",
    "strike.axis('off')\n",
    "strike.axis('square')\n",
    "p1 = gradient_data[:, [0, 1]] - t\n",
    "p2 = gradient_data[:, [0, 1]] + t\n",
    "strike.plot([p1[:, 0], p2[:, 0]], [p1[:, 1]+1.2, p2[:, 1]+1.2], 'black')\n",
    "strike.annotate(\"Fault trace\", xy=[0,1.2],xytext=[1.5,.75+1.2],annotation_clip=False,size=32)\n",
    "\n",
    "# f3_ax1.set_title('gs[0, :-2]')\n",
    "geol_map_ax = fig3.add_subplot(gs[:, 1:nx+1])\n",
    "model_map_ax = fig3.add_subplot(gs[:, nx+1:])\n",
    "\n",
    "column_ax = StratigraphicColumnView(model,legend)\n",
    "geol_map_ax.imshow(geo_map.reshape(mapview.nsteps).T,extent=extent,cmap=cmap,origin='lower',vmin=0,vmax=7)\n",
    "# proj.config.faults_clip.plot(ax=geol_map_ax)\n",
    "mapview2 = MapView(model,ax=model_map_ax)\n",
    "mapview2.nsteps=(200,200)\n",
    "mapview2.add_data(model.features[-1],dip=False,val=False,symb_scale=600)\n",
    "mapview2.add_model(cmap=cmap,z=dtm)\n",
    "mapview2.add_faults(colors='black')\n",
    "model_map_ax.imshow(model_map.reshape(mapview.nsteps).T,extent=extent,cmap=cmap,origin='lower',vmin=0,vmax=7)\n",
    "\n",
    "legend.set_title('A. Legend',fontsize=40)\n",
    "geol_map_ax.set_title('B. Geological map',fontsize=40)\n",
    "model_map_ax.set_title('C. Interpolated geological map',fontsize=40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}