-
Notifications
You must be signed in to change notification settings - Fork 0
/
Plot AHP.py
101 lines (79 loc) · 2.85 KB
/
Plot AHP.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 5 15:22:39 2022
@author: qvudinh
"""
# -*- coding: utf-8 -*-
"""
Created on Thu Mar 24 11:51:06 2022
@author: quang
"""
import pandas as pd
import string, os, sys
os.environ['PROJ_LIB'] = 'C:/Users/qvudinh/Anaconda3/Lib/site-packages/mpl_toolkits/basemap'
import netCDF4
import matplotlib.pyplot as plt
from wrf import to_np, getvar, smooth2d, get_basemap, latlon_coords
import numpy as np
from mpl_toolkits.basemap import Basemap
import xarray as xr
sys.path.append('../common/')
#from common import *
import cartopy.crs as ccrs
import matplotlib as mpl
import shapefile
def findIndexnear (lat_grid,lon_grid,lat_point, lon_point):
abslat = np.abs(lat_grid-lat_point)
abslon= np.abs(lon_grid-lon_point)
lat_id = np.argmin(abslat)
lon_id = np.argmin(abslon)
return lat_id, lon_id
dco = xr.open_dataset('AHP_150m.nc')
Ws = dco.AHP.values/1000
x, y = dco.lon.values, dco.lat.values
fsize = (6,7)
fig = plt.figure(figsize=fsize)
clevs = np.arange(16,56,2)
ax = plt.axes([0.1,0.18,0.8,0.8])
cmap = plt.get_cmap("rainbow")
norm = mpl.colors.BoundaryNorm(clevs, cmap.N)
m = Basemap(projection='merc',llcrnrlat=50.2,urcrnrlat=56.3,\
llcrnrlon=-12,urcrnrlon=-3.5,lat_ts=20,resolution='i')
#m.readshapefile(shp,'hanoi')
xx, yy = np.meshgrid(dco.lon.values, dco.lat.values)
x, y = m(xx, yy)
ax.contourf(x, y, Ws,levels = clevs,norm = mpl.colors.BoundaryNorm(clevs, cmap.N),cmap= cmap)
#m.readshapefile('haidao/haidao','haidao')
m.drawcoastlines(linewidth=0.25)
m.drawcountries(linewidth=0.25)
m.fillcontinents(color='paleturquoise',lake_color='1')
# draw parallels
m.drawparallels(np.arange(50,57,1),linewidth=0.,labels=[1,0,0,0],size=10)
# draw meridians
m.drawmeridians(np.arange(-12,-3,2),linewidth=0.,labels=[0,0,0,1],size=10)
df = pd.read_csv('sample point.csv')
lat_grid = dco.lat.values
lon_grid = dco.lon.values
j = 1
for i in range(len(df)):
latst = df.lat[i];
lonst = df.lon[i];
lat_id, lon_id = findIndexnear(lat_grid,lon_grid,latst,lonst)
x, y = m(lon_grid[lon_id], lat_grid[lat_id])
plt.scatter(x, y,s=50 , alpha=1,zorder=5, color='r')
plt.text(x,y,'Point '+str(j), fontsize = 10, color = 'k',weight = 'bold',zorder=5)
j = j+1
cb = mpl.colorbar.ColorbarBase(plt.axes([0.1, 0.11, 0.8,0.02]),
cmap=cmap,
norm = mpl.colors.BoundaryNorm(clevs, cmap.N),
orientation='horizontal',
#label='Wind speed (m/s)',fontsize=16,
#extend='both',
format='%d')
cb.ax.tick_params(labelsize=14)
cb.set_label(label='AHP (kt)',fontsize=15)
plt.show()
#if not os.path.exists('fig'): os.makedirs('fig')
#fig.savefig('AHP.png')
fig.savefig('AHP.svg', format = 'svg')
#fig.savefig('Figures/CF_datacorrection_2080_2099.png')