updated the mpas plotting scripts

rrfs-test/IODA/offline_domain_check.py

@@ -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')
@@ -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')
 

rrfs-test/ush/mpasjedi_increment_fulldom.py

@@ -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":
@@ -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":
@@ -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
@@ -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)}")