#!/usr/bin/env python
# coding: utf-8

# # OpendTect points to seismic subimages
# 
# 
# ## Read OdT data
# 
# First we want to read the OdT pointset files.

# In[1]:


fname = "/home/matt/Downloads/salt.pck"

with open(fname, 'r') as f:
    data = f.read()
    
print(data[:300])


# This seems to be:
# 
# - UTMx, UTMy, TWTT, something, something, something

# In[4]:


import numpy as np

def read_odt_pts(fname):
    data = np.loadtxt(fname, skiprows=6, comments='!', usecols=[0,1,2])
    return data


# In[3]:


# Example
d = read_odt_pts(fname)
d


# ## Read Petrel data
# 
# I wrote this a while back; it's from [`gio`](https://github.com/agile-geoscience/gio)

# In[4]:


from io import StringIO

def read_petrel_points(filename):
    """
    Read a Petrel points file. Return an array.
    TODO:
        Do something with Comments and Fields.
    """
    with open(filename) as f:

        comments = []
        fields = []
        in_header = False

        while True:
            line = f.readline().strip()

            if line.startswith('#'):
                comments.append(line.strip('# '))
            elif line.startswith('VERSION'):
                # version = line.split()[-1]
                pass
            elif line.startswith('BEGIN'):
                in_header = True
            elif line.startswith('END'):
                in_header = False
                break
            elif in_header:
                fields.append(line.strip())
            else:
                break

        d = f.read()
        s = StringIO(d)

        return np.loadtxt(s)


# In[5]:


d = read_petrel_points('/home/matt/Downloads/fault.pts')


# In[6]:


d


# ## Transform coordinates from UTM to inline/crossline
# 
# There's an affine transformation class in Agile's [`bruges`](https://github.com/agile-geoscience/bruges). 
# 
# You'll have to do this in your environment before you can import it:
# 
#     pip install bruges

# In[7]:


f3_corners_xy = np.array([
    [605835.5, 6073556.5],
    [629576.25, 6074220.0],
    [629122.5, 6090463.2]
   ])

f3_corners_ix = np.array([[0,  0],
                          [0, 950],
                          [650, 950]
                         ])

import bruges as bg

transform = bg.transform.CoordTransform(f3_corners_ix, f3_corners_xy)


# In[8]:


# Example:
transform.reverse(d[0, :2])


# In[9]:


# Put it all together
def pts_to_ixt(fname, dt=0.004):
    data = read_odt_pts(fname)
    transform = bg.transform.CoordTransform(f3_corners_ix, f3_corners_xy)
    ix = [transform.reverse(r) for r in data[:, :2]]
    t = (data[:, 2][:, None] / dt).astype(int)
    ixt = np.hstack([ix, t])
    return ixt

ixt = pts_to_ixt("/home/matt/Downloads/salt.pck")


# In[10]:


ixt


# These are the (inline, crossline, timeslice) coordinates of the (UTMx, UTMy, TWTT) points in the pointset.

# ## Extract data
# 
# Let's start with a small example of random data.

# In[11]:


# This is 278MB in memory - random integers in (0, 255).
f3 = np.random.randint(0, 256, size=(650, 950, 450), dtype=np.uint8)


# In[12]:


samples = []
for (i, x, t) in ixt:
    sample = f3[i-3:i+3, x-3:x+3, t]
    samples.append(sample)


# In[13]:


get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt

fig, axs = plt.subplots(nrows=2, ncols=5, figsize=(15, 6))
for ax, sample in zip(axs.ravel(), samples[:10]):
    ax.imshow(sample)


# ## Unsupervised dimensionality reduction on seismic data
# 
# In case the full volume is not available (eg from [Data Underground](https://dataunderground.org)), we can use a smaller version of F3 from the Geocomputing class.

# In[14]:


import numpy as np


# f3 = np.load("/home/matt/Dropbox/dev/geocomp-19/data/F3_volume_3x3_16bit.npy")


# In[15]:


import segyio

with segyio.open("/home/matt/Downloads/F3_entire.segy") as s:
    f3 = segyio.cube(s)


# In[17]:


f3.shape


# In[16]:


# Normalize to 1
f3 = f3 / np.amax(np.abs(f3))


# This volume is over 1GB:

# In[33]:


f3.nbytes/1e9


# In[35]:


np.save('/home/matt/F3_entire.npy', f3)


# In[25]:


ma = np.percentile(f3[100], 99.8)

plt.figure(figsize=(15, 8))
plt.imshow(f3[10, :, :].T, aspect='auto', vmin=-ma, vmax=ma)


# Let's do a quick experiment on a small subvolume. We're going to divide it into 3 × 3 subcubes, so the sides have to be divisible by 3.

# In[27]:


f3_ = f3[10:31, 10:31, 15:450]
f3_.shape


# 'Tile' (in 3d) the volume. From https://stackoverflow.com/questions/42297115/numpy-split-cube-into-cubes

# In[28]:


import numpy as np

def cubify(arr, newshape):
    oldshape = np.array(arr.shape)
    repeats = (oldshape / newshape).astype(int)
    tmpshape = np.column_stack([repeats, newshape]).ravel()
    order = np.arange(len(tmpshape))
    order = np.concatenate([order[::2], order[1::2]])
    # newshape must divide oldshape evenly or else ValueError will be raised
    return arr.reshape(tmpshape).transpose(order).reshape(-1, *newshape)


# In[29]:


X = cubify(f3_, (3, 3, 3)).reshape(-1, 27)


# For loading into https://projector.tensorflow.org/ for example...

# In[30]:


np.savetxt('cubelets.tsv', X, delimiter='\t')


# ## UMAP
# 
# Let's try UMAP first:
# 
#     pip install umap-learn

# In[31]:


import umap

reducer = umap.UMAP()

embedding = reducer.fit_transform(X)


# In[32]:


get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt

plt.figure(figsize=(15, 10))
plt.scatter(*embedding.T)


# ## tSNE
# 
# This is part of `sklearn`:

# In[84]:


import numpy as np
from sklearn.manifold import TSNE

embedding_tsne = TSNE().fit_transform(X)
embedding_tsne.shape


# In[88]:


plt.figure(figsize=(15, 10))
plt.scatter(*embedding_tsne.T)


# In[ ]: