# -*- coding: utf-8 -*-
"""
Read a directory of MISMIP3d grounding line file,
and create a netcdf file similar to the MISMIP+ format but
containing only the x and y co-oordinates of the GL point

assumes the subdirectory mismip3d_dma6 exists and
contains the mimsip3d/DMA6*b05.dat files at 10 year intervals

Created on Tue Sep  8 09:29:09 2015

@author: s.l.cornford@bris.ac.uk
"""


import numpy as np
from netCDF4 import Dataset
import glob
import csv
import os
import matplotlib.pyplot as plt

def readm3d(path):
    """
    read a mismip 3d grounding line file.
    return numpy arrays x and y co-ordinates of all points

    The MIMSIP3d grounding line format is a plain text file,
    containing the data for one point on each line as
    space separated values:

    x y thickness ubase vbase usurface vsurface umean vmeam butressing_factor tangential_stress
    """

    csv.register_dialect('dat', delimiter=' ', quoting=csv.QUOTE_NONE)

    #first pass to work out number of rows
    nrow = 0
    with open(path) as csvfile:
        reader = csv.reader(csvfile, delimiter=' ', quotechar='"')
        for row in reader:
            nrow = nrow  + 1

    bufx = np.empty(nrow)
    bufy = np.empty(nrow)

    irow = 0
    with open(path) as csvfile:
        reader = csv.reader(csvfile, delimiter=' ', quotechar='"')
        for row in reader:
            bufx[irow] = np.float32(row[0])
            bufy[irow] = np.float32(row[1])
            irow = irow + 1

    return bufx, bufy

def parsetime(path):
    """
    work out the time from a mismip3d file name.
    """
    def flnm(s):
        try:
            return float(s)
        except ValueError:
            return -1.0

    def swtime(exp):
        return {
            'Stnd' : 0.0,
            'P75S' : 0.0 + flnm(fnam[8:11]),
            'P75R' : 100.0 + flnm(fnam[8:11]),
        }.get(exp, -1.0)

    fnam = os.path.basename(path)

    return swtime(fnam[4:8])


nc = Dataset('mismip3d_example.nc', 'w')

glptdim = nc.createDimension('nPointGL', size=None)
timedim = nc.createDimension('nTime', size=21) # every 10 years
xgl = nc.createVariable('xGL', 'f4', ('nPointGL', 'nTime'))
ygl = nc.createVariable('yGL', 'f4', ('nPointGL', 'nTime'))
time = nc.createVariable('time', 'f4', ('nTime'))

itime = 0
pathlist = sorted(glob.glob("mismip3d_dma6/DMA6Stnd*"))
pathlist.extend(sorted(glob.glob("mismip3d_dma6/DMA6P75S*")))
pathlist.extend(sorted(glob.glob("mismip3d_dma6/DMA6P75R*")))


for p in pathlist:

    x, y = readm3d(p)
    xgl[0:len(x), itime] = x
    ygl[0:len(x), itime] = y
    time[itime] = parsetime(p)
    itime = itime + 1


nc.close()

#test the file : open it and plot the gl

def glplot(ncfile, times, colora, label):
    """
    add a plot of grounding line points to current axes.
    makes use of the numpy.ma.MaskedArray when reading xGL,yGL
    """
    ncid = Dataset(ncfile, 'r')
    ltime = ncid.variables["time"][:]
    lxmax = 0.0
    lxmin = 800.0
    for i in range(0, len(times)):
        seq = (ltime == times[i])
        xGL = ncid.variables["xGL"][:, seq]*1e-3
        lxmax = max(np.max(xGL), lxmax)
        lxmin = min(np.min(xGL), lxmin)
        yGL = ncid.variables["yGL"][:, seq]*1e-3
        plt.plot(xGL, yGL, 's', ms=3, mfc=colora[i],
                 mec=colora[i], label=label + ', t = ' + format(times[i]))
    ncid.close()
    return lxmin, lxmax


xmin, xmax = glplot('mismip3d_example.nc', [0, 20, 100, 120, 200],
                    ['black', 'pink', 'red', 'purple', 'blue'], 'P75')

plt.legend(loc='upper right', ncol=1, frameon=True,
           borderaxespad=0, fontsize='small')

plt.savefig("mismip3d_plot_example.pdf")

