#!python3
# plot eke from a vtk file
import sys

prfx='/home/abozec/PYTHON/'
sys.path.append(prfx)
import myenv as my
import meshio
from scipy.interpolate import griddata

# define bb86 20km grid
res=0.20 # resolution of bb86 20km
idm=101 ; jdm=101
ini_lon=10.0 ; ini_lat=35.0
lon_bb86=my.np.zeros([jdm,idm])
lat_bb86=my.np.zeros([jdm,idm])

lon_bb86[:, 0] = ini_lon
for i in my.np.arange(idm-1)+1:
   lon_bb86[:,i] = lon_bb86[:, i-1] + res

 ## Latitude
lat_bb86[0, :] = ini_lat
for j in my.np.arange(jdm-1)+1:
   lat_bb86[j, :] = lat_bb86[j-1, :] + res

# directory
io='/nexsan/people/abozec/NUMO/MICHAL/VTK_REPORT2/'
mesh=meshio.read(io+'EKE_20km_new.vtk')

# get list of attibutes of the object
#vars(mesh) #or dir(mesh)
# eke in point_data
# get x,y coordinate
points=mesh.points[:,0:2]
n=points.shape[0]
##get coordinate of input grid in degrees for griddata
points[:,0]=points[:,0]/(100. * 1.0e3)+ini_lon
points[:,1]=points[:,1]/(100. * 1.0e3)+ini_lat

# get EKE from dictionary
eke= mesh.point_data.get('EKE')

## interpolate model to bb86
eke_bb86=griddata(points,eke[:,0],(lon_bb86,lat_bb86),method='linear')


# colormap
cmapeke=my.mygc.get_cmapeke()
N=cmapeke.colors.shape[0]
cmap=cmapeke

vmin=0.;vmax=300.
levels_to_draw=my.np.linspace(vmin,vmax,N+1)
# plot the field
fig, axs = my.plot.subplots(nrows=1,axwidth='4.5in')
axs[0].format(title='GNuME EKE (cm2/s2) 20km  Y:10')
m = axs[0].contourf(eke_bb86, cmap=cmap, extend='neither',\
                    levels=levels_to_draw)
axs[0].colorbar(m,loc='b',ticks=50.)

ps_dir='/nexsan/people/abozec/NUMO/MICHAL/VTK_REPORT2/PS/'
file_ps='GNuME_eke_20km_y10.pdf'
print(file_ps)
fig.savefig(ps_dir+file_ps,dpi=150,\
            facecolor='w', edgecolor='w',transparent=False) ## .pdf,.eps

my.plot.show() ## interactive (block=False), otherwise ()