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)
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='')
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)
gallp = gall * g1
hvplot.save(gallp, '1pointlocation.html')
# 左岸杭座標
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] )
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)
gallp = gall * g2
hvplot.save(gallp, 'multipointlocation.html')