Merge pull request #39 from NOC-MSM/master

Bringing tide branch up to date

.gitignore

@@ -5,7 +5,8 @@ __pycache__/
 
 # C extensions
 *.so
-
+*.nc
+*.xml
 # Distribution / packaging
 .Python
 build/
@@ -14,6 +15,7 @@ dist/
 downloads/
 eggs/
 .eggs/
+.idea/
 lib/
 lib64/
 parts/
@@ -24,6 +26,8 @@ wheels/
 .installed.cfg
 *.egg
 MANIFEST
+outputs
+.DS_Store
 
 # PyInstaller
 #  Usually these files are written by a python script from a template

README.md

-# PyNEMO
- 
-To be udated soon. This work springboards from the [PyNEMO](http://pynemo.readthedocs.io/en/latest/index.html) Project.
-
-## What is this repository for? ##
-
-## How do I get set up? ##
-
-## Contribution guidelines ##
-
-## Who do I talk to? ##
-
-* Repo owner or admin
-
-jdha
-
-* Other community or team contact
-
-
-For more information regarding the use and development of PyNEMO see: [PyNEMO Wiki](https://github.com/jdha/PyNEMO/wiki)

README.rst

+PyNEMO
+======
+
+To be updated soon. This work springboards from the `PyNEMO Project <http://pynemo.readthedocs.io/en/latest/index.html/>`_.
+
+What is this repository for?
+----------------------------
+
+How do I get set up?
+--------------------
+
+Steps to take to install PyNEMO, creating a specific conda virtual environment is highly recommended. 
+`click here for more about virtual enviroments <https://docs.conda.io/projects/conda/en/latest/user-guide/tasks/manage-environments.html/>`_
+
+- Install Git (outside scope of this readme)
+- Clone PyNEMO repository::
+    
+    $ git clone https://github.com/NOC-MSM/PyNEMO.git 
+    
+- Install Conda, either Anaconda or Miniconda (outside scope of this readme)
+- Create conda environment for PyNEMO::
+
+    $ cd to/PyNEMO/directory
+    $ conda env create -f environment_pynemo.yml
+
+- Activate the new virtual environment::
+
+   $ source activate pynemo_env
+
+- Install Jave JRE (outside scope of this readme) and link libjvm.dylib to LD_LIBRARY_PATH variable::
+
+    $ export LD_LIBRARY_PATH=/path/to/java/library/folder/containing/libjvm.dylib:$LD_LIBARY_PATH # see notes below
+
+- Install PyNEMO::
+  
+    $ cd /location/of/pynemo/repo 
+    $ python setup.py build
+    $ python setup.py install
+
+This should result in PyNEMO being installed in the virtual environment, and can be checked by entering::  
+
+    $ pynemo -v
+
+Resulting in a help usage prompt::
+ 
+    $ usage: pynemo -g -s <namelist.bdy> 
+
+The virtual environment can be deactivated to return you to the normal prompt by typing::  
+    
+$ conda deactivate
+To reactivate, the following needs to be typed::
+
+    $ source activate pynemo_env
+
+
+To use PyNEMO, the following command is entered: (the example will run an benchmarking test)::
+
+    $ pynemo -s /path/to/namelist/file (e.g. PyNEMO/inputs/namelist_remote.bdy)
+
+**Additional NOTES** 
+
+For Macbook Pro 2015, macOS Mojave and Java SDK 13 and JRE 8 the following path for the libjvm library should be correct:: 
+
+    /Library/Java/JavaVirtualMachines/jdk-13.0.1.jdk/Contents/Home/lib/server
+
+Resulting in the following command: (this will be different for different java versions and operating systems)::
+
+    $ export LD_LIBRARY_PATH=/Library/Java/JavaVirtualMachines/jdk-13.0.1.jdk/Contents/Home/lib/server:$LD_LIBRARY_PATH
+
+For an iMac 2013, macOS Catalina and JRE 8 only the followinng path was found to be correct::
+    
+   /Library/Internet\ Plug-Ins/JavaAppletPlugin.plugin/Contents/Home/lib/server
+   
+With the following command being required to set the environment variable::
+
+    $ export LD_LIBRARY_PATH=/Library/Internet\ Plug-Ins/JavaAppletPlugin.plugin/Contents/Home/lib/server:$LD_LIBRARY_PATH
+
+The conda environment creation command has not yet been tested. The yml document (can be opened using text editor) gives a list of all the modules and their versions that are required for PyNEMO so a environment can be constructed using this document as reference (or if you use pip!)
+
+**Update** conda environment yaml file has been tested (and works!) on a Macbook Pro 2015 and iMac 2013 running Anaconda 3.7 and Miniconda 3.7 respectively. 
+
+Contribution guidelines
+-----------------------
+
+Bench Marking Tests
+-------------------
+
+The PyNEMO module can be tested using the bench marking namelist bdy file in the inputs folder. To check the outputs of the benchmark test, these can be visualised using the plotting script within the test_scripts folder. The following steps are required,
+
+- Run PyNEMO using the namelist file in the inputs folder (namelist_remote.bdy) e.g.::
+
+    $ pynemo -s /path/to/namelist/file
+
+- This will create two output files coordinates.bdy.nc and NNA_R12_bdyT_y1979)m11.nc in an outputs folder
+
+- To check the coordinates.bdy.nc has the correct boundary points, the script bdy_coords_plot.py will plot the domain boundaries and shown the different locations of the rim width (increasing number should go inwards) This script is located in the test_scripts folder.
+
+- The result should look like this (if using the current benchmark data)
+
+.. image:: /screenshots/example_bdy_coords.png
+  :width: 800
+  :alt: Example BDY coords output
+
+Who do I talk to?
+-----------------
+
+* Repo owner or admin
+
+  jdha
+
+* Other community or team contact
+
+
+For more information regarding the use and development of PyNEMO see: [PyNEMO Wiki](https://github.com/jdha/PyNEMO/wiki)

README_markdown.md

+# PyNEMO
+ 
+To be udated soon. This work springboards from the [PyNEMO](http://pynemo.readthedocs.io/en/latest/index.html) Project.
+
+## What is this repository for? ##
+
+## How do I get set up? ##
+
+Steps to take to install PyNEMO, creating a specific conda virtual environment is highly recommended. [click here for more about virtual enviroments](https://docs.conda.io/projects/conda/en/latest/user-guide/tasks/manage-environments.html)
+
+- Install Conda, either Anaconda or Miniconda (outside of this readme)
+- Create conda environment for PyNEMO
+```
+$ cd to/PyNEMO/directory
+```
+```
+$ conda env create -f environment_pynemo.yml
+```
+- Activate the new virtual environment
+```
+$ source activate pynemo_env
+```
+- Install Jave JRE (outside this readme) and link libjvm.dylib to LD_LIBRARY_PATH variable
+```
+$ export LD_LIBRARY_PATH=/path/to/java/library/folder/containing/libjvm.dylib:$LD_LIBARY_PATH # see notes below
+```
+- Install Git (outside this readme)
+```
+$ git clone https://github.com/NOC-MSM/PyNEMO.git
+```
+- Install PyNEMO:  
+```
+$ cd /location/of/pynemo/repo 
+```
+```
+$ python setup.py build
+```
+```
+$ python setup.py install
+```
+
+This should result in PyNEMO being installed in the virtual environment, and can be checked by entering:  
+```
+$ pynemo -v
+```
+
+Resulting in a help usage prompt:
+``` 
+$ usage: pynemo -g -s <namelist.bdy> 
+```
+
+The virtual environment can be deactivated to return you to the normal prompt by typing:  
+```
+$ conda deactivate
+```
+
+To reactivate, the following needs to be typed
+```
+$ source activate pynemo_env
+```
+
+To use PyNEMO, the following command is entered: (the example will run an benchmarking test)
+```
+$ pynemo -s /path/to/namelist/file (e.g. PyNEMO/inputs/namelist_remote.bdy)
+```
+
+**Additional NOTES** 
+
+for MacOs and Java SDK 13 and JRE 8 the following path should be correct - /Library/Java/JavaVirtualMachines/jdk-13.0.1.jdk/Contents/Home/lib/server
+
+Resulting in the following command: (this will be different for different java versions and operating systems)
+```
+$ export LD_LIBRARY_PATH=/Library/Java/JavaVirtualMachines/jdk-13.0.1.jdk/Contents/Home/lib/server:$LD_LIBRARY_PATH
+```
+The conda environment creation command has not yet been tested. The yml document (can be opened using text editor) gives a list of all the modules and their versions that are required for PyNEMO so a environment can be constructed using this document as reference (or if you use pip!)
+
+## Contribution guidelines ##
+
+## Bench Marking Tests ##
+
+The PyNEMO module can be tested using the bench marking namelist bdy file in the inputs folder. To check the outputs of the benchmark test, these can be visualised using the plotting script within the test_scripts folder. The following steps are required,
+
+- Run PyNEMO using the namelist file in the inputs folder (namelist_remote.bdy) e.g.
+
+    -   $ pynemo -s /path/to/namelist/file
+
+- This will create two output files coordinates.bdy.nc and NNA_R12_bdyT_y1979)m11.nc in an outputs folder
+
+- To check the coordinates.bdy.nc has the correct boundary points, the script bdy_coords_plot.py will plot the domain boundaries and shown the different locations of the rim width (increasing number should go inwards) This script is located in the test_scripts folder.
+
+- The result should look like this (if using the current benchmark data)
+
+![Example BDY coords output](/screenshots/example_bdy_coords.png)
+
+## Who do I talk to? ##
+
+* Repo owner or admin
+
+jdha
+
+* Other community or team contact
+
+
+For more information regarding the use and development of PyNEMO see: [PyNEMO Wiki](https://github.com/jdha/PyNEMO/wiki)

environment_pynemo.yml

+name: pynemo_env
+channels:
+  - conda-forge
+  - anaconda
+  - srikanthnagella
+  - defaults
+dependencies:
+  - alabaster=0.7.12=py27_0
+  - appnope=0.1.0=py27hb466136_0
+  - asn1crypto=1.0.1=py27_0
+  - astroid=1.6.5=py27_0
+  - attrs=19.2.0=py_0
+  - babel=2.7.0=py_0
+  - backports=1.0=py_2
+  - backports.functools_lru_cache=1.5=py_2
+  - backports.shutil_get_terminal_size=1.0.0=py27_2
+  - backports_abc=0.5=py27h6972548_0
+  - basemap=1.2.0=py27h0acbc05_0
+  - blas=1.0=mkl
+  - bleach=3.1.0=py27_0
+  - bzip2=1.0.8=h1de35cc_0
+  - ca-certificates=2019.10.16=0
+  - cartopy=0.16.0=py27h9263bd1_0
+  - certifi=2019.9.11=py27_0
+  - cffi=1.12.3=py27hb5b8e2f_0
+  - cftime=1.0.3.4=py27h1d22016_1001
+  - chardet=3.0.4=py27_1003
+  - cloudpickle=1.2.2=py_0
+  - configparser=4.0.2=py27_0
+  - cryptography=2.7=py27ha12b0ac_0
+  - curl=7.65.3=ha441bb4_0
+  - cycler=0.10.0=py27hfc73c78_0
+  - cython=0.29.13=py27h0a44026_0
+  - dbus=1.13.6=h90a0687_0
+  - decorator=4.4.0=py27_1
+  - defusedxml=0.6.0=py_0
+  - docutils=0.15.2=py27_0
+  - entrypoints=0.3=py27_0
+  - enum34=1.1.6=py27_1
+  - expat=2.2.6=h0a44026_0
+  - freetype=2.9.1=hb4e5f40_0
+  - functools32=3.2.3.2=py27_1
+  - futures=3.3.0=py27_0
+  - geos=3.6.2=h5470d99_2
+  - gettext=0.19.8.1=h15daf44_3
+  - glib=2.56.2=hd9629dc_0
+  - hdf4=4.2.13=h39711bb_2
+  - hdf5=1.10.1=ha036c08_1
+  - icu=58.2=h4b95b61_1
+  - idna=2.8=py27_0
+  - imagesize=1.1.0=py27_0
+  - intel-openmp=2019.4=233
+  - ipaddress=1.0.22=py27_0
+  - ipykernel=4.10.0=py27_0
+  - ipython=5.8.0=py27_0
+  - ipython_genutils=0.2.0=py27h8b9a179_0
+  - ipywidgets=7.5.1=py_0
+  - isort=4.3.21=py27_0
+  - jedi=0.15.1=py27_0
+  - jinja2=2.10.3=py_0
+  - jpeg=9b=he5867d9_2
+  - jsonschema=3.0.2=py27_0
+  - jupyter=1.0.0=py27_7
+  - jupyter_client=5.3.3=py27_1
+  - jupyter_console=5.2.0=py27_1
+  - jupyter_core=4.5.0=py_0
+  - keyring=18.0.0=py27_0
+  - kiwisolver=1.1.0=py27h0a44026_0
+  - krb5=1.16.1=hddcf347_7
+  - lazy-object-proxy=1.4.2=py27h1de35cc_0
+  - libcurl=7.65.3=h051b688_0
+  - libcxx=4.0.1=hcfea43d_1
+  - libcxxabi=4.0.1=hcfea43d_1
+  - libedit=3.1.20181209=hb402a30_0
+  - libffi=3.2.1=h475c297_4
+  - libgfortran=3.0.1=h93005f0_2
+  - libiconv=1.15=hdd342a3_7
+  - libnetcdf=4.5.0=h42fd751_7
+  - libpng=1.6.37=ha441bb4_0
+  - libsodium=1.0.16=h3efe00b_0
+  - libssh2=1.8.2=ha12b0ac_0
+  - libtiff=4.0.10=hcb84e12_2
+  - libxml2=2.9.9=hf6e021a_1
+  - libxslt=1.1.33=h33a18ac_0
+  - lxml=4.4.1=py27hef8c89e_0
+  - markupsafe=1.1.1=py27h1de35cc_0
+  - matplotlib=2.2.3=py27h54f8f79_0
+  - mccabe=0.6.1=py27_1
+  - mistune=0.8.4=py27h1de35cc_0
+  - mkl=2019.4=233
+  - mkl-service=2.3.0=py27hfbe908c_0
+  - mkl_fft=1.0.14=py27h5e564d8_0
+  - mkl_random=1.1.0=py27ha771720_0
+  - nbconvert=5.6.0=py27_1
+  - nbformat=4.4.0=py27hddc86d0_0
+  - ncurses=6.1=h0a44026_1
+  - netcdf4=1.3.1=py27he3ffdca_2
+  - notebook=5.7.8=py27_0
+  - numpy=1.16.5=py27hacdab7b_0
+  - numpy-base=1.16.5=py27h6575580_0
+  - numpydoc=0.9.1=py_0
+  - olefile=0.46=py27_0
+  - openssl=1.1.1d=h1de35cc_3
+  - owslib=0.18.0=py_0
+  - packaging=19.2=py_0
+  - pandoc=2.2.3.2=0
+  - pandocfilters=1.4.2=py27_1
+  - parso=0.5.1=py_0
+  - pathlib2=2.3.5=py27_0
+  - pcre=8.43=h0a44026_0
+  - pexpect=4.7.0=py27_0
+  - pickleshare=0.7.5=py27_0
+  - pillow=6.2.0=py27hb68e598_0
+  - pip=19.2.3=py27_0
+  - proj4=5.0.1=h1de35cc_0
+  - prometheus_client=0.7.1=py_0
+  - prompt_toolkit=1.0.15=py27h4a7b9c2_0
+  - psutil=5.6.3=py27h1de35cc_0
+  - ptyprocess=0.6.0=py27_0
+  - pycodestyle=2.5.0=py27_0
+  - pycparser=2.19=py27_0
+  - pyepsg=0.4.0=py27_0
+  - pyflakes=2.1.1=py27_0
+  - pygments=2.4.2=py_0
+  - pyjnius=1.4=py27_0
+  - pylint=1.9.2=py27_0
+  - pyopenssl=19.0.0=py27_0
+  - pyparsing=2.4.2=py_0
+  - pyproj=1.9.5.1=py27h833a5d7_1
+  - pyqt=4.11.4=py27_4
+  - pyrsistent=0.15.4=py27h1de35cc_0
+  - pyshp=2.1.0=py_0
+  - pysocks=1.7.1=py27_0
+  - python=2.7.16=h97142e2_7
+  - python-dateutil=2.8.0=py27_0
+  - python.app=2=py27_9
+  - pytz=2019.3=py_0
+  - pyzmq=18.1.0=py27h0a44026_0
+  - qt=4.8.7=1
+  - qtawesome=0.6.0=py_0
+  - qtconsole=4.5.5=py_0
+  - qtpy=1.9.0=py_0
+  - readline=7.0=h1de35cc_5
+  - requests=2.22.0=py27_0
+  - rope=0.14.0=py_0
+  - scandir=1.10.0=py27h1de35cc_0
+  - scipy=1.2.1=py27h1410ff5_0
+  - seawater=3.3.4=py_1
+  - send2trash=1.5.0=py27_0
+  - setuptools=41.4.0=py27_0
+  - shapely=1.6.4=py27h20de77a_0
+  - simplegeneric=0.8.1=py27_2
+  - singledispatch=3.4.0.3=py27he22c18d_0
+  - sip=4.18=py27_0
+  - six=1.12.0=py27_0
+  - snowballstemmer=2.0.0=py_0
+  - sphinx=1.8.5=py27_0
+  - sphinxcontrib=1.0=py27_1
+  - sphinxcontrib-websupport=1.1.2=py_0
+  - spyder=3.2.8=py27_0
+  - spyder-kernels=1.4.0=py27_0
+  - sqlite=3.30.0=ha441bb4_0
+  - subprocess32=3.5.4=py27h1de35cc_0
+  - terminado=0.8.2=py27_0
+  - testpath=0.4.2=py27_0
+  - thredds_crawler=1.0.0=py27_0
+  - tk=8.6.8=ha441bb4_0
+  - tornado=5.1.1=py27h1de35cc_0
+  - traitlets=4.3.3=py27_0
+  - typing=3.7.4.1=py27_0
+  - urllib3=1.24.2=py27_0
+  - wcwidth=0.1.7=py27h817c265_0
+  - webencodings=0.5.1=py27_1
+  - wheel=0.33.6=py27_0
+  - widgetsnbextension=3.5.1=py27_0
+  - wrapt=1.11.2=py27h1de35cc_0
+  - wurlitzer=1.0.3=py27_0
+  - xz=5.2.4=h1de35cc_4
+  - zeromq=4.3.1=h0a44026_3
+  - zlib=1.2.11=h1de35cc_3
+  - zstd=1.3.7=h5bba6e5_0

inputs/namelist_remote.bdy

@@ -30,18 +30,18 @@
 !------------------------------------------------------------------------------
 !  grid information 
 !------------------------------------------------------------------------------
-   sn_src_hgr = './benchmark/grid_low_res_C/mesh_hgr.nc' 
-   sn_src_zgr = './benchmark/grid_low_res_C/mesh_zgr.nc'
-   sn_dst_hgr = './benchmark/grid_C/mesh_hgr_zps.nc'
-   sn_dst_zgr = './benchmark/grid_C/mesh_zgr_zps.nc'
-   sn_src_msk = './benchmark/grid_low_res_C/mask.nc'
-   sn_bathy   = './benchmark/grid_C/NNA_R12_bathy_meter_bench.nc'
+   sn_src_hgr = 'http://opendap4gws.jasmin.ac.uk/thredds/noc_msm/dodsC/pynemo_grid_low_res_C/mesh_hgr.nc' 
+   sn_src_zgr = 'http://opendap4gws.jasmin.ac.uk/thredds/noc_msm/dodsC/pynemo_grid_low_res_C/mesh_zgr.nc'
+   sn_dst_hgr = 'http://opendap4gws.jasmin.ac.uk/thredds/noc_msm/dodsC/pynemo_grid_C/mesh_hgr_zps.nc'
+   sn_dst_zgr = 'http://opendap4gws.jasmin.ac.uk/thredds/noc_msm/dodsC/pynemo_grid_C/mesh_zgr_zps.nc'
+   sn_src_msk = 'http://opendap4gws.jasmin.ac.uk/thredds/noc_msm/dodsC/pynemo_grid_low_res_C/mask.nc'
+   sn_bathy   = 'http://opendap4gws.jasmin.ac.uk/thredds/noc_msm/dodsC/pynemo_grid_C/NNA_R12_bathy_meter_bench.nc'
 
 !------------------------------------------------------------------------------
 !  I/O 
 !------------------------------------------------------------------------------
-   sn_src_dir = '/Users/jdha/Projects/GitHub/PyNEMO/inputs/src_data.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

pynemo/nemo_bdy_zgrv2.py

@@ -68,6 +68,7 @@ class Depth:
         v_ind = sub2ind(mbathy.shape, bdy_v[:,0], bdy_v[:,1])
         v_ind2 = sub2ind(mbathy.shape, bdy_v[:,0], bdy_v[:,1] + 1) 
    
+    	[tmp_zt, tmp_zw] = e3_to_depth(np.squeeze(nc['e3t'][:,:,:,:]), np.squeeze(nc['e3w'][:,:,:,:]), nz)
         # This is very slow
         self.logger.debug( 'starting nc reads loop' )
         for k in range(nz):
@@ -85,7 +86,8 @@ class Depth:
 		# jelt: replace 'load gdep[wt] with load e3[tw] and compute gdep[tw]
                 #wrk1 = nc['gdept'][0,k,:,:]#nc.variables['gdept'][0,k,:,:]
                 #wrk2 = nc['gdepw'][0,k,:,:]#nc.variables['gdepw'][0,k,:,:]
-		[wrk1, wrk2] = e3_to_depth(nc['e3t'][0,k,:,:], nc['e3w'][0,k,:,:], nz)
+                #print 'e3t shape: ', nc['e3t_0'][:].shape
+                [wrk1, wrk2] = tmp_zt[k,:,:], tmp_zw[k,:,:] 
 
             # Replace deep levels that are not used with NaN
             wrk2[mbathy + 1 < k + 1] = np.NaN

pynemo/reader/ncml.py

@@ -225,7 +225,7 @@ class Variable(object):
             retval = data_array.copyToNDJavaArray() 
             #TODO: copy it into numpy instead of Java array and then convert to numpy
             # convert to numpy array
-            retval = np.asarray(retval)
+            retval = np.asarray(retval, dtype='float')
             self.logger.info(retval.shape)
             if np_input: #if an array is passed as selection
                 ret_dim_list = ()

pynemo/tests/bdy_coords.py

+'''
+###
+## 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

test_scripts/bdy_coords_plot.py

+#!/usr/bin/env python2
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Oct 24 11:28:04 2019
+
+@author: thopri
+"""
+
+# Hard file/folder paths in namelist file and in process_bdy call below need to be updated to suit user
+
+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')
+
+plt.show()
\ No newline at end of file

test_scripts/bdy_var_plot.py

+#!/usr/bin/env python2
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import scipy.spatial as sp
+from netCDF4 import Dataset
+import matplotlib.pyplot as plt
+from matplotlib.collections import PatchCollection
+from matplotlib.patches import Polygon
+
+
+def nemo_bdy_order(fname):
+    """
+    Determine the ordering and breaks in BDY files to aid plotting.
+    This function takes the i/j coordinates from BDY files and orders them sequentially
+    making it easier to visualise sections along the open boundary. Breaks in the open
+    boundary are also determined (i.e. where the distance between two points > 2**0.5)
+    Args:
+        fname     (str) : filename of BDY file
+    Returns:
+        bdy_ind   (dict): re-ordered indices
+        bdy_dst   (dict): distance (in model coords) between points
+        bdy_brk   (dict): location of the break points in the open boundary
+    """
+
+    # open file pointer and extract data
+    rootgrp = Dataset(fname, "r", format="NETCDF4")
+    nbi = np.squeeze(rootgrp.variables['nbidta'][:, :]) - 1  # subtract 1 for python indexing
+    nbj = np.squeeze(rootgrp.variables['nbjdta'][:, :]) - 1
+    nbr = np.squeeze(rootgrp.variables['nbrdta'][:, :]) - 1
+    lon = np.squeeze(rootgrp.variables['nav_lon'][:, :])
+    lat = np.squeeze(rootgrp.variables['nav_lat'][:, :])
+    rootgrp.close()
+
+    rw = np.amax(nbr) + 1
+    bdy_ind = {}
+    bdy_brk = {}
+    bdy_dst = {}
+    nbdy = []
+
+    for r in range(rw):
+        nbdy.append(np.sum(nbr == r))
+
+    # TODO: deal with domain that spans wrap
+
+    # start with outer rim and work in
+
+    for r in range(rw):
+
+        # set initial constants
+
+        ind = nbr == r
+        nbi_tmp = nbi[ind]
+        nbj_tmp = nbj[ind]
+
+        count = 1
+        id_order = np.zeros((nbdy[r], 1), dtype=int) - 1
+        id_order[0,] = 0
+        flag = False
+        mark = 0
+        source_tree = sp.cKDTree(zip(nbi_tmp, nbj_tmp), balanced_tree=False, compact_nodes=False)
+
+        # order bdy entries
+
+        while count < nbdy[r]:
+
+            nn_dist, nn_id = source_tree.query(zip(nbi_tmp[id_order[count - 1]], nbj_tmp[id_order[count - 1]]),
+                                               k=3, distance_upper_bound=2.9)
+            if np.sum(id_order == nn_id[0, 1]) == 1:  # is the nearest point already in the list?
+                if np.sum(id_order == nn_id[0, 2]) == 1:  # is the 2nd nearest point already in the list?
+                    if flag == False:  # then we've found an end point and we can start the search in earnest!
+
+                        flag = True
+                        id_order[mark] = id_order[count - 1]  # reset values
+                        id_order[mark + 1:] = -1  # reset values
+                        count = mark + 1  # reset counter
+                    else:
+                        i = 0  # should this be zero?
+                        t = count
+                        while count == t:
+                            if np.sum(id_order == i) == 1:
+                                i += 1
+                            else:
+                                id_order[count] = i
+                                flag = False
+                                mark = count
+                                count += 1
+                elif nn_id[0, 2] == nbdy[r]:
+                    i = 0
+                    t = count
+                    while count == t:
+                        if np.sum(id_order == i) == 1:
+                            i += 1
+                        else:
+                            id_order[count] = i
+                            flag = False
+                            mark = count
+                            count += 1
+                else:
+                    id_order[count] = nn_id[0, 2]
+                    count += 1
+            else:
+                id_order[count] = nn_id[0, 1]
+                count += 1
+
+        bdy_ind[r] = id_order
+        bdy_dst[r] = np.sqrt((np.diff(np.hstack((nbi_tmp[id_order], nbj_tmp[id_order])), axis=0) ** 2).sum(axis=1))
+        bdy_brk[r], = np.where(bdy_dst[r] > 2 ** 0.5)
+        bdy_brk[r] += 1
+        bdy_brk[r] = np.insert(bdy_brk[r], 0, 0)  # insert start point
+        bdy_brk[r] = np.insert(bdy_brk[r], len(bdy_brk[r]), len(id_order))  # insert end point
+        bdy_dst[r] = np.insert(np.cumsum(bdy_dst[r]), 0, 0)
+
+        if r == 2:  # change to a valid rw number to get a visual output (outer most rw is zero)
+            f, ax = plt.subplots(nrows=1, ncols=1, figsize=(10, 10))
+            plt.scatter(nbi_tmp[bdy_ind[r][:]], nbj_tmp[bdy_ind[r][:]], s=1, marker='.')
+            for t in np.arange(0, len(bdy_ind[r]), 20):
+                plt.text(nbi_tmp[bdy_ind[r][t]], nbj_tmp[bdy_ind[r][t]], t)
+
+    return bdy_ind, bdy_dst, bdy_brk
+
+
+def plot_bdy(fname, bdy_ind, bdy_dst, bdy_brk, varnam, t, rw):
+    """
+    Determine the ordering and breaks in BDY files to aid plotting.
+    This function takes the i/j coordinates from BDY files and orders them sequentially
+    making it easier to visualise sections along the open boundary. Breaks in the open
+    boundary are also determined (i.e. where the distance between two points > 2**0.5)
+    Args:
+        fname     (str) : filename of BDY file
+    Returns:
+        bdy_ind   (dict): re-ordered indices
+        bdy_dst   (dict): distance (in model coords) between points
+        bdy_brk   (dict): location of the break points in the open boundary
+    """
+
+    # need to write in a check that either i or j are single values
+
+    rootgrp = Dataset(fname, "r", format="NETCDF4")
+    var = np.squeeze(rootgrp.variables[varnam][t, :])
+    nbr = np.squeeze(rootgrp.variables['nbrdta'][:, :]) - 1
+    var = var[:, nbr == rw]
+
+    # let us use gdept as the central depth irrespective of whether t, u or v
+
+    try:
+        gdep = np.squeeze(rootgrp.variables['deptht'][:, :, :])
+    except KeyError:
+        try:
+            gdep = np.squeeze(rootgrp.variables['gdept'][:, :, :])
+        except KeyError:
+            try:
+                gdep = np.squeeze(rootgrp.variables['depthu'][:, :, :])
+            except KeyError:
+                try:
+                    gdep = np.squeeze(rootgrp.variables['depthv'][:, :, :])
+                except KeyError:
+                    print 'depth variable not found'
+
+    rootgrp.close()
+
+    jpk = len(gdep[:, 0])
+    nsc = len(bdy_brk[rw][:]) - 1
+
+    dep = {}
+    dta = {}
+
+    # divide data up into sections and re-order
+
+    for n in range(nsc):
+        dta[n] = np.squeeze(var[:, bdy_ind[rw][bdy_brk[rw][n]:bdy_brk[rw][n + 1]]])
+        dep[n] = np.squeeze(gdep[:, bdy_ind[rw][bdy_brk[rw][n]:bdy_brk[rw][n + 1]]])
+
+    # loop over number of sections and plot
+
+    f, ax = plt.subplots(nrows=1, ncols=1, figsize=(11, 4))
+    ax.plot(dta[0][10, :])
+    ax.set_title('BDY points along section: 1, depth level: 10')
+
+    f, ax = plt.subplots(nrows=nsc, ncols=1, figsize=(14, 10 * nsc))
+
+    for n in range(nsc):
+        print('NSC '+str(n))
+        plt.sca(ax[n])
+
+        # from gdep create some pseudo w points
+
+        gdept = dep[n][:, :]
+        coords = np.arange(0, len(gdept[0, :]))
+        gdepw = np.zeros((len(gdept[:, 0]) + 1, len(gdept[0, :])))
+        for z in range(jpk):
+            gdepw[z + 1, :] = gdept[z, :] + (gdept[z, :] - gdepw[z, :])
+
+        gdepvw = np.zeros((len(gdept[:, 0]) + 1, len(gdept[0, :]) + 1))
+
+        # TODO: put in an adjustment for zps levels
+
+        gdepvw[:, 1:-1] = (gdepw[:, :-1] + gdepw[:, 1:]) / 2
+        gdepvw[:, 0] = gdepvw[:, 1]
+        gdepvw[:, -1] = gdepvw[:, -2]
+
+        # create a pseudo bathymetry from the depth data
+
+        bathy = np.zeros_like(coords)
+        mbath = np.sum(dta[n].mask == 0, axis=0)
+
+        for i in range(len(coords)):
+            bathy[i] = gdepw[mbath[i], i]
+
+        bathy_patch = Polygon(np.vstack((np.hstack((coords[0], coords, coords[-1])),
+                                         np.hstack((np.amax(bathy[:]), bathy, np.amax(bathy[:]))))).T,
+                              closed=True,
+                              facecolor=(0.8, 0.8, 0), alpha=0, edgecolor=None)
+
+        # Add patch to axes
+        ax[n].add_patch(bathy_patch)
+        ax[n].set_title('BDY points along section: ' + str(n))
+        patches = []
+        colors = []
+
+        for i in range(len(coords)):
+
+            for k in np.arange(jpk - 2, -1, -1):
+
+                if dta[n][k, i] > -10:
+                    x = [coords[i] - 0.5, coords[i], coords[i] + 0.5,
+                         coords[i] + 0.5, coords[i], coords[i] - 0.5, coords[i] - 0.5]
+                    y = [gdepvw[k + 1, i], gdepw[k + 1, i], gdepvw[k + 1, i + 1],
+                         gdepvw[k, i + 1], gdepw[k, i], gdepvw[k, i], gdepvw[k + 1, i]]
+
+                    polygon = Polygon(np.vstack((x, y)).T, True)
+                    patches.append(polygon)
+                    colors = np.append(colors, dta[n][k, i])
+
+        # for i in range(len(coords)):
+        #     #print(i)
+        #     for k in np.arange(jpk - 2, -1, -1):
+        #         x = [coords[i] - 0.5, coords[i], coords[i] + 0.5,
+        #              coords[i] + 0.5, coords[i], coords[i] - 0.5, coords[i] - 0.5]
+        #         y = [gdepvw[k + 1, i], gdepw[k + 1, i], gdepvw[k + 1, i + 1],
+        #              gdepvw[k, i + 1], gdepw[k, i], gdepvw[k, i], gdepvw[k + 1, i]]
+        #         plt.plot(x, y, 'k-', linewidth=0.1)
+        #         plt.plot(coords[i], gdept[k, i], 'k.', markersize=1)
+
+        plt.plot(coords, bathy, '-', color=(0.4, 0, 0))
+        p = PatchCollection(patches, alpha=0.8, edgecolor='none')
+        p.set_array(np.array(colors))
+        ax[n].add_collection(p)
+        f.colorbar(p, ax=ax[n])
+        ax[n].set_ylim((0, np.max(bathy)))
+        ax[n].invert_yaxis()
+
+    return f
+
+fname = '/Users/thopri/Projects/PyNEMO/outputs/NNA_R12_bdyT_y1979m11.nc'
+ind, dst, brk = nemo_bdy_order(fname)
+f = plot_bdy(fname, ind, dst, brk, 'votemper', 0, 0)
+
+plt.show()
\ No newline at end of file