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

# # import module

# In[1]:


import holoviews as hv
import geoviews as gv
import numpy as np
import pandas as pd
import pyproj
import bokeh
import xarray as xr
from cartopy import crs

hv.extension('bokeh')
gv.extension('bokeh')


# In[2]:


print( 'holoviews-version  ' + hv.__version__ )
print( 'geoviews-version  ' + gv.__version__ )
print( 'bokeh-version  ' + bokeh.__version__ )


# # reference
# - QuadMesh : http://holoviews.org/reference/elements/matplotlib/QuadMesh.html
# 
# - Curvilinear Grids : http://geoviews.org/user_guide/Resampling_Grids.html
# 
# - about xarray(in japanese) :
# https://qiita.com/fujiisoup/items/0e9b2e74902551aa2855
# 
# - xarray of irregular mesh :
# http://xarray.pydata.org/en/stable/plotting.html
# 
# - netCDF export:http://xarray.pydata.org/en/stable/io.html

# # make geomap
#  - WMTS maps of Geospatial Information Authority of Japan

# In[3]:


outcrs = crs.epsg(6671)
url = 'https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/{Z}/{X}/{Y}.jpg'
geomap = gv.WMTS(url, crs=outcrs).options(width=600, height=400)


# # make  vector of $\xi$ axis and $\eta$ axis
# 
#  - origin coordinate : longitude, latitude = 132.8143, 35.4159
#  - An angle of xi axis from horizontal axis is 26 degrees
#  - $\xi$ axis length is 1000m, $\eta$ axis length is 90m
#  - coordinate system is EPSG:6668 in geographic, or EPSG:6671 in projected

# In[4]:


p1 = [132.8143, 35.4159]

EPSGin = pyproj.Proj("+init=EPSG:6668")
EPSGout = pyproj.Proj("+init=EPSG:6671")
p1x = np.array( [pyproj.transform( EPSGin, EPSGout, *p1)] )

def mkunitvecotr(th):
    ex = np.array( [[1, 0]] )
    th = th * np.pi / 180
    r = np.array( [[np.cos(th), -np.sin(th)],[np.sin(th), np.cos(th)]])
    return np.dot( r, ex.T)

thx = 26
exi = mkunitvecotr(thx)
eeta = mkunitvecotr(thx+90)

Lx, Ly = 1000, 90
pathxi = np.c_[p1x.T, p1x.T + exi*Lx]
patheta = np.c_[p1x.T, p1x.T + eeta*Ly]

# drawing using geoviews
l1 = gv.Path(tuple(pathxi), crs=outcrs).options(line_width=2)
l2 = gv.Path(tuple(patheta), crs=outcrs).options(line_width=2)

out = geomap*l1*l2
gv.renderer('bokeh').save(out, 'set_axis')
out


# # make mesh points 
# - set coordinates of x_index and yindex 
# - ndarray [x_index, y_index, [x_coordinate, y_coordinate]]
# - mesh size is dx:5m, dy:5m

# In[5]:


dx, dy = 5, 5
mps = np.empty(( int(Lx/dx)+1, int(Ly/dy)+1, 2 ))
for i, _ in enumerate(mps[:,0]):
    for j, _ in enumerate(mps[0]):
        mps[i,j] = (p1x.T + exi*dx*i + eeta*dy*j).flatten()


# # setting elevation data a mesh point each
# 
# - using libgsidem2el https://github.com/computational-sediment-hyd/libgsidem2el

# In[6]:


import libgsidem2el as gsi
dem = gsi.libgsidem2el('DEM5A')

EPSGout = pyproj.Proj("+init=EPSG:6668")
EPSGin = pyproj.Proj("+init=EPSG:6671")

Z = np.empty_like( mps[:,:,0] ) 
nx, ny = len(Z), len(Z[0])
for j in range(ny) :
    for i in range(nx) :
        p = mps[i, j]
        lon, lat = pyproj.transform( EPSGin, EPSGout, *p)
        Z[i,j] = dem.getEL(lon, lat, zoom = 15)


# # draw mesh elevation data

# In[7]:


Qx = mps[:,:,0]
Qy = mps[:,:,1]
qmesh = gv.QuadMesh((Qx, Qy, Z), crs=outcrs).options( colorbar=True, cmap='jet',alpha=0.5)
out = geomap*qmesh
gv.renderer('bokeh').save(out, 'mesh_elevataion')
out


# # set normalized elevation by elevation average bed slope
#  - averaged bed slop is 1/670

# In[8]:


org = []
for i in range(len(Z[0])): 
    org.append(hv.Curve((Z[:,i]) ,vdims='Z'))

ZZ = np.empty_like(Z)
xarr = np.arange(len(Z[:,i]))*dx
d = 5.5 - xarr * 1/670
for i in range(len(Z[0])): ZZ[:,i] = d
ZZZ = Z - ZZ

norm = []
for i in range(len(ZZZ[0])):
    norm.append(hv.Curve((ZZZ[:,i]), vdims='dz') )

out = (hv.Overlay(org).relabel('Elevation') +  hv.Overlay(norm).relabel('Normalized_Elevation')).options(width=600, height=400)
gv.renderer('bokeh').save(out, 'normalized')
out


# In[9]:


Qx = mps[:,:,0]
Qy = mps[:,:,1]
qmesh = gv.QuadMesh((Qx, Qy, ZZZ), crs=outcrs).options( colorbar=True, cmap='jet',alpha=0.5)
out = geomap*qmesh
out


# # make xarray and export by NetCDF format

# In[10]:


ds = xr.Dataset({'dz': (['x','y'], ZZZ)}, coords={'xc': (('x', 'y'), Qx), 'yc': (('x', 'y'), Qy)})
ds.to_netcdf('dz.nc')