<jupyter_start><jupyter_text>Interactive exploration of AEM log dataAdapted from [voila-gpx-viewer](https://github.com/jtpio/voila-gpx-viewer), found in the [Voila Gallery](https://voila-gallery.org/services/gallery/)Preparations:```bashconda activate ELAconda install -c conda-forge bqplot ipyleafletjupyter labextension install jupyter-leaflet bqplotjupyter labextension listcd /path/to/parent/of/thisfilejupyter-lab .```<jupyter_code>import os
import sys
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import rasterio
from rasterio.plot import show
import geopandas as gpd

import datetime
import json

from io import StringIO
from statistics import mean

from bqplot import Axis, Figure, Lines, LinearScale, LogScale
from bqplot.interacts import IndexSelector
from ipyleaflet import basemaps, FullScreenControl, LayerGroup, Map, MeasureControl, Polyline, Marker, CircleMarker, WidgetControl
from ipywidgets import Button, HTML, HBox, VBox, Checkbox, FileUpload, Label, Output, IntSlider, Layout, Image, link
# Only set to True for co-dev of ela from this use case:
ela_from_source = False
ela_from_source = True
if ela_from_source:
    if ('ELA_SRC' in os.environ):
        root_src_dir = os.environ['ELA_SRC']
    elif sys.platform == 'win32':
        root_src_dir = r'C:\src\github_jm\pyela'
    else:
        username = os.environ['USER']
        root_src_dir = os.path.join('/home', username, 'src/ela/pyela')
    pkg_src_dir = root_src_dir
    sys.path.insert(0, pkg_src_dir)

from ela.textproc import *
from ela.utils import *
from ela.classification import *
from ela.visual import *
from ela.spatial import SliceOperation<jupyter_output><empty_output><jupyter_text>Importing dataThere are two main sets of information we need: the borehole lithology logs, and the spatial information in the surface elevation (DEM) and geolocation of a subset of bores around Bungendore.<jupyter_code>data_path = None<jupyter_output><empty_output><jupyter_text>You probably want to explicitly set `data_path` to the location where you put the folder(s) e.g:<jupyter_code>#data_path = '/home/myusername/data' # On Linux, if you now have the folder /home/myusername/data/Bungendore
#data_path = r'C:\data\Lithology'  # windows, if you have C:\data\Lithology\Bungendore<jupyter_output><empty_output><jupyter_text>Otherwise a fallback for the pyela developer(s)<jupyter_code>if data_path is None:
    if ('ELA_DATA' in os.environ):
        data_path = os.environ['ELA_DATA']
    elif sys.platform == 'win32':
        data_path = r'C:\data\Lithology'
    else:
        username = os.environ['USER']
        data_path = os.path.join('/home', username, 'data')
data_path
data_path = '/home/per202/data/Lithology'
aem_datadir = os.path.join(data_path, 'AEM')
swan_datadir = os.path.join(data_path, 'swan_coastal')
scp_datadir = os.path.join(aem_datadir, 'Swan_coastal_plains')
scp_grids_datadir = os.path.join(scp_datadir, 'grids')
ngis_datadir = os.path.join(data_path, 'NGIS')
scp_shp_datadir = os.path.join(data_path, 'NGIS/swan_coastal')<jupyter_output><empty_output><jupyter_text>DEMLet's have a look at the DEM provided as part of the AEM package<jupyter_code>dem = rasterio.open(os.path.join(swan_datadir,'Swan_DEM/CLIP.tif'))
cnd_slice_dir = os.path.join(scp_grids_datadir,'cnd')

cnd_000_005 = rasterio.open(os.path.join(cnd_slice_dir,'Swan_Coastal_Plain_CND_000m_to_005m_Final.ers'))
bore_locations_raw = gpd.read_file(os.path.join(scp_shp_datadir, 'scp.shp'))<jupyter_output><empty_output><jupyter_text>The DEM raster and the bore location shapefile do not use the same projection (coordinate reference system) so we reproject one of them. We choose the raster's UTM.<jupyter_code>bore_locations = bore_locations_raw.to_crs(dem.crs)
import aseg_gdf2
gdf = aseg_gdf2.read( os.path.join(scp_datadir, 'located_data/Swan_Coastal_Plain_Final_CND'))
gdf
gdf.field_names()
df = gdf.df()
df.head()<jupyter_output><empty_output><jupyter_text>ViewerDerived from [voila-gpx-viewer](https://github.com/jtpio/voila-gpx-viewer)<jupyter_code># create the output widget to place the results
out = Output()
def plot_map(points):
    """
    Plot the GPS trace on a map
    """
    mean_lat = mean(p.Latitude for p in points)
    mean_lng = mean(p.Longitude for p in points)

    # create the map
    m = Map(center=(mean_lat, mean_lng), zoom=12, basemap=basemaps.Stamen.Terrain)

    # show trace
    line = Polyline(locations=[[[p.Latitude, p.Longitude] for p in points],],
                    color = "red", fill=False)
    m.add_layer(line)
    m.add_control(FullScreenControl())
    return m
ind = range(len(df))
points = [df.iloc[i] for i in ind if (i % 10 == 0)]
len(points)
gdf.get_field_definition('Date'), gdf.get_field_definition('Time')
p = points[123]
p.Date, p.Time
p = points[len(points) - 123]
p.Date, p.Time
set(p.Date for p in points)
points = [p for p in points if (p.Date == 20130512.0)]
points = sorted(points, key=lambda p: p.Time) 
def plot_elevation(gpx):
    px = [p.Time for p in points]
    py = [p.CND_011 for p in points]

    x_scale, y_scale = LinearScale(), LogScale()
    x_scale.allow_padding = False
    x_ax = Axis(label='Time (s)', scale=x_scale)
    y_ax = Axis(label='CND 011(?)', scale=y_scale, orientation='vertical')

    lines = Lines(x=px, y=py, scales={'x': x_scale, 'y': y_scale})

    elevation = Figure(title='CND 011 Chart', axes=[x_ax, y_ax], marks=[lines])
    elevation.layout.width = 'auto'
    elevation.layout.height = 'auto'
    elevation.layout.min_height = '500px'

    elevation.interaction = IndexSelector(scale=x_scale)

    return elevation

def link_trace_elevation(trace, elevation, points):
    """
    Link the trace the elevation graph.
    Changing the selection on the elevation will update the
    marker on the map
    """
    times = np.asarray([p.Time for p in points])

    def find_point(time):
        """
        Find a point given the time
        """
        dist_1 = abs(times - time)
        pos = np.argmin(dist_1)
        return points[pos]
    
    # add a checkbox to auto center
    autocenter = Checkbox(value=False, description='Auto Center')
    autocenter_control = WidgetControl(widget=autocenter, position='bottomright')
    trace.add_control(autocenter_control)
    # mark the current position on the map
    start = points[0]
    marker = CircleMarker(visible=False, location=(start.Latitude, start.Longitude),
                          radius=10, color="green", fill_color="green")
    trace.add_layer(marker)
    brushintsel = elevation.interaction
    def update_range(change):
        """
        Update the position on the map when the elevation
        graph selector changes
        """
        if brushintsel.selected.shape != (1,):
            return
        marker.visible = True
        selected = brushintsel.selected # time stamp in seconds for a day
        point = find_point(selected)
        marker.location = (point.Latitude, point.Longitude)
        if autocenter.value:
            trace.center = marker.location
        #position = max(0, int((selected / distance_from_start) * len(points)))
    brushintsel.observe(update_range, 'selected')

def plot_gpx(points):
    trace = plot_map(points)
    elevation = plot_elevation(points)
    debug = Label(value='')
    display(trace)
    display(elevation)
    display(debug)
    link_trace_elevation(trace, elevation, points)
plot_gpx(points)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup