Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

updated the mpas plotting scripts #174

Merged
merged 3 commits into from
Sep 18, 2024
Merged

updated the mpas plotting scripts #174

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
e467655
updated the mpas plotting scripts (the domain check also has plotting…
delippi Sep 17, 2024
f58b58f
remove tri.simplices and make triangulation part of scripts more simi…
delippi Sep 17, 2024
efbcbec
revert to using restart.nc instead of mpasout.nc until the mpas case …
delippi Sep 18, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 18 additions & 5 deletions rrfs-test/IODA/offline_domain_check.py
View file
Open in desktop
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
106 changes: 44 additions & 62 deletions rrfs-test/ush/mpasjedi_increment_fulldom.py
View file
Open in desktop
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,32 @@

# 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/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 +70,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 +98,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
tri = Triangulation(lons, lats)
tri_mask = TriAnalyzer(tri).get_flat_tri_mask(min_circle_ratio=0.1)
tri.set_mask(tri_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(tri, jedi_inc, clevs, cmap=cm, extend='both')

# 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
Loading