From e467655e2777524ad4d351425fb787223756ac72 Mon Sep 17 00:00:00 2001 From: "Donald.E.Lippi" Date: Tue, 17 Sep 2024 18:01:09 +0000 Subject: [PATCH 1/3] updated the mpas plotting scripts (the domain check also has plotting script built in) to correctly plot mpas domain edges. Prior to this update the concave boundaries would be filled by the plotting routines incorrectly leading to potential misinterpretation of data. --- rrfs-test/IODA/offline_domain_check.py | 23 +++- rrfs-test/ush/mpasjedi_increment_fulldom.py | 106 ++++++++----------- rrfs-test/ush/mpasjedi_increment_singleob.py | 100 ++++++++--------- rrfs-test/ush/mpasjedi_spread.py | 90 ++++++++++------ 4 files changed, 162 insertions(+), 157 deletions(-) mode change 100644 => 100755 rrfs-test/ush/mpasjedi_increment_singleob.py diff --git a/rrfs-test/IODA/offline_domain_check.py b/rrfs-test/IODA/offline_domain_check.py index 3d995d8..1be7992 100755 --- a/rrfs-test/IODA/offline_domain_check.py +++ b/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') diff --git a/rrfs-test/ush/mpasjedi_increment_fulldom.py b/rrfs-test/ush/mpasjedi_increment_fulldom.py index d906ee3..edefb34 100755 --- a/rrfs-test/ush/mpasjedi_increment_fulldom.py +++ b/rrfs-test/ush/mpasjedi_increment_fulldom.py @@ -1,10 +1,9 @@ +#!/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 @@ -12,21 +11,21 @@ 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,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": @@ -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 @@ -101,17 +98,12 @@ 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 @@ -119,23 +111,14 @@ 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) @@ -143,16 +126,15 @@ def plot_T_inc(var_n, clevmax): 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)}") diff --git a/rrfs-test/ush/mpasjedi_increment_singleob.py b/rrfs-test/ush/mpasjedi_increment_singleob.py old mode 100644 new mode 100755 index bada944..5264c27 --- a/rrfs-test/ush/mpasjedi_increment_singleob.py +++ b/rrfs-test/ush/mpasjedi_increment_singleob.py @@ -1,10 +1,9 @@ +#!/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 @@ -12,18 +11,16 @@ 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 = 8. # define size of plot domain (units: lat/lon degrees) @@ -34,12 +31,12 @@ # JEDI data datapath = "./" -jstatic = "./data/restart.2024-05-27_00.00.00.nc" # to load the MPAS lat/lon -jdiag = "adpupa_hofx.nc4" # obs diag file +jstatic = "./data/static.nc" # to load the MPAS lat/lon +jdiag = "adpupa_hofx.nc4" # obs diag file # 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) +jbackgrnd = "./data/mpasout.2024-05-27_00.00.00.nc" # background file (control member) # FOR LETKF #janalysis = "./ana.2024-05-27_00.00.00.nc" # analysis file @@ -50,15 +47,15 @@ # 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' print(f"{jdiag}") # Do a quick omf and hofx value check from GSI diag file # FROM JEDI diag jncdiag = Dataset(jdiag, mode='r') -if 'restart' in jbackgrnd: +if 'mpasout' in jbackgrnd: # Assuming Var run oberr_input = jncdiag.groups["EffectiveError0"].variables[f"{variable}"][:][0] oberr_final = jncdiag.groups["EffectiveError2"].variables[f"{variable}"][:][0] @@ -80,7 +77,7 @@ print(f"JEDI:") print(f" oberr_input: {oberr_input}") print(f" oberr_final: {oberr_final}") -subtitle1_hofx = f" ob: {job}\n omf: {jomf}\n hofx: {jhofx}\n oberr_final: {oberr_final}\n" +subtitle1_hofx = f"ob: {job}\nomf: {jomf}\nhofx: {jhofx}\noberr_final: {oberr_final}" print(subtitle1_hofx) # Grab single ob lat/lon and double check they're the same. @@ -98,21 +95,21 @@ print(f" singleob_lon: {singleob_lon}\n") print(f"Creating {plot_var} Plot...\n") -lev = lev - 1 # JEDI #################################### 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": @@ -123,16 +120,20 @@ 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=(3,3)) +fig = plt.figure(figsize=(3, 3)) m1 = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree(central_longitude=0)) # Determine extent for plot domain @@ -147,17 +148,12 @@ 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 @@ -165,23 +161,14 @@ 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) +# 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') # Scatter the single ob location m1.scatter(singleob_lon, singleob_lat, color='g', marker='o', s=2) @@ -192,16 +179,17 @@ def plot_T_inc(var_n, clevmax): 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)}" +subtitle1 = f"{subtitle1_hofx}\n{subtitle1_minmax}" +m1.text(left*0.995, bot*1.005, f"{subtitle1}", 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)}") diff --git a/rrfs-test/ush/mpasjedi_spread.py b/rrfs-test/ush/mpasjedi_spread.py index 277ddd6..2973bbc 100755 --- a/rrfs-test/ush/mpasjedi_spread.py +++ b/rrfs-test/ush/mpasjedi_spread.py @@ -1,3 +1,4 @@ +#!/usr/bin/env python from netCDF4 import Dataset import matplotlib matplotlib.use('agg') @@ -5,82 +6,103 @@ import cartopy import cartopy.crs as ccrs import cartopy.feature as cfeature -from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER import numpy as np -import time import sys, os import warnings +from matplotlib.tri import Triangulation, TriAnalyzer +from scipy.spatial import Delaunay warnings.filterwarnings('ignore') ############ USER INPUT ########################################################## # MPAS info -jstatic = "data/static.nc" # to load the MPAS lat/lon -ens_file = "restart.2024-05-27_00.00.00.nc" +jstatic = "data/static.nc" # to load the MPAS lat/lon +ens_file = "mpasout.2024-05-27_00.00.00.nc" nmems = 30 # Plotting options -variable = "airTemperature" # currently only working for [airTemperature] +variable = "airTemperature" # currently only working for [airTemperature] lev = 23 # 1=sfc, 55=toa -cen_lat = 35.44 -cen_lon = -97.47 -plot_box_width = 60. # 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 ################################################################################### # Set cartopy shapefile path platform = os.getenv('HOSTNAME').upper() if 'ORION' in platform: - 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']='/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' # Get domain info 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) # Now read the var you want bg_all = [] -for imem in range(1, nmems+1): - infile = 'data/ens/mem%s/%s' % (str(imem).zfill(2), ens_file) - nc = Dataset(infile, 'r') - if variable == 'airTemperature': - bg = nc.variables['theta'][0,:,lev].astype(np.float64) - pres = (nc.variables['pressure_p'][0,:,lev] + nc.variables['pressure_base'][0,:,lev])/100.0 - dividend = (1000.0/pres)**(0.286) - bg = bg / dividend - else: - sys.exit("Error: requested variable not yet implemented.") - bg_all.append(bg) +for imem in range(1, nmems + 1): + infile = 'data/ens/mem%s/%s' % (str(imem).zfill(3), ens_file) + nc = Dataset(infile, 'r') + if variable == 'airTemperature': + bg = nc.variables['theta'][0, :, lev].astype(np.float64) + pres = (nc.variables['pressure_p'][0, :, lev] + nc.variables['pressure_base'][0, :, lev]) / 100.0 + dividend = (1000.0 / pres) ** (0.286) + bg = bg / dividend + else: + sys.exit("Error: requested variable not yet implemented.") + bg_all.append(bg) # Compute variance bg_all = np.asarray(bg_all) var = np.nanvar(bg_all, axis=0) -# CREATE PLOT ############################## -fig = plt.figure(figsize=(7,4)) -m1 = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree(central_longitude=0)) +# Construct Delaunay triangulation +points = np.vstack((lons, lats)).T +tri = Delaunay(points) + +# Define plot domain bounds half = plot_box_width / 2. left = cen_lon - half right = cen_lon + half half = plot_box_height / 2. bot = cen_lat - half top = cen_lat + half -m1.set_extent([left, right, top, bot]) + +# Create plot ############################## +fig = plt.figure(figsize=(7, 4)) +m1 = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree(central_longitude=0)) + +m1.set_extent([left, right, bot, top]) + m1.add_feature(cfeature.COASTLINE) m1.add_feature(cfeature.BORDERS) m1.add_feature(cfeature.STATES) -units=r"$K^2$" -levs = np.arange(0, 2.1, 0.1) -c1 = plt.tricontourf(lons, lats, var, levs, cmap = 'viridis', extend='both') + +# Define colormap +units = r"$K^2$" +cmap = plt.get_cmap('viridis') +cmap.set_under('gray') + +# Plot using Delaunay triangulation +triang = Triangulation(lons, lats, tri.simplices) +mask = TriAnalyzer(triang).get_flat_tri_mask(min_circle_ratio=0.1) +triang.set_mask(mask) + +# Plot the variance data on the triangulation +c1 = m1.tricontourf(triang, var, cmap=cmap, levels=100, 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) + plt.tight_layout() plt.savefig(f"./spread_{variable}.png", dpi=250, bbox_inches='tight') plt.close() + From f58b58f3e3f969da018fed96a109529b649723c4 Mon Sep 17 00:00:00 2001 From: "Donald.E.Lippi" Date: Tue, 17 Sep 2024 19:05:06 +0000 Subject: [PATCH 2/3] remove tri.simplices and make triangulation part of scripts more similar. --- rrfs-test/ush/mpasjedi_increment_fulldom.py | 8 ++++---- rrfs-test/ush/mpasjedi_increment_singleob.py | 8 ++++---- rrfs-test/ush/mpasjedi_spread.py | 7 +------ 3 files changed, 9 insertions(+), 14 deletions(-) diff --git a/rrfs-test/ush/mpasjedi_increment_fulldom.py b/rrfs-test/ush/mpasjedi_increment_fulldom.py index edefb34..031d7b5 100755 --- a/rrfs-test/ush/mpasjedi_increment_fulldom.py +++ b/rrfs-test/ush/mpasjedi_increment_fulldom.py @@ -112,13 +112,13 @@ def plot_T_inc(var_n, clevmax): 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) +triang = Triangulation(lons, lats) +mask = TriAnalyzer(triang).get_flat_tri_mask(min_circle_ratio=0.1) +triang.set_mask(mask) # 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') +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) diff --git a/rrfs-test/ush/mpasjedi_increment_singleob.py b/rrfs-test/ush/mpasjedi_increment_singleob.py index 5264c27..c193a5c 100755 --- a/rrfs-test/ush/mpasjedi_increment_singleob.py +++ b/rrfs-test/ush/mpasjedi_increment_singleob.py @@ -162,13 +162,13 @@ def plot_T_inc(var_n, clevmax): 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) +triang = Triangulation(lons, lats) +mask = TriAnalyzer(triang).get_flat_tri_mask(min_circle_ratio=0.1) +triang.set_mask(mask) # 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') +c1 = m1.tricontourf(triang, jedi_inc, clevs, cmap=cm, extend='both', transform=ccrs.PlateCarree()) # Scatter the single ob location m1.scatter(singleob_lon, singleob_lat, color='g', marker='o', s=2) diff --git a/rrfs-test/ush/mpasjedi_spread.py b/rrfs-test/ush/mpasjedi_spread.py index 2973bbc..4961264 100755 --- a/rrfs-test/ush/mpasjedi_spread.py +++ b/rrfs-test/ush/mpasjedi_spread.py @@ -10,7 +10,6 @@ import sys, os import warnings from matplotlib.tri import Triangulation, TriAnalyzer -from scipy.spatial import Delaunay warnings.filterwarnings('ignore') @@ -62,10 +61,6 @@ bg_all = np.asarray(bg_all) var = np.nanvar(bg_all, axis=0) -# Construct Delaunay triangulation -points = np.vstack((lons, lats)).T -tri = Delaunay(points) - # Define plot domain bounds half = plot_box_width / 2. left = cen_lon - half @@ -90,7 +85,7 @@ cmap.set_under('gray') # Plot using Delaunay triangulation -triang = Triangulation(lons, lats, tri.simplices) +triang = Triangulation(lons, lats) mask = TriAnalyzer(triang).get_flat_tri_mask(min_circle_ratio=0.1) triang.set_mask(mask) From efbcbec61f1b3ca2262c446b1231726cb3753585 Mon Sep 17 00:00:00 2001 From: "Donald.E.Lippi" Date: Wed, 18 Sep 2024 15:25:01 +0000 Subject: [PATCH 3/3] revert to using restart.nc instead of mpasout.nc until the mpas case has been fully updated. --- rrfs-test/ush/mpasjedi_increment_fulldom.py | 3 ++- rrfs-test/ush/mpasjedi_increment_singleob.py | 17 +++++++++-------- rrfs-test/ush/mpasjedi_spread.py | 13 +++++++++---- 3 files changed, 20 insertions(+), 13 deletions(-) diff --git a/rrfs-test/ush/mpasjedi_increment_fulldom.py b/rrfs-test/ush/mpasjedi_increment_fulldom.py index 031d7b5..80c6992 100755 --- a/rrfs-test/ush/mpasjedi_increment_fulldom.py +++ b/rrfs-test/ush/mpasjedi_increment_fulldom.py @@ -37,7 +37,8 @@ jstatic = "./data/static.nc" # to load the MPAS lat/lon janalysis = "./ana.2024-05-27_00.00.00.nc" # analysis file -jbackgrnd = "./data/mpasout.2024-05-27_00.00.00.nc" # background 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 diff --git a/rrfs-test/ush/mpasjedi_increment_singleob.py b/rrfs-test/ush/mpasjedi_increment_singleob.py index c193a5c..292d092 100755 --- a/rrfs-test/ush/mpasjedi_increment_singleob.py +++ b/rrfs-test/ush/mpasjedi_increment_singleob.py @@ -35,12 +35,13 @@ jdiag = "adpupa_hofx.nc4" # obs diag file # FOR HYBRID (or ENVAR) -janalysis = "./ana.2024-05-27_00.00.00.nc" # analysis file -jbackgrnd = "./data/mpasout.2024-05-27_00.00.00.nc" # background file (control member) +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) +#jbackgrnd = "./data/mpasout.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 = "./bkg.2024-05-27_00.00.00.nc" # background file (ensmean) ################################################################################### @@ -55,7 +56,7 @@ # FROM JEDI diag jncdiag = Dataset(jdiag, mode='r') -if 'mpasout' in jbackgrnd: +if 'mpasout' or 'restart' in jbackgrnd: # Assuming Var run oberr_input = jncdiag.groups["EffectiveError0"].variables[f"{variable}"][:][0] oberr_final = jncdiag.groups["EffectiveError2"].variables[f"{variable}"][:][0] @@ -65,9 +66,9 @@ oberr_input = jncdiag.groups["ObsError"].variables[f"{variable}"][:][0] oberr_final = jncdiag.groups["ObsError"].variables[f"{variable}"][:][0] ob = jncdiag.groups["ObsValue"].variables[f"{variable}"][:][0] + offset -omf= jncdiag.groups["ombg"].variables[f"{variable}"][:][0] -fmo= -1*omf -hofx= fmo+ob +omf = jncdiag.groups["ombg"].variables[f"{variable}"][:][0] +fmo = -1 * omf +hofx = fmo + ob oberr_input = np.around(oberr_input.astype(np.float64), decimals) oberr_final = np.around(oberr_final.astype(np.float64), decimals) diff --git a/rrfs-test/ush/mpasjedi_spread.py b/rrfs-test/ush/mpasjedi_spread.py index 4961264..36f86fc 100755 --- a/rrfs-test/ush/mpasjedi_spread.py +++ b/rrfs-test/ush/mpasjedi_spread.py @@ -17,7 +17,8 @@ # MPAS info jstatic = "data/static.nc" # to load the MPAS lat/lon -ens_file = "mpasout.2024-05-27_00.00.00.nc" +ens_file = "restart.2024-05-27_00.00.00.nc" +#ens_file = "mpasout.2024-05-27_00.00.00.nc" nmems = 30 # Plotting options @@ -43,11 +44,15 @@ lons0 = np.array(f_latlon.variables['lonCell'][:]) * 180.0 / np.pi lons = np.where(lons0 > 180.0, lons0 - 360.0, lons0) -# Now read the var you want +# Now read the var you want bg_all = [] for imem in range(1, nmems + 1): - infile = 'data/ens/mem%s/%s' % (str(imem).zfill(3), ens_file) - nc = Dataset(infile, 'r') + try: + infile = 'data/ens/mem%s/%s' % (str(imem).zfill(3), ens_file) + nc = Dataset(infile, 'r') + except FileNotFoundError: + infile = 'data/ens/mem%s/%s' % (str(imem).zfill(2), ens_file) + nc = Dataset(infile, 'r') if variable == 'airTemperature': bg = nc.variables['theta'][0, :, lev].astype(np.float64) pres = (nc.variables['pressure_p'][0, :, lev] + nc.variables['pressure_base'][0, :, lev]) / 100.0