From 6614403b3f45b2f2a532ee3ae00f614b9e266701 Mon Sep 17 00:00:00 2001 From: Micael Oliveira Date: Wed, 7 Sep 2022 15:44:35 +1000 Subject: [PATCH] Add two new tools: one automatically fixes non-advective cells found in a mask file; the other combines two masks. Also modified the topog2mask.py utility to take treat cells with a fraction of ocean lower than a user-provided value as land. --- combine_masks.py | 60 +++++++++++++++++++++++++++++++++ fix_nonadvective_mask.py | 72 ++++++++++++++++++++++++++++++++++++++++ topog2mask.py | 19 ++++++----- 3 files changed, 143 insertions(+), 8 deletions(-) create mode 100755 combine_masks.py create mode 100755 fix_nonadvective_mask.py diff --git a/combine_masks.py b/combine_masks.py new file mode 100755 index 0000000..8d0e319 --- /dev/null +++ b/combine_masks.py @@ -0,0 +1,60 @@ +#!/usr/bin/env python + +import sys +import os +import numpy as np +import argparse +import netCDF4 as nc + +def fix_mask(mask1_filename, mask2_filename, grid_filename, lat, output_mask_filename): + + lat = -60 # we will fix the mask south of this latitude + + with nc.Dataset(mask1_filename, 'r') as mf1, \ + nc.Dataset(mask2_filename, 'r') as mf2, \ + nc.Dataset(grid_filename, 'r') as gr, \ + nc.Dataset(output_mask_filename, 'w') as mf_out: + + # Sanity checks + num_lons = mf1.dimensions['nx'].size + num_lats = mf1.dimensions['ny'].size + if 2*num_lons != gr.dimensions['nx'].size or 2*num_lats != gr.dimensions['ny'].size: + print('Error: dimensions of {} and {} are not compatible'.format(mask1_filename, grid_filename), file=sys.stderr) + return 1 + + if num_lons != mf2.dimensions['nx'].size or num_lats != mf2.dimensions['ny'].size: + print('Error: dimensions of {} and {} are not compatible'.format(mask1_filename, mask2_filename), file=sys.stderr) + return 1 + + # Get the coordinates of the T points from the supergrid: + yt = gr.variables['y'][1::2,1::2] + + # Combine the two masks, taking mask 1 north of 'lat' and mask 2 south of 'lat'. + combined = np.where(yt > lat, mf1.variables['mask'][:], mf2.variables['mask'][:]) + + # Output mask + mf_out.createDimension('nx', num_lons) + mf_out.createDimension('ny', num_lats) + mask_out = mf_out.createVariable('mask', 'f8', dimensions=('ny', 'nx')) + mask_out[:] = combined[:] + +def main(): + + parser = argparse.ArgumentParser() + parser.add_argument('mask1', help='First mask to combine.') + parser.add_argument('mask2', help='Second mask to combine.') + parser.add_argument('grid', help='Ocean super-grid.') + parser.add_argument('output_mask', help='The combined mask.') + parser.add_argument('latitude', help='Take mask1 north of this value and mask2 otherwise.') + + args = parser.parse_args() + + if not os.path.exists(args.mask1) or not os.path.exists(args.mask2) or not os.path.exists(args.grid): + print('Error: one of {} {} was not found'.format(args.topog, args.input_mask, args.grid), file=sys.stderr) + parser.print_help() + return 1 + + fix_mask(args.mask1, args.mask2, args.grid, args.latitude, args.output_mask) + +if __name__ == '__main__': + sys.exit(main()) diff --git a/fix_nonadvective_mask.py b/fix_nonadvective_mask.py new file mode 100755 index 0000000..086a6e9 --- /dev/null +++ b/fix_nonadvective_mask.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python + +import sys +import os +import numpy as np +import argparse +import netCDF4 as nc + +def fix_mask(input_mask_filename, output_mask_filename): + + with nc.Dataset(input_mask_filename, 'r') as mf_in, \ + nc.Dataset(output_mask_filename, 'w') as mf_out: + + mask = mf_in.variables['mask'][:,:] + num_lons = mf_in.dimensions['nx'].size + num_lats = mf_in.dimensions['ny'].size + print('mask dimensions {:11d} {:11d}'.format(num_lons, num_lats)) + + # Create mask with halo (includes copied points along LH edge and RH edge and points along tripole seam) + mask_halo = np.zeros((num_lats+1,num_lons+2), dtype=int) + mask_halo[:-1,1:-1] = mask[:,:] + + mask_halo[:,-1] = mask_halo[:,1] + mask_halo[:,0] = mask_halo[:,-2] + mask_halo[-1, :] = mask_halo[-2, ::-1] + + # Find non-advective cells + done = False + n_total = 0 + print('non-advective cells:') + while (not done): + sw = (mask_halo[1:-1, :-2] == 0) | (mask_halo[ :-2, :-2] == 0) | (mask_halo[ :-2, 1:-1] == 0) + se = (mask_halo[ :-2, 1:-1] == 0) | (mask_halo[ :-2, 2: ] == 0) | (mask_halo[1:-1, 2: ] == 0) + ne = (mask_halo[1:-1, 2: ] == 0) | (mask_halo[2: , 2: ] == 0) | (mask_halo[2: , 1:-1] == 0) + nw = (mask_halo[2: , 1:-1] == 0) | (mask_halo[2: , :-2] == 0) | (mask_halo[1:-1, :-2] == 0) + advective_cells = sw & se & ne & nw & (mask_halo[1:-1,1:-1] == 1) + + jj,ii = np.nonzero(advective_cells) + for ij in list(zip(ii,jj+1)): + print('{:11d} {:11d}'.format(ij[0], ij[1])) + + mask_halo[jj+1,ii+1] = 0 + + n_new = np.count_nonzero(advective_cells) + done = n_new == 0 + n_total += n_new + + print('Found {} non-advective cells'.format(n_total)) + + mf_out.createDimension('nx', num_lons) + mf_out.createDimension('ny', num_lats) + + mask_out = mf_out.createVariable('mask', 'f8', dimensions=('ny', 'nx')) + mask_out[:,:] = mask_halo[:-1,1:-1] + +def main(): + + parser = argparse.ArgumentParser() + parser.add_argument('input_mask', help='Mask to fix.') + parser.add_argument('output_mask', help='Output mask.') + + args = parser.parse_args() + + if not os.path.exists(args.input_mask): + print('Error: {} was not found'.format(args.input_mask), file=sys.stderr) + parser.print_help() + return 1 + + fix_mask(args.input_mask, args.output_mask) + +if __name__ == '__main__': + sys.exit(main()) diff --git a/topog2mask.py b/topog2mask.py index 9e3bd83..23661ab 100755 --- a/topog2mask.py +++ b/topog2mask.py @@ -8,7 +8,7 @@ import argparse import netCDF4 as nc -def make_mask(topog_filename, mask_filename, ice=False): +def make_mask(topog_filename, mask_filename, frac, ice=False): with nc.Dataset(topog_filename, 'r') as tf, \ nc.Dataset(mask_filename, 'w') as mf: @@ -26,10 +26,7 @@ def make_mask(topog_filename, mask_filename, ice=False): else: mask = mf.createVariable('mask', 'f8', dimensions=('ny', 'nx')) # CICE and MOM use 0 as masked - m = np.zeros_like(mask) - m[:] = 1.0 - m[np.where(tf.variables['depth'][:] <= 0.0)] = 0.0 - mask[:] = m[:] + mask[:] = np.where((tf.variables['frac'][:] < frac) | (tf.variables['depth'][:] <= 0.0), 0, 1) def main(): @@ -38,7 +35,7 @@ def main(): parser.add_argument('topog', help='The topog file.') parser.add_argument('cice_mask', help='The new CICE mask file.') parser.add_argument('mom_mask', help='The new MOM mask file.') - + parser.add_argument('fraction', type=float, nargs='?', default=0.0, help='Cells with a fraction of ocean lower than this value will be treated as land') args = parser.parse_args() if not os.path.exists(args.topog): @@ -53,8 +50,14 @@ def main(): parser.print_help() return 1 - make_mask(args.topog, args.cice_mask, ice=True) - make_mask(args.topog, args.mom_mask) + if args.fraction > 1.0 or args.fraction < 0.0: + print('Error: fraction must be between 0 and 1, but it is {}'.format(args.fraction), + file=sys.stderr) + parser.print_help() + return 1 + + make_mask(args.topog, args.cice_mask, args.fraction, ice=True) + make_mask(args.topog, args.mom_mask, args.fraction) if __name__ == '__main__': sys.exit(main())