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

# ### background

# In[1]:


import numpy as np
import geopandas as gpd
import xarray as xr
import cartopy.crs as crs
import hvplot.pandas
import hvplot.xarray
import geoviews as gv
from shapely.geometry import Point

gv.extension('bokeh', logo=False)


# In[2]:


ds = xr.open_rasterio('04gg542_0.5g.tif')

z = ds.values
z = np.where(z==-9999, np.nan, z)
ds.values = z/100

nepsg = int(ds.attrs['crs'][-4:])

g = ds.isel(band=0).hvplot.image(x='x',y='y', cmap='terrain',geo=True
                    , dynamic=False, rasterize=True, colorbar=False, alpha=1).options(clipping_colors={'NaN':'transparent'})

back = gv.WMTS('https://mt1.google.com/vt/lyrs=s&x={X}&y={Y}&z={Z}', name="GoogleMapsImagery" 
               , crs=crs.epsg(nepsg))
# .redim( Longitude={'range':(ds.x.min(), ds.x.max())}, Latitude={'range':(ds.y.min(), ds.y.max())})

gall = (back*g).options(xlabel='', ylabel='')


# ### plot 1 point

# In[3]:


pl = [98891.63181796, 104393.37756069]
gdf = gpd.GeoDataFrame(geometry=[Point(pl[0], pl[1])], crs={'init': 'epsg:' + str(nepsg)})

g1 = gdf.hvplot.points(dynamic=False, geo=True, project=True, crs=crs.epsg(nepsg),size=50, color='red', hover=False)


# In[4]:


gallp = gall * g1
hvplot.save(gallp, '1pointlocation.html')


# ### plot multi points

# In[5]:


# 左岸杭座標
pl = [98891.63181796, 104393.37756069]
# 右岸杭座標
pr = [ 99438.53777557, 104214.89229045]

pl = np.array(pl)
pr = np.array(pr)

# 測線の長さ
L = np.linalg.norm(pr - pl)
# 単位ベクトルe
e = (pr - pl)/ L

# 杭の外側（堤内地側）にLex延伸する。
Lex = 50
# 測線上にdelta間隔の点群を発生させる。
delta = 0.5
dl = np.arange(-Lex, L+Lex, delta)
points = np.array( [pl+e*dlp for dlp in dl] )


# In[6]:


gdf = gpd.GeoDataFrame(geometry=[Point(p[0],p[1]) for p in points], crs={'init': 'epsg:' + str(nepsg)})
g2 = gdf.hvplot.points(dynamic=False, geo=True, project=True, crs=crs.epsg(nepsg),size=50, color='red', hover=False)


# In[7]:


gallp = gall * g2
hvplot.save(gallp, 'multipointlocation.html')