-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathzoom_korea_c48.py
84 lines (65 loc) · 2.22 KB
/
zoom_korea_c48.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
import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import os
datadir = '/work/noaa/da/weiwli/c48/gfs_land_c48/mem001/'
griddir = '/work/noaa/da/weiwli/c48/gfs_land_c48/grid/'
# one tile for now
tile=1
fname = os.path.join(datadir, f"20191215.180000.sfc_data.tile{tile}.nc")
oroname = os.path.join(griddir, f"C48_oro_data.tile{tile}.nc")
# read data and grid
ncdf = nc.Dataset(fname)
ncgf = nc.Dataset(oroname)
lat = ncgf.variables['geolat'][:]
lon = ncgf.variables['geolon'][:]
lon[np.where(lon>180)] = lon[np.where(lon>180)]-360
q2m = ncdf.variables['q2m'][0,:,:]
# set up map
#ax = plt.axes(projection=ccrs.PlateCarree())
# define grid
#gnomonic = ccrs.Gnomonic(central_latitude=lat[384,384], central_longitude=lon[384,384])
# plot data
#ax.pcolormesh(lon, lat, t127)
##ax.pcolormesh(lon, lat, t127, transform=gnomonic)
#ax.coastlines(zorder=10)
#ax.set_global()
#plt.show()
# loop and try all six and get the min/max values
arr3d=np.empty((6,48,48))
lat3d=np.empty((6,48,48))
lon3d=np.empty((6,48,48))
tilenum = [2,3,4]
for tile in tilenum:
#for tile in range(1,7,2):
print(tile)
fname = os.path.join(datadir, f"20191215.180000.sfc_data.tile{tile}.nc")
oroname = os.path.join(griddir, f"C48_oro_data.tile{tile}.nc")
# read data and grid
ncdf = nc.Dataset(fname)
ncgf = nc.Dataset(oroname)
lat = ncgf.variables['geolat'][0:48,0:48]
lon = ncgf.variables['geolon'][0:48,0:48]
lon[np.where(lon>180)] = lon[np.where(lon>180)]-360
q2m = ncdf.variables['q2m'][0,:,:]
# put these 2D arrays in the 3D arrays
arr3d[tile-1,...] = q2m
lat3d[tile-1,...] = lat
lon3d[tile-1,...] = lon
# get min/max for plotting
minval = np.nanmin(arr3d)
maxval = np.nanmax(arr3d)
# set up map
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([120, 135, 30, 45])
# loop and plot data
cmaps=['jet', 'brg', 'terrain', 'gnuplot', 'coolwarm', 'RdGy']
itilenum = [1,2,3]
for itile in itilenum:
print(itile)
#ax.pcolormesh(lon3d[itile], lat3d[itile], arr3d[itile], vmin=minval, vmax=maxval)
ax.pcolormesh(lon3d[itile], lat3d[itile], arr3d[itile], vmin=minval, vmax=maxval, cmap=cmaps[itile])
ax.coastlines()
plt.show()