diff --git a/.gitignore b/.gitignore index 78f64563e5784a02a66632b07f387220c7a9c16d..6c5fd0777cd0321f61523b5e77e2fb55c3fecff8 100644 --- a/.gitignore +++ b/.gitignore @@ -24,6 +24,7 @@ wheels/ .installed.cfg *.egg MANIFEST +outputs # PyInstaller # Usually these files are written by a python script from a template diff --git a/inputs/namelist_remote.bdy b/inputs/namelist_remote.bdy index ef6bf5d76feac030beb46e965938f45f416800d7..72b58c848840f7985bcd296a2fc4f6b475442701 100644 --- a/inputs/namelist_remote.bdy +++ b/inputs/namelist_remote.bdy @@ -40,8 +40,8 @@ !------------------------------------------------------------------------------ ! I/O !------------------------------------------------------------------------------ - sn_src_dir = './PyNEMO/inputs/src_data_remote.ncml' ! src_files/' - sn_dst_dir = './outputs' + sn_src_dir = '/Users/thopri/Projects/PyNEMO/inputs/src_data_remote.ncml' ! src_files/' + sn_dst_dir = '/Users/thopri/Projects/PyNEMO/outputs' sn_fn = 'NNA_R12' ! prefix for output files nn_fv = -1e20 ! set fill value for output files nn_src_time_adj = 0 ! src time adjustment diff --git a/pynemo/tests/bdy_coords.py b/pynemo/tests/bdy_coords.py new file mode 100755 index 0000000000000000000000000000000000000000..5f801fffa9bbd0b20850bc610c951db156828237 --- /dev/null +++ b/pynemo/tests/bdy_coords.py @@ -0,0 +1,157 @@ +''' +### +## This is a subset of pynemo for Robinson. It should be able to read a mask +## file and produce a coordinates.bdy.nc file (jdha@noc.ac.uk) +### +''' + +# pylint: disable=E1103 +# pylint: disable=no-name-in-module + +#External imports +from time import clock +import numpy as np +import logging + +#local imports +from pynemo import nemo_bdy_setup as setup +from pynemo import nemo_bdy_gen_c as gen_grid +from pynemo import nemo_coord_gen_pop as coord + +from pynemo import nemo_bdy_source_coord as source_coord +from pynemo import nemo_bdy_dst_coord as dst_coord +from pynemo import nemo_bdy_ice +from pynemo import pynemo_settings_editor +from pynemo.utils import Constants + +from pynemo.gui.nemo_bdy_mask import Mask as Mask_File +from PyQt4.QtGui import QMessageBox +#import pickle + + +logger = logging.getLogger(__name__) +def process_bdy(setup_filepath=0, mask_gui=False): + """ Main entry to the processing of the bdy + Keyword arguments: + setup_filepath -- file path to bdy file + mask_gui -- whether gui to select the mask file needs to be poped up + """ + #Logger + logger.info('START') + start = clock() + SourceCoord = source_coord.SourceCoord() + DstCoord = dst_coord.DstCoord() + + logger.info(clock() - start) + start = clock() + Setup = setup.Setup(setup_filepath) # default settings file + settings = Setup.settings + + logger.info(clock() - start) + ice = settings['ice'] + logger.info('ice = %s', ice) + + logger.info('Done Setup') + # default file, region settingas + + start = clock() + bdy_msk = _get_mask(Setup, mask_gui) + logger.info(clock() - start) + logger.info('Done Mask') + + DstCoord.bdy_msk = bdy_msk == 1 + reload(gen_grid) + start = clock() + logger.info('start bdy_t') + grid_t = gen_grid.Boundary(bdy_msk, settings, 't') + + logger.info(clock() - start) + start = clock() + logger.info('start bdy_u') + grid_u = gen_grid.Boundary(bdy_msk, settings, 'u') + logger.info('start bdy_v') + + logger.info(clock() - start) + start = clock() + grid_v = gen_grid.Boundary(bdy_msk, settings, 'v') + logger.info('start bdy_f') + + logger.info(clock() - start) + start = clock() + + grid_f = gen_grid.Boundary(bdy_msk, settings, 'f') + logger.info('done bdy t,u,v,f') + + logger.info(clock() - start) + start = clock() + if ice: + grid_ice = nemo_bdy_ice.BoundaryIce() + grid_ice.grid_type = 't' + grid_ice.bdy_r = grid_t.bdy_r + + bdy_ind = {'t': grid_t, 'u': grid_u, 'v': grid_v, 'f': grid_f} + + for k in bdy_ind.keys(): + logger.info('bdy_ind %s %s %s', k, bdy_ind[k].bdy_i.shape, bdy_ind[k].bdy_r.shape) + + start = clock() + co_set = coord.Coord(settings['dst_dir']+'/coordinates.bdy.nc', bdy_ind) + logger.info('done coord gen') + logger.info(clock() - start) + start = clock() + logger.info(settings['dst_hgr']) + co_set.populate(settings['dst_hgr']) + logger.info('done coord pop') + + logger.info(clock() - start) + + # may need to rethink grid info + # tracer 3d frs over rw + # tracer 2d frs over rw (i.e. ice) + # dyn 2d over 1st rim of T grid (i.e. ssh) + # dyn 2d over 1st rim + # dyn 2d frs over rw + # dyn 3d over 1st rim + # dyn 3d frs over rw + +def _get_mask(Setup, mask_gui): + """ This method reads the mask information from the netcdf file or opens a gui + to create a mask depending on the mask_gui input. return the mask data. The default mask + data is using bathymetry and applying a 1px halo + Keyword arguments: + Setup -- settings for bdy + mask_gui -- boolean to open mask gui. + """ + bdy_msk = None + if mask_gui: + #Open the gui to create a mask + _, mask = pynemo_settings_editor.open_settings_dialog(Setup) + bdy_msk = mask.data + Setup.refresh() + else: + try: + #mask filename and mask file flag is set + if Setup.bool_settings['mask_file'] and Setup.settings['mask_file'] is not None: + mask = Mask_File(mask_file=Setup.settings['mask_file']) + bdy_msk = mask.data + elif Setup.bool_settings['mask_file']: + logger.error("Mask file is not given") + return + else: #no mask file specified then use default 1px halo mask + logger.warning("Using default mask with bathymetry!!!!") + mask = Mask_File(Setup.settings['bathy']) + mask.apply_border_mask(Constants.DEFAULT_MASK_PIXELS) + bdy_msk = mask.data + except ValueError: # why is this except here? as there is an else: statement TODO + print 'something wrong?' + return + if np.amin(bdy_msk) == 0: + # Mask is not set throw a warning message and set border to 1px. + logger.warning("Setting the mask to 1px border") + QMessageBox.warning(None,"pyNEMO", "Mask is not set, setting a 1 pixel border mask") + if bdy_msk is not None and 1 < bdy_msk.shape[0] and 1 < bdy_msk.shape[1]: + tmp = np.ones(bdy_msk.shape, dtype=bool) + tmp[1:-1, 1:-1] = False + bdy_msk[tmp] = -1 + + return bdy_msk diff --git a/test_scripts/bdy_coords_plot.py b/test_scripts/bdy_coords_plot.py new file mode 100755 index 0000000000000000000000000000000000000000..09bbf7c42b5408902b368c6955a4e068dbc4739e --- /dev/null +++ b/test_scripts/bdy_coords_plot.py @@ -0,0 +1,61 @@ +#!/usr/bin/env python2 +# -*- coding: utf-8 -*- +""" +Created on Thu Oct 24 11:28:04 2019 + +@author: thopri +""" + +# Must be run from the PyNEMO main folder. Otherwise permission errors abound + +import matplotlib.pyplot as plt +from netCDF4 import Dataset +import numpy as np +from mpl_toolkits.basemap import Basemap + +from pynemo.tests import bdy_coords as bdc + +bdc.process_bdy('/Users/thopri/Projects/PyNEMO/inputs/namelist_remote.bdy',False) + +rootgrp = Dataset('/Users/thopri/Projects/PyNEMO/outputs/NNA_R12_bdyT_y1979m11.nc', "r", format="NETCDF4") +bdy_msk = np.squeeze(rootgrp.variables['bdy_msk'][:]) - 1 +bdy_lon = np.squeeze(rootgrp.variables['nav_lon'][:]) +bdy_lat = np.squeeze(rootgrp.variables['nav_lat'][:]) +rootgrp.close() + +rootgrp = Dataset('/Users/thopri/Projects/PyNEMO/outputs/coordinates.bdy.nc', "r", format="NETCDF4") +bdy_it = np.squeeze(rootgrp.variables['nbit'][:]) +bdy_jt = np.squeeze(rootgrp.variables['nbjt'][:]) +bdy_rt = np.squeeze(rootgrp.variables['nbrt'][:]) +rootgrp.close() + +bdy_msk = np.ma.masked_where(bdy_msk < 0, bdy_msk) +for f in range(len(bdy_it)): + bdy_msk[bdy_jt[f,],bdy_it[f,]] = bdy_rt[f,] + +# Plot to check output +fig=plt.figure(figsize=(12,12)) +ax=fig.add_subplot(111) + +map = Basemap(llcrnrlat=41.,urcrnrlat=65.,\ + llcrnrlon=-22.,urcrnrlon=25.,\ + rsphere=(6378137.00,6356752.3142),\ + resolution='l',projection='lcc',\ + lat_1=57.,lon_0=-12.5) +map.drawcoastlines() +map.drawcountries() +map.fillcontinents(color='grey') +map.drawmeridians(np.arange(-25.,25.,2),labels=[0,0,0,1]) +map.drawparallels(np.arange(40.,66.,2),labels=[1,0,0,0]) + +cmap = plt.cm.get_cmap('jet',10) +cmaplist = [cmap(i) for i in range(cmap.N)] +cmaplist[0] = (.7,.7,.7,1.0) +cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N) + +x1,y1 = map(bdy_lon[:,:],bdy_lat[:,:]) +im = map.pcolormesh(x1, y1, bdy_msk, cmap=cmap, vmin=-0.5, vmax=9.5) +cb = plt.colorbar(orientation='vertical', shrink=0.75, aspect=30, fraction=0.1,pad=0.05) +cb.set_label('RimWidth Number') +cb.set_ticks(np.arange(10)) +th=plt.title(('BDY Points'),fontweight='bold') \ No newline at end of file