import numpy as np
import xarray as xr
import hvplot.xarray
import cartopy.crs as crs, geoviews as gv
import holoviews as hv
import glob
gv.extension('bokeh', logo=False)
hv.extension('bokeh', logo=False)
ds = xr.open_dataset('out00000024.nc')
ds.depth.values[ ds.depth.values < 0.01 ] = np.nan
outcrs = crs.epsg(6677)
url = 'https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/{Z}/{X}/{Y}.jpg'
geomap = gv.WMTS(url, crs=outcrs) #.options(width=600, height=400)
out1 = ds.hvplot.quadmesh( 'xc','yc',z='depth'
, crs=outcrs, rasterize=True, dynamic=True, cmap='jet'
,width=600, height=800,colorbar=True ).redim.range(depth=(0.2,1.5))
x0 = ds.xc.values[0,0]
x1 = ds.xc.values[1,0]
y0 = ds.yc.values[0,0]
y1 = ds.yc.values[1,0]
vc = [x1-x0, y1-y0]
th = np.arccos( vc[0]/np.linalg.norm(vc) )
U = ds.u.values[::3,::3]
V = ds.v.values[::3,::3]
X = ds.xc.values[::3,::3]
Y = ds.yc.values[::3,::3]
mag = np.sqrt(U**2 + V**2) + 0.00000001
angle = (np.pi/2.) - np.arctan2(U/mag, V/mag) - th
vectorfield = gv.VectorField((X, Y, angle, mag), crs=outcrs).options(magnitude='Magnitude', scale=0.5, color='white')
g = out1.opts(alpha=0.5) * vectorfield * geomap
hv.save(g,'out.html')