Skip to content

Commit

Permalink
updated the mpas plotting scripts (NOAA-EMC#174)
Browse files Browse the repository at this point in the history 
Updated the mpas plotting scripts (the domain check also has plotting
script built in) to correctly plot mpas domain edges.
  • Loading branch information
@delippi
delippi authored Sep 18, 2024
1 parent bf67744 commit 22c7702
Show file tree
Hide file tree
Showing 4 changed files with 168 additions and 161 deletions.
23 changes: 18 additions & 5 deletions rrfs-test/IODA/offline_domain_check.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,23 +59,26 @@ def toc(tic=tic, label=""):
parser.add_argument('-g', '--grid', type=str, help='grid file', required=True)
parser.add_argument('-o', '--obs', type=str, help='ioda observation file', required=True)
parser.add_argument('-f', '--fig', action='store_true', help='disable figure (default is False)', required=False)
parser.add_argument('-s', '--shrink', type=float, help='hull shrink factor', required=True)
args = parser.parse_args()

# Assign filenames
obs_filename = args.obs
grid_filename = args.grid # see note above.
make_fig = args.fig
hull_shrink_factor = args.shrink

print(f"Obs file: {obs_filename}")
print(f"Grid file: {grid_filename}")
print(f"Figure flag: {args.fig}")
print(f"Hull shrink factor: {hull_shrink_factor}")

# Plotting options
plot_box_width = 100. # define size of plot domain (units: lat/lon degrees)
plot_box_height = 50
cen_lat = 34.5
cen_lon = -97.5
hull_shrink_factor = 0.04 #4% was found to work fairly well.
#hull_shrink_factor = 0.10 #10% was found to work fairly well.

grid_ds = nc.Dataset(grid_filename, 'r')
obs_ds = nc.Dataset(obs_filename, 'r')
Expand Down Expand Up @@ -234,24 +237,34 @@ def shrink_boundary(points, centroid, factor=0.04):
gl1.ylabel_style = {'size': 5, 'color': 'gray'}

# Plot the domain and the observations
m1.fill(adjusted_lon.flatten(), grid_lat.flatten(), color='b', label='Domain Boundary', zorder=1, transform=ccrs.PlateCarree())
#m1.fill(adjusted_lon.flatten(), grid_lat.flatten(), color='b', label='Domain Boundary', zorder=1, transform=ccrs.PlateCarree())
m1.scatter(adjusted_lon.flatten(), grid_lat.flatten(), c='b', s=1, label='Domain Boundary', zorder=2)
m1.plot(adjusted_shrunken_points[:, 0], shrunken_points[:, 1], 'r-', label='Convex Hull', zorder=10, transform=ccrs.PlateCarree())

# Plot included observations
included_lat = obs_lat[inside_indices]
included_lon = obs_lon[inside_indices]
plt.scatter(included_lon, included_lat, c='g', s=2, label='Included Observations', zorder=3, transform=ccrs.PlateCarree())
included_count = len(included_lat)
plt.scatter(included_lon, included_lat, c='g', s=2, label=f'Included Observations ({included_count})', zorder=3, transform=ccrs.PlateCarree())

# Plot excluded observations
excluded_indices = np.setdiff1d(np.arange(len(obs_lat)), inside_indices)
excluded_lat = obs_lat[excluded_indices]
excluded_lon = obs_lon[excluded_indices]
plt.scatter(excluded_lon, excluded_lat, c='r', s=2, label='Excluded Observations', zorder=4, transform=ccrs.PlateCarree())

excluded_count = len(excluded_lat)
total_count = len(obs_lat)

print(f"Ob counts:")
print(f" Excluded: {excluded_count}")
print(f" Included: {included_count}")
print(f" Total: {total_count}")
plt.scatter(excluded_lon, excluded_lat, c='r', s=2, label=f'Excluded Observations ({excluded_count})', zorder=4, transform=ccrs.PlateCarree())

plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.legend(loc='upper right')
plt.title(f'{dycore} Domain and Observations')
plt.title(f'{dycore} Domain and Observations ({hull_shrink_factor*100}%)')
plt.tight_layout()
plt.savefig(f'./domain_check_{dycore}.png')

Expand Down
107 changes: 45 additions & 62 deletions rrfs-test/ush/mpasjedi_increment_fulldom.py
View file
Original file line number Diff line number Diff line change
@@ -1,32 +1,31 @@
#!/usr/bin/env python
from netCDF4 import Dataset
import pdb
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import cartopy.geodesic
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
import numpy as np
import colormap
import time
import sys, os
import shapely.geometry
import os
import warnings
from scipy.spatial.distance import cdist
from matplotlib.tri import Triangulation, TriAnalyzer

warnings.filterwarnings('ignore')

############ USER INPUT ##########################################################
plot_var = "Increment"
lev = 23 # 1=sfc, 55=toa
clevmax_incr = 2.0 # max contour level for colorbar increment plots
lev = 1 # 1=sfc, 55=toa
clevmax_incr = 2.0 # max contour level for colorbar increment plots
decimals = 2 # number of decimals to round for text boxes
plot_box_width = 75. # define size of plot domain (units: lat/lon degrees)
plot_box_height = 30
plot_box_width = 100. # define size of plot domain (units: lat/lon degrees)
plot_box_height = 50.
cen_lat = 34.5
cen_lon = -97.5

variable = "airTemperature"
if variable == "airTemperature":
Expand All @@ -35,36 +34,33 @@

# JEDI data
datapath = "./"
jstatic = "./data/restart.2024-05-27_00.00.00.nc" # to load the MPAS lat/lon
jstatic = "./data/static.nc" # to load the MPAS lat/lon

# FOR HYBRID (or ENVAR)
janalysis = "./ana.2024-05-27_00.00.00.nc" # analysis file
jbackgrnd = "./data/restart.2024-05-27_00.00.00.nc" # background file (control member)

# FOR LETKF
#janalysis = "./ana.2024-05-27_00.00.00.nc" # analysis file
#jbackgrnd = "./bkg.2024-05-27_00.00.00.nc" # background file (ensmean)
janalysis = "./ana.2024-05-27_00.00.00.nc" # analysis file
jbackgrnd = "./data/restart.2024-05-27_00.00.00.nc" # background file
#jbackgrnd = "./data/mpasout.2024-05-27_00.00.00.nc" # background file

###################################################################################
# Set cartopy shapefile path
platform = os.getenv('HOSTNAME').upper()
if 'ORION' in platform:
cartopy.config['data_dir']='/work/noaa/fv3-cam/sdegelia/cartopy'
cartopy.config['data_dir']='/work/noaa/fv3-cam/sdegelia/cartopy'
elif 'H' in platform: # Will need to improve this once Hercules is supported
cartopy.config['data_dir']='/home/Donald.E.Lippi/cartopy'
cartopy.config['data_dir']='/home/Donald.E.Lippi/cartopy'

f_latlon = Dataset(jstatic, "r")
lats = np.array( f_latlon.variables['latCell'][:] ) * 180.0 / np.pi
lons0 = np.array( f_latlon.variables['lonCell'][:] ) * 180.0 / np.pi
lons = np.where(lons0>180.0,lons0-360.0,lons0)
lats = np.array(f_latlon.variables['latCell'][:]) * 180.0 / np.pi
lons0 = np.array(f_latlon.variables['lonCell'][:]) * 180.0 / np.pi
lons = np.where(lons0 > 180.0, lons0 - 360.0, lons0)

# Open NETCDF4 dataset for reading
nc_a = Dataset(janalysis, mode='r')
nc_b = Dataset(jbackgrnd, mode='r')

# Read data and get dimensions
jedi_a = nc_a.variables["theta"][0,:,lev].astype(np.float64)
jedi_b = nc_b.variables["theta"][0,:,lev].astype(np.float64)
lev = lev - 1 # subtract 1 because python uses indices starting from 0
jedi_a = nc_a.variables["theta"][0, :, lev].astype(np.float64)
jedi_b = nc_b.variables["theta"][0, :, lev].astype(np.float64)

# Convert to temperature
if variable == "airTemperature":
Expand All @@ -75,21 +71,23 @@
jedi_a = jedi_a / dividend_a
jedi_b = jedi_b / dividend_b

# compute increment
# Compute increment
jedi_inc = jedi_a - jedi_b

if plot_var == "Increment":
title1 = "JEDI"
jedi = jedi_inc
clevmax = clevmax_incr
def plot_T_inc(var_n, clevmax):
"""Temperature increment/diff [K]"""
longname = "airTemperature"
units = "K"
inc = 0.05 * clevmax
clevs = np.arange(-1.0 * clevmax, 1.0 * clevmax + inc, inc)
cm = colormap.diff_colormap(clevs)
return clevs, cm, units, longname

# CREATE PLOT ##############################
fig = plt.figure(figsize=(7,4))
fig = plt.figure(figsize=(7, 4))
m1 = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree(central_longitude=0))

# Determine extent for plot domain
cen_lat = 35.5
cen_lon = -99.0
half = plot_box_width / 2.
left = cen_lon - half
right = cen_lon + half
Expand All @@ -101,58 +99,43 @@
m1.set_extent([left, right, top, bot])

# Add features to the subplots
#m1.add_feature(cfeature.GSHHSFeature(scale='low'))
m1.add_feature(cfeature.COASTLINE)
m1.add_feature(cfeature.BORDERS)
m1.add_feature(cfeature.STATES)
#m.add_feature(cfeature.OCEAN)
#m.add_feature(cfeature.LAND)
#m.add_feature(cfeature.LAKES)

# Gridlines for the subplots
gl1 = m1.gridlines(crs = ccrs.PlateCarree(), draw_labels = True, linewidth = 0.5, color = 'k', alpha = 0.25, linestyle = '-')
gl1.xlocator = mticker.FixedLocator([])
gl1 = m1.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.5, color='k', alpha=0.25, linestyle='-')
gl1.xlocator = mticker.FixedLocator(np.arange(-180., 181., 5.))
gl1.ylocator = mticker.FixedLocator(np.arange(-80., 91., 5.))
gl1.xformatter = LONGITUDE_FORMATTER
gl1.yformatter = LATITUDE_FORMATTER
gl1.xlabel_style = {'size': 5, 'color': 'gray'}
gl1.ylabel_style = {'size': 5, 'color': 'gray'}

# Create triangulation and mask
triang = Triangulation(lons, lats)
mask = TriAnalyzer(triang).get_flat_tri_mask(min_circle_ratio=0.1)
triang.set_mask(mask)

def plot_T_inc(var_n, clevmax):
"""Temperature increment/diff [K]"""
longname = "airTemperature"
units="K"
inc = 0.05*clevmax
clevs = np.arange(-1.0*clevmax, 1.0*clevmax+inc, inc)
cm = colormap.diff_colormap(clevs)
return(clevs, cm, units, longname)

# Plot the data
if variable == "airTemperature":
clevs, cm, units, longname = plot_T_inc(jedi_inc, clevmax)

units="K"
#c1 = m1.contourf(lons, lats, jedi, clevs, cmap = cm)
c1 = plt.tricontourf(lons, lats, jedi, clevs, cmap = cm, extend='both')
# Plot the data using triangulation
clevs, cm, units, longname = plot_T_inc(jedi_inc, clevmax_incr)
c1 = m1.tricontourf(triang, jedi_inc, clevs, cmap=cm, extend='both', transform=ccrs.PlateCarree())

# Add colorbar
cbar1 = fig.colorbar(c1, orientation="horizontal", fraction=0.046, pad=0.07)
cbar1.set_label(units, size=8)
cbar1.ax.tick_params(labelsize=5, rotation=30)

# Add titles, text, and save the figure
#plt.suptitle(f"Temperature {plot_var} at Level: {lev+1}\nobtype: {longname}", fontsize = 9, y = 1.05)
#m1.set_title(f"{title1}", fontsize = 9, y = 0.98)
#subtitle1_minmax = f"min: {np.around(np.min(jedi), decimals)}\nmax: {np.around(np.max(jedi), decimals)}"
#m1.text(left, top, f"{subtitle1_minmax}", fontsize = 6, ha='left', va='bottom')
#m1.text(left, bot, f"{subtitle1_hofx}", fontsize = 6, ha='left', va='bottom')
# Add 1 to final lev since indicies start from 0
plt.suptitle(f"Temperature {plot_var} at Level: {lev+1}", fontsize=9, y=0.95)
subtitle1_minmax = f"min: {np.around(np.min(jedi_inc), decimals)}\nmax: {np.around(np.max(jedi_inc), decimals)}"
m1.text(left * 0.99, bot * 1.01, f"{subtitle1_minmax}", fontsize=6, ha='left', va='bottom')
if plot_var == "Increment":
plt.tight_layout()
plt.savefig(f"./increment_{variable}.png", dpi=250, bbox_inches='tight')

# Print some final stats
print(f"Stats:")
print(f" {title1} max: {np.around(np.max(jedi), decimals)}")
print(f" {title1} min: {np.around(np.min(jedi), decimals)}")
print(f" {longname} max: {np.around(np.max(jedi_inc), decimals)}")
print(f" {longname} min: {np.around(np.min(jedi_inc), decimals)}")
Loading

0 comments on commit 22c7702

Please sign in to comment.