#!/usr/bin/env python
import numpy as np
import pyproj
import cst

# parameters
path = ''
meta = cst.util.load( 'ws-meta.py' )
bounds = (0.0, 125.0), (0.0, 250.0)
extent = meta.extent
proj = pyproj.Proj( **meta.projection )
proj = cst.coord.Transform( proj, **meta.transform )

# hypocenter
x, y, z = np.loadtxt( 'src/hypocenter.txt' ).T
z *= 0.001
x, y = proj( x, y )
np.savetxt( path + 'hypocenter-xyz.txt', np.array( [[x,y,z]] ) )

# topo
topo, extent = cst.data.topo( extent, scale=0.001 )
ddeg = 0.5 / 60.0

# stations
f = np.loadtxt( 'src/station-list.txt', 'S8, f, f, f' )
x, y = proj( f['f2'], f['f1'] )
z = cst.coord.interp2( extent, topo, (x, y) )
np.savetxt( path + 'stations-xyz.txt', np.array( [x,y,z] ).T )

# map data
x, y = cst.data.mapdata( 'coastlines', 'high', extent, 10.0 )
z = cst.coord.interp2( extent, topo, (x, y) )
x, y = proj( x, y )
x, y, i = cst.data.clipdata( x, y, bounds )
z = z[i]
np.savetxt( path + 'coastlines-xyz.txt', np.array( [x,y,z] ).T )

# mesh
n = meta.x_shape[:2]
n = n[1], n[0], 2
x, y = np.loadtxt( 'src/grid-1.0x1.0-gp1_h0.125.ll.gz', usecols=(0,1) ).reshape( n ).T
z = cst.coord.interp2( extent, topo, (x, y) )
x.astype( 'f' ).T.tofile( path + 'lon.bin' )
y.astype( 'f' ).T.tofile( path + 'lat.bin' )
z.astype( 'f' ).T.tofile( path + 'topo.bin' )

# limits
corners = (
    (x[0,0], x[-1,0], x[-1,-1], x[0,-1]),
    (y[0,0], y[-1,0], y[-1,-1], y[0,-1]),
)
print( 'corners = %r, %r' % corners )

# basins
v = 2500.0,
z.fill( 1000.0 )
z = cst.cvm.extract( x, y, z, 'vs', rundir='run/cvm' )
x, y = proj( x, y )
bounds = (
    (x[0].mean(), x[-1].mean()),
    (y[:,0].mean(), y[:,-1].mean()),
)
print( 'bounds = %r, %r' % bounds )
x, y = cst.plt.contour( x, y, z, v )[0]
z = np.empty_like( x )
z.fill( -1000.0 )
np.savetxt( path + 'basins-xyz.txt', np.array( [x,y,z] ).T )