Earth example and interactive plotting



In [1]:

    
%load_ext autoreload
%autoreload 2
%matplotlib inline



In [2]:

    
from os import path, environ
import fiona
import numpy as N
import attitude
from attitude import Orientation, ReconstructedPlane, create_groups
from attitude.display import plot_interactive, init_notebook_mode
from attitude.plot import plot_aligned
from json import dumps
from sys import argv
import palettable as P



In [3]:

    
BASE = path.join(path.dirname(attitude.__file__),'..','js-frontend')
DATA = path.join(BASE,"ui-test/data")

Two views of the stratigraphy



In [4]:

    
from IPython.display import Image
Image(filename=path.join(DATA,'DJI_0454.PNG'), width=500)









    Out[4]:



In [5]:

    
Image(filename=path.join(DATA,'DJI_0455.PNG'), width=500)









    Out[5]:

We're now focusing on the two cliffs in the bottom of the image

Step 1. Import drone-measured elevations

We exported points laid onto the 3d model from photoscan



In [6]:

    
fn = path.join(DATA,'elevations_v2.dxf')

orientations = []
with fiona.open(fn) as ds:
    for i,item in ds.items():
        coords = N.array(item['geometry']['coordinates'])
        try:
            orientations.append(Orientation(coords))
        except ValueError:
            continue

disabled = ["afb1b38f","ebab5ecb","7b12f5d4","72417f52","8565a276","939d9b9a",
            "e884eac1", "c8902997","15026bc4","730a6124","eb31c4b0","d3793eb9"]
orientations = [o for o in orientations if o.hash not in disabled]

groups = (
    ["3df53944","e98e12a9"],
    ["493a240a","04987fe8"],
    ["c2122b2a","de785eb5"]
)

groupedOrientations = create_groups(orientations, *groups,
                             same_plane=False)

collection = [a.to_mapping(color='#ff0000', type='remote')
              for a in groupedOrientations]

Step 2. Import field-measured orientations



In [72]:

    
fn = path.join(DATA,'field-orientations.geojson')

with fiona.open(fn) as ds:
    for i,item in ds.items():
        p = item['properties']
        if p['planeType'].strip() != 'Bedding':
            continue

        asm = p.get("aster_smoothed")
        alt = asm

        alt -= 40 # Global datum is higher than local
        center = (*item['geometry']['coordinates'],alt)

        err = 0.1*N.pi/180
        a = ReconstructedPlane(p['strike'], p['dip'],0,err,err)
        orientation = a.to_mapping(
            center=center,
            color='#444', type='in-situ')
        collection.append(orientation)

removedUIDs = ["89636280","6031fd6f"]
collection = [c for c in collection if 1600 < c['center'][2] < 1680]
collection = [c for c in collection if c['uid'] not in removedUIDs]



In [73]:

    
cmap = P.mycarta.CubeYF_7
cmap.show_as_blocks()



In [74]:

    
cmap.hex_colors









    Out[74]:





['#7B0290', '#7D48EA', '#5A89EE', '#39BAB1', '#51D866', '#8AEB51', '#D1EB5B']



In [75]:

    
heights = N.array([o['center'][2] for o in collection])
rng = [heights.min(),heights.max()]

for o in collection:
    ix = N.interp(o['center'][2], rng, [0,6])
    o['color'] = cmap.hex_colors[6-int(ix)]



In [77]:

    
init_notebook_mode()
plot_interactive(collection)



In [ ]: