mpas plotting#
House keeping and load modules#
%%time
# autoload external python modules if they changed
%load_ext autoreload
%autoreload 2
import sys, os
pyDAmonitor_ROOT=os.getenv("pyDAmonitor_ROOT")
if pyDAmonitor_ROOT is None:
print("!!! pyDAmonitor_ROOT is NOT set. Run `source ush/load_pyDAmonitor.sh`")
else:
print(f"pyDAmonitor_ROOT={pyDAmonitor_ROOT}\n")
sys.path.insert(0, pyDAmonitor_ROOT)
# import modules
import warnings
import cartopy.crs as ccrs
import geoviews.feature as gf
import holoviews as hv
import hvplot.xarray
from holoviews import opts
import uxarray as ux
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
hv.extension("bokeh")
# hv.extension("matplotlib")
import geoviews as gv
import cartopy.crs as ccrs
import geopandas as gp
pyDAmonitor_ROOT=/home/Guoqing.Ge/pyDAmonitor
CPU times: user 3.88 s, sys: 1.19 s, total: 5.07 s
Wall time: 6.25 s
from DAmonitor.shapes import coast_lines, state_lines, county_lines
from DAmonitor.mpas import hcross_contour, hcross_contour0, vcross_contour
load MPAS data#
%%time
grid_file=f"{pyDAmonitor_ROOT}/data/mpasjedi/invariant.nc"
bkg_file=f"{pyDAmonitor_ROOT}/data/mpasjedi/bkg.nc"
ana_file=f"{pyDAmonitor_ROOT}/data/mpasjedi/ana.nc"
# conus12km example:
# grid_file="../data/mpasjedi/conus12km/conus12km.invariant.nc_L60_GFS"
# bkg_file="../data/mpasjedi/conus12km/bkg/mpasout.2024-05-06_01.00.00.nc"
# ana_file="../data/mpasjedi/conus12km/ana/mpasout.2024-05-06_01.00.00.nc"
# haidao's cris-fsr_n20 run (55 levels)
# grid_file="../data/mpasjedi/cris-fsr_n20/invariant.nc"
# bkg_file="../data/mpasjedi/cris-fsr_n20/mpasout.2024-05-27_00.00.00.nc"
# ana_file="../data/mpasjedi/cris-fsr_n20/ana.2024-05-27_00.00.00.nc"
uxds_a = ux.open_dataset(grid_file, ana_file)
uxds_b = ux.open_dataset(grid_file, bkg_file)
CPU times: user 1.88 s, sys: 3.52 s, total: 5.4 s
Wall time: 5.64 s
make a horizontal temperature plot#
%%time
uxvar0 = uxds_a['theta'] - 273.15
uxvar = uxvar0
nt = 0 # time dimension
lev = 0 # vertical level
## plot = hcross_contour0(uxvar.isel(Time=nt, nVertLevels=lev), title=f'lev={lev}', cincr=5) # mostly for difference plots
plot = hcross_contour(uxvar.isel(Time=nt, nVertLevels=lev), title=f'lev={lev}') #, symmetric_cmap=True)
plot * coast_lines * state_lines
CPU times: user 7.62 s, sys: 1.47 s, total: 9.09 s
Wall time: 9.22 s
setting up a subdomain box#
lon_center = -105.03
lat_center = 39.0
lon_incr = 5 # degree
lat_incr = 3 # degree
lon_bounds = (lon_center - lon_incr, lon_center + lon_incr)
lat_bounds = (lat_center - lat_incr, lat_center + lat_incr)
plot temperature analysis increments at different levels#
%%time
var_name = "theta"
uxdiff0 = uxds_a[var_name] - uxds_b[var_name]
uxvar = uxdiff0
nt=0 # time dimension
plot_levels = [30] # [0, 29, 42] # [0, 19, 29, 39, 49, 58]
# Create the colormap
# colors = ['blue', 'white', 'red']
# cmap = LinearSegmentedColormap.from_list('blue_white_red', colors)
# zero_shift = 0.02
cmap = plt.get_cmap('coolwarm')
zero_shift = 0.0
plots = []
for lev in plot_levels:
### use hcross_contour0() which uses 0 divide the cool and warm colors in the plot by default
#
# tmp = hcross_contour0(
# uxvar.isel(Time=nt, nVertLevels=lev),
# title=f'lev={lev}',
# cmap=cmap,
# zero_shift=zero_shift,
# clevs_multiplier=1,
# ) # for the whole domain
# hcross_contour() does not dvide the cool and warm colors at 0 by default
# But it can be achieved by setting symmetric_cmap=True which will set symmetric cmax/cmin automatically,
# or mannualy setting symmetric cmax/cmin
#
tmp = hcross_contour(
uxvar.isel(Time=nt, nVertLevels=lev),
title=f'lev={lev}',
symmetric_cmap=True,
#clevs=20,
) # for the whole domain
plots.append(tmp * coast_lines * state_lines)
# plots share one toolbar, which facilitates doing sync'ed zoom-in/out
# hv.Layout(plots).cols(1)
# each plot has its own toolbar, which facilitates controlling each plot individually
for p in plots:
display(p)
CPU times: user 2.52 s, sys: 711 ms, total: 3.23 s
Wall time: 3.25 s
plot temperature increments zoomed into the Colorado domain#
%%time
var_name = "theta"
uxdiff0 = uxds_a[var_name] - uxds_b[var_name]
### subset to a small domain
uxdiff1 = uxdiff0.subset.bounding_box(lon_bounds, lat_bounds,)
uxvar = uxdiff1
nt=0 # time dimension
plot_levels = [42] # [0, 29, 42] # [0, 19, 29, 39, 49, 58]
plots = []
for lev in plot_levels:
tmp = hcross_contour(uxvar.isel(Time=nt, nVertLevels=lev), title=f'lev={lev}', width=700, height=500) # for the subdomain
# overlay state_lines
#plots.append(tmp * coast_lines * state_lines)
# overlay county lines, this takes longer time to render
plots.append(tmp * coast_lines * county_lines.opts(xlim=(lon_bounds[0], lon_bounds[1]), ylim=(lat_bounds[0], lat_bounds[1])))
# plots share one toolbar, which facilitates doing sync'ed zoom-in/out
# hv.Layout(plots).cols(1)
# each plot has its own toolbar, which facilitates controlling each plot individually
for p in plots:
display(p)