holoviz-topics / EarthSim

Tools for working with and visualizing environmental simulations.
https://earthsim.holoviz.org
BSD 3-Clause "New" or "Revised" License
65 stars 21 forks source link

vectror plot and colorbar issue #321

Closed sameerCoder closed 4 years ago

sameerCoder commented 4 years ago

datamat = loadmat('IOEC_ECM_noDA_20190703_masked.mat') Xp = datamat['Xp'] Yp = datamat['Yp'] strt = datetime.datetime(2019, 7, 3, 0, 0) end = datetime.datetime(2019, 1, 13, 0, 0)

def perdelta(strt, end, delta): curr = strt while curr < end: yield curr curr += delta

Read element file

data = pd.read_table('fort.ele',delim_whitespace=True,names=('A','B','C','D'))

tri_new = pd.read_csv('fort.ele', delim_whitespace=True, names=('A', 'B', 'C', 'D'), usecols=[1, 2, 3], skiprows=1, dtype={'D': np.int})

data1=data[['B','C','D']]

tri=data1[1:]

dateList = [] for result in perdelta(strt, strt + timedelta(days=1), timedelta(hours=3)): dat = result

print(result)

datetostr=result.strftime("%Y%b%d%H")
dt = parse(str(dat))
yr = dt.year
mn = dt.month
d = dt.day
hr = dt.hour
mi = dt.minute
# print(y,mn,d,hr,mi)
if hr < 10:
    # d='0'+str(d)
    hr = '0' + str(hr)
else:
    d = str(d)
    hr = str(hr)
if int(d) < 10:
    d = '0' + str(d)
else:
    d = str(d)
varname = 'Hsig_' + str(yr) + '0' + str(mn) + str(d) + '_' + hr + '0000'
print(varname)

x = Xp.flatten()
y = Yp.flatten()
z = datamat[varname]
z = z.flatten()
Longitude=x
Latitude=y
HS=z
#z=z.round()
print(z)
#break

pts = np.stack((x,y,z)).T

verts = pd.DataFrame(np.stack((x,y,z)).T, columns=['x', 'y' , 'z'])

pts = np.stack((Longitude,Latitude,HS)).T
verts = pd.DataFrame(np.stack((Longitude,Latitude,HS)).T, columns=['Longitude', 'Latitude' , 'HS'])

pkname = 'PkDir_' + str(yr) + '0' + str(mn) + str(d) + '_' + hr + '0000'

pkvalue = datamat[pkname]
pksub=pkvalue[1:200000]
d1 = pkvalue * (math.pi / 180)
#d1 = pksub * (math.pi / 180)

print(d1)
print(len(d1))

u = np.cos(d1) * 1
v = np.sin(d1) * 1

u = u.flatten()
v = v.flatten()

print("length of u ", len(u))
print("lenght of v ", len(v))

um, vm = np.meshgrid(u, v, sparse=True)
print(um)
um = um.flatten()

vm = vm.flatten()

# Below opts line is for plot image size , wmts change size .
#%opts WMTS[width=900 height=900]

tiles=gv.WMTS('https://maps.wikimedia.org/osm-intl/{Z}/{X}/{Y}@2x.png')

#openStreet Background.
tri_sub = tri_new.apply(lambda x: x - 1)

levels=[0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.6,1.8,2,2.2,2.5,3,]
colorbar_opts={
    'major_label_overrides': {
        0:'0',

        0.5:'0.5',

        1:'1',

        1.5:'1.5',

        2:'2',

        2.5:'2.5',
        3:'3',

        3.5:'3.5',
        3.8:' ',
        3.9:'>3.9',

    },
    'major_label_text_align': 'left','major_label_text_font_style':'bold',}

def plot_limits(plot, element):
    plot.handles['x_range'].min_interval = 100
    plot.handles['x_range'].max_interval = 55000000
    #3000000
    plot.handles['y_range'].min_interval = 500
    plot.handles['y_range'].max_interval = 900000
opts = dict(width=700, height=700, logz=False, tools=['hover','save','wheel_zoom'],hooks=[plot_limits],colorbar=True,color_levels=15,
            colorbar_opts=colorbar_opts,cmap= ['#000080',
                                               '#0000cd',

'#0008ff', '#004cff', '#0090ff', '#00d4ff', '#29ffce', '#60ff97', '#97ff60', '#ceff29', '#ffe600', '#ffa700',

'#ff6800',

                                               '#ff2900',
                                               '#cd0000',
                                               '#800000',

],clim=(0,3.76),title=" \t\t\t\tSignificant Wave Height (m) ",fontsize={'title': 26, 'xlabel':15, 'ylabel':15, 'ticks':12} )

pkname = 'PkDir_' + str(yr) + '0' + str(mn) + str(d) + '_' + hr + '0000'

#pkvalue = datamat[pkname]
print(pkname)
pkvalue = datamat[pkname]
pkvalue=pkvalue.flatten()
d=pkvalue*(math.pi/180)

#z=z[:,:5]
#break

# target grid to interpolate to
xt = np.arange(76.937,92.008,0.1)
yt = np.arange(1.482,22.461,0.1)
xi,yi = np.meshgrid(xt,yt)
di = griddata((Longitude,Latitude),d,(xi,yi))

dfcoast = pd.read_csv('ECpolygonTwoDegreeOffsetBorder.txt', delim_whitespace=True, names=('X', 'Y'))
dfcoast['geometry'] = dfcoast.apply(lambda row: Point(row.X, row.Y), axis=1)
poly = Polygon([(p.x, p.y)  for p in  dfcoast.geometry])

point = Point(87.3148, 15.6066)

poly.contains(point)

arr=np.zeros((len(yt),len(xt)))
for i in range(len(xt)):
    for j in range(len(yt)):
        point=Point(xt[i],yt[j])
        arr[j,i]=poly.contains(point)

#di[~arr] = np.nan

mask = (xi > 79.7817) & (xi < 81.2718) & (yi > 7.6951) & (yi < 9.7406)
di[mask] = np.nan
di[arr==False]=np.nan
U=np.cos(di)
V=np.sin(di)

mag = np.sqrt(U**2 + V**2)
angle = (np.pi/2.) - np.arctan2(U/mag, V/mag)
tiles = gv.tile_sources.Wikipedia
#vectorfield = gv.VectorField((xi[::5,::5], yi[::5,::5], angle[::5,::5], mag[::5,::5]))

ggpoints=gv.Points(verts,vdims=['HS'])

gg=tiles * rasterize(gv.TriMesh((tri_sub,gv.Points(verts)) )).options(**opts)*gv.VectorField((xi[::15,::15], yi[::15,::15], angle[::15,::15], mag[::15,::15]))
#gg=tiles * rasterize(gv.TriMesh((tri_sub,gv.Points(verts)) ), aggregator=ds.mean('z')).options(**opts)

display(gg)
gv.save(gg,datetostr+"aggregator=ds.mean_out8july.html")

show(gg)

print("End2")
print("End")

break