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')
print( 'holoviews-version ' + hv.__version__ )
print( 'geoviews-version ' + gv.__version__ )
print( 'bokeh-version ' + bokeh.__version__ )
holoviews-version 1.11.2 geoviews-version 1.6.2 bokeh-version 1.0.4
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
http://xarray.pydata.org/en/stable/plotting.html
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)
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
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()
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)
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
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
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
ds = xr.Dataset({'dz': (['x','y'], ZZZ)}, coords={'xc': (('x', 'y'), Qx), 'yc': (('x', 'y'), Qy)})
ds.to_netcdf('dz.nc')