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..80c6992 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,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,17 +99,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 +112,14 @@ 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) @@ -143,16 +127,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..292d092 --- 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,31 +31,32 @@ # 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) +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) ################################################################################### # 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' 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] @@ -68,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) @@ -80,7 +78,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 +96,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 +121,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 +149,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 +162,14 @@ 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) +# 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()) # Scatter the single ob location m1.scatter(singleob_lon, singleob_lat, color='g', marker='o', s=2) @@ -192,16 +180,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..36f86fc 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 warnings.filterwarnings('ignore') ############ USER INPUT ########################################################## # MPAS info -jstatic = "data/static.nc" # to load the MPAS lat/lon +jstatic = "data/static.nc" # to load the MPAS lat/lon ens_file = "restart.2024-05-27_00.00.00.nc" +#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 +# Now read the var you want bg_all = [] -for imem in range(1, nmems+1): +for imem in range(1, nmems + 1): + 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 - dividend = (1000.0/pres)**(0.286) - bg = bg / dividend - else: - sys.exit("Error: requested variable not yet implemented.") - bg_all.append(bg) + 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)) +# 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) +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() +