From b79e26467e4525ac1a643a880c867635cf1fc43a Mon Sep 17 00:00:00 2001
From: thopri <thopri@144-148.noc.soton.ac.uk>
Date: Mon, 16 Mar 2020 13:49:27 +0000
Subject: [PATCH] added bathy and mask generation tools and got running with
orthoganal test grid
---
pynemo/nemo_bdy_extr_tm3.py | 2 +-
pynemo/reader/ncml.py | 2 +-
unit_tests/gen_tools.py | 210 ++++++++++++++++++++++------
unit_tests/namelist_unit_test.bdy | 10 +-
unit_tests/src_data_unit_tests.ncml | 6 -
unit_tests/test_gen.py | 24 ++--
6 files changed, 192 insertions(+), 62 deletions(-)
diff --git a/pynemo/nemo_bdy_extr_tm3.py b/pynemo/nemo_bdy_extr_tm3.py
index 9c457b7..4d7c745 100644
--- a/pynemo/nemo_bdy_extr_tm3.py
+++ b/pynemo/nemo_bdy_extr_tm3.py
@@ -831,7 +831,7 @@ class Extract:
del_t = time_counter[1] - time_counter[0]
dstep = 86400 / np.int(del_t)
-
+ dstep = int(dstep)
# TODO: put in a test to check all deltaT are the same otherwise throw
# an exception
diff --git a/pynemo/reader/ncml.py b/pynemo/reader/ncml.py
index 6a7b323..aba198d 100644
--- a/pynemo/reader/ncml.py
+++ b/pynemo/reader/ncml.py
@@ -280,7 +280,7 @@ class Variable(object):
attr = dvar.findAttributeIgnoreCase(attr_name)
if attr is not None:
retval = attr.getValue(0)
- return retval
+ return retval
except KeyError:
self.logger.error('Cannot find the requested variable '+self.variable)
return None
diff --git a/unit_tests/gen_tools.py b/unit_tests/gen_tools.py
index 80e1abe..38ba686 100644
--- a/unit_tests/gen_tools.py
+++ b/unit_tests/gen_tools.py
@@ -2,20 +2,17 @@
"""
Created on Wed Mar 11 11:16:57 2020
-example grid set up
+unit test grid generation tools
-@author: jdha
-"""
+@author: thopri and jdha (for example code)
-import numpy as np
+"""
from netCDF4 import Dataset
from math import cos,sin,radians
import numpy as np
-#import xarray as xr
import matplotlib.pyplot as plt
from math import radians
-#import cartopy
-#import cartopy.crs as ccrs
+
def set_hgrid(dx, dy, jpi, jpj, zoffx=0, zoffy=0,sf=1):
if sf > 1:
@@ -79,6 +76,8 @@ def set_zgrid(grid_h,jpk,max_dep,min_dep,z_dim):
dept_1d = np.linspace(min_dep,max_dep,jpk)
e3t_1d = np.linspace(1.0,z_dim,jpk)
e3w_1d = np.linspace(1.0,z_dim,jpk)
+ gdept_0 = np.linspace(min_dep,max_dep,jpk)
+ gdepw_0 = np.linspace(min_dep,max_dep,jpk)
e3t = np.zeros((jpi, jpj, len(e3t_1d)))
e3u = np.zeros((jpi, jpj, len(e3t_1d)))
@@ -87,6 +86,8 @@ def set_zgrid(grid_h,jpk,max_dep,min_dep,z_dim):
e3f = np.zeros((jpi, jpj, len(e3t_1d)))
e3uw = np.zeros((jpi, jpj, len(e3t_1d)))
e3vw = np.zeros((jpi, jpj, len(e3t_1d)))
+ gdept = np.zeros((jpi,jpj,len(gdept_0)))
+ gdepw = np.zeros((jpi,jpj,len(gdepw_0)))
e3t[:] = e3t_1d
e3u[:] = e3t_1d
@@ -95,9 +96,11 @@ def set_zgrid(grid_h,jpk,max_dep,min_dep,z_dim):
e3f[:] = e3t_1d
e3uw[:] = e3t_1d
e3vw[:] = e3t_1d
+ gdept[:] = gdept_0
+ gdepw[:] = gdepw_0
grid_z = {'dept_1d':dept_1d,'e3t_1d':e3t_1d,'e3w_1d':e3w_1d,'e3t':e3t,'e3u':e3u,'e3v':e3v,'e3w':e3w,'e3f':e3f, \
- 'e3uw':e3uw,'e3vw':e3vw}
+ 'e3uw':e3uw,'e3vw':e3vw,'gdept_0':gdept_0,'gdept':gdept,'gdepw_0':gdepw_0,'gdepw':gdepw}
return grid_z
@@ -121,33 +124,37 @@ def write_coord_H(fileout, grid_h):
# Get input size and create appropriate dimensions
# TODO: add some sort of error handling
nx, ny = np.shape(grid_h['lont'])
+ nt = 1
dataset.createDimension('x', nx)
dataset.createDimension('y', ny)
+ dataset.createDimension('t', nt)
# Create Variables
nav_lon = dataset.createVariable('nav_lon', np.float32, ('y', 'x'))
nav_lat = dataset.createVariable('nav_lat', np.float32, ('y', 'x'))
-
- glamt = dataset.createVariable('glamt', np.float64, ('y', 'x'))
- glamu = dataset.createVariable('glamu', np.float64, ('y', 'x'))
- glamv = dataset.createVariable('glamv', np.float64, ('y', 'x'))
- glamf = dataset.createVariable('glamf', np.float64, ('y', 'x'))
- gphit = dataset.createVariable('gphit', np.float64, ('y', 'x'))
- gphiu = dataset.createVariable('gphiu', np.float64, ('y', 'x'))
- gphiv = dataset.createVariable('gphiv', np.float64, ('y', 'x'))
- gphif = dataset.createVariable('gphif', np.float64, ('y', 'x'))
-
- ge1t = dataset.createVariable('e1t', np.float64, ('y', 'x'))
- ge1u = dataset.createVariable('e1u', np.float64, ('y', 'x'))
- ge1v = dataset.createVariable('e1v', np.float64, ('y', 'x'))
- ge1f = dataset.createVariable('e1f', np.float64, ('y', 'x'))
- ge2t = dataset.createVariable('e2t', np.float64, ('y', 'x'))
- ge2u = dataset.createVariable('e2u', np.float64, ('y', 'x'))
- ge2v = dataset.createVariable('e2v', np.float64, ('y', 'x'))
- ge2f = dataset.createVariable('e2f', np.float64, ('y', 'x'))
+ time_counter = dataset.createVariable('time_counter',np.float32,('t'))
+
+ glamt = dataset.createVariable('glamt', np.float64, ('t','y', 'x'))
+ glamu = dataset.createVariable('glamu', np.float64, ('t','y', 'x'))
+ glamv = dataset.createVariable('glamv', np.float64, ('t','y', 'x'))
+ glamf = dataset.createVariable('glamf', np.float64, ('t','y', 'x'))
+ gphit = dataset.createVariable('gphit', np.float64, ('t','y', 'x'))
+ gphiu = dataset.createVariable('gphiu', np.float64, ('t','y', 'x'))
+ gphiv = dataset.createVariable('gphiv', np.float64, ('t','y', 'x'))
+ gphif = dataset.createVariable('gphif', np.float64, ('t','y', 'x'))
+
+ ge1t = dataset.createVariable('e1t', np.float64, ('t','y', 'x'))
+ ge1u = dataset.createVariable('e1u', np.float64, ('t','y', 'x'))
+ ge1v = dataset.createVariable('e1v', np.float64, ('t','y', 'x'))
+ ge1f = dataset.createVariable('e1f', np.float64, ('t','y', 'x'))
+ ge2t = dataset.createVariable('e2t', np.float64, ('t','y', 'x'))
+ ge2u = dataset.createVariable('e2u', np.float64, ('t','y', 'x'))
+ ge2v = dataset.createVariable('e2v', np.float64, ('t','y', 'x'))
+ ge2f = dataset.createVariable('e2f', np.float64, ('t','y', 'x'))
nav_lon.units, nav_lon.long_name = 'km', 'X'
nav_lat.units, nav_lat.long_name = 'km', 'Y'
+ time_counter.units, time_counter.long_name = 'seconds','time_counter'
# Populate file with input data
# TODO: do we need to transpose?
@@ -192,29 +199,38 @@ def write_coord_Z(fileout, grid_h,grid_z):
# Get input size and create appropriate dimensions
# TODO: add some sort of error handling
nx, ny, nz = np.shape(grid_z['e3t'])
+ nt = 1
dataset.createDimension('x', nx)
dataset.createDimension('y', ny)
dataset.createDimension('z', nz)
+ dataset.createDimension('t', nt)
# Create Variables
nav_lon = dataset.createVariable('nav_lon', np.float32, ('y', 'x'))
nav_lat = dataset.createVariable('nav_lat', np.float32, ('y', 'x'))
nav_lev = dataset.createVariable('nav_lev', np.float32, 'z')
-
- ge3t1d = dataset.createVariable('e3t_1d', np.float64, 'z')
- ge3w1d = dataset.createVariable('e3w_1d', np.float64, 'z')
- gbat = dataset.createVariable('mbathy', np.float64, ('y', 'x'))
- ge3t = dataset.createVariable('e3t', np.float64, ('z','y', 'x'))
- ge3u = dataset.createVariable('e3u', np.float64, ('z','y', 'x'))
- ge3v = dataset.createVariable('e3v', np.float64, ('z','y', 'x'))
- ge3f = dataset.createVariable('e3f', np.float64, ('z','y', 'x'))
- ge3uw = dataset.createVariable('e3uw', np.float64, ('z', 'y', 'x'))
- ge3vw = dataset.createVariable('e3vw', np.float64, ('z', 'y', 'x'))
+ time_counter = dataset.createVariable('time_counter', np.float32, ('t'))
+
+ ge3t1d = dataset.createVariable('e3t_0', np.float64, ('t','z'))
+ ge3w1d = dataset.createVariable('e3w_0', np.float64, ('t','z'))
+ gdept_0 = dataset.createVariable('gdept_0', np.float64, ('t','z'))
+ gdepw_0 = dataset.createVariable('gdepw_0', np.float64, ('t','z'))
+ gbat = dataset.createVariable('mbathy', np.float64, ('t','y', 'x'))
+ ge3t = dataset.createVariable('e3t', np.float64, ('t','z','y', 'x'))
+ ge3u = dataset.createVariable('e3u', np.float64, ('t','z','y', 'x'))
+ ge3v = dataset.createVariable('e3v', np.float64, ('t','z','y', 'x'))
+ ge3w = dataset.createVariable('e3w', np.float64, ('t', 'z', 'y', 'x'))
+ ge3f = dataset.createVariable('e3f', np.float64, ('t','z','y', 'x'))
+ ge3uw = dataset.createVariable('e3uw', np.float64, ('t','z', 'y', 'x'))
+ ge3vw = dataset.createVariable('e3vw', np.float64, ('t','z', 'y', 'x'))
+ gdept = dataset.createVariable('gdept', np.float64, ('t','z', 'y', 'x'))
+ gdepw = dataset.createVariable('gdepw', np.float64, ('t','z', 'y', 'x'))
nav_lon.units, nav_lon.long_name = 'km', 'X'
nav_lat.units, nav_lat.long_name = 'km', 'Y'
nav_lev.units, nav_lev.long_name = 'm', 'Z'
+ time_counter.units, time_counter.long_name = 'seconds', 'time_counter'
# Populate file with input data
# TODO: do we need to transpose?
@@ -226,12 +242,17 @@ def write_coord_Z(fileout, grid_h,grid_z):
ge3u[:, :, :] = grid_z['e3u'].T
ge3v[:, :, :] = grid_z['e3v'].T
ge3f[:, :, :] = grid_z['e3f'].T
+ ge3w[:, :, :] = grid_z['e3w'].T
ge3uw[:, :, :] = grid_z['e3uw'].T
ge3vw[:, :, :] = grid_z['e3vw'].T
ge3t1d[:] = grid_z['e3t_1d']
ge3w1d[:] = grid_z['e3w_1d']
gbat[:,:] = grid_h['batt'].T
+ gdept[:,:,:] = grid_z['gdept'].T
+ gdept_0[:] = grid_z['gdept_0']
+ gdepw[:,:,:] = grid_z['gdepw'].T
+ gdepw_0[:] = grid_z['gdepw_0']
# Close off pointer
@@ -446,11 +467,73 @@ def plot_grids(lat_in,lon_in,new_lat,new_lon,off_lat,off_lon,src_lat,src_lon):
ax1.plot(src_lon,src_lat, color='green',marker='o',linestyle="")
ax2.plot(off_lon,off_lat, color='yellow',marker='.',linestyle="")
ax2.plot(src_lon, src_lat, color='green', marker='o', linestyle="")
- # tweak margins of subplots as tight layout doesn't work
+ # tweak margins of subplots as tight layout doesn't work for some reason?
plt.subplots_adjust(left=0.01, right=1, top=0.9, bottom=0.05,wspace=0.01)
plt.show()
+def write_mask(fileout,grid_h,grid_z):
+ '''
+ Writes out a
+
+ Args:
+
+ Returns:
+ '''
+ # Open pointer to netcdf file
+ dataset = Dataset(fileout, 'w', format='NETCDF4_CLASSIC')
+
+ # Get input size and create appropriate dimensions
+ # TODO: add some sort of error handling
+ nt = 1
+ nx, ny, nz = np.shape(grid_z['e3t'])
+ dataset.createDimension('x', nx)
+ dataset.createDimension('y', ny)
+ dataset.createDimension('z', nz)
+ dataset.createDimension('t', nt)
+ # Create Variables
+ nav_lon = dataset.createVariable('nav_lon', np.float32, ('y', 'x'))
+ nav_lat = dataset.createVariable('nav_lat', np.float32, ('y', 'x'))
+ nav_lev = dataset.createVariable('nav_lev', np.float32, 'z')
+ time_counter = dataset.createVariable('time_counter', np.float32, ('t'))
+
+ fmaskutil = dataset.createVariable('fmaskutil', np.float64, ('t','y', 'x'))
+ tmaskutil = dataset.createVariable('tmaskutil', np.float64, ('t','y', 'x'))
+ umaskutil = dataset.createVariable('umaskutil', np.float64, ('t','y', 'x'))
+ vmaskutil = dataset.createVariable('vmaskutil', np.float64, ('t','y', 'x'))
+
+ fmask = dataset.createVariable('fmask', np.float64, ('t','z','y', 'x'))
+ tmask = dataset.createVariable('tmask', np.float64, ('t','z','y', 'x'))
+ umask = dataset.createVariable('umask', np.float64, ('t','z','y', 'x'))
+ vmask = dataset.createVariable('vmask', np.float64, ('t','z','y', 'x'))
+
+ nav_lon.units, nav_lon.long_name = 'km', 'X'
+ nav_lat.units, nav_lat.long_name = 'km', 'Y'
+ nav_lev.units, nav_lev.long_name = 'm', 'Z'
+ time_counter.units, time_counter.long_name = 'seconds', 'time_counter'
+
+ # Populate file with input data
+ # TODO: do we need to transpose?
+ nav_lon[:, :] = grid_h['lont'].T
+ nav_lat[:, :] = grid_h['latt'].T
+ nav_lev[:] = grid_z['dept_1d']
+
+ threeD_mask = np.zeros(np.shape(grid_z['e3t']))
+ twoD_mask = np.zeros(np.shape(grid_h['e1t']))
+
+ fmask[:,:,:] = threeD_mask.T
+ fmaskutil[:,:] = twoD_mask.T
+ tmask[:,:,:] = threeD_mask.T
+ tmaskutil[:,:] = twoD_mask.T
+ umask[:,:,:] = threeD_mask.T
+ umaskutil[:,:] = twoD_mask.T
+ vmask[:,:,:] = threeD_mask.T
+ vmaskutil[:,:] = twoD_mask.T
+
+ dataset.close()
+
+ return 0
+
def write_parameter(fileout, grid_h,grid_z,params):
'''
Writes out a
@@ -460,31 +543,41 @@ def write_parameter(fileout, grid_h,grid_z,params):
Returns:
'''
- #TODO: implement multiple parameters using dicts.
+ #TODO: implement multiple parameters using params dict.
# Open pointer to netcdf file
dataset = Dataset(fileout, 'w', format='NETCDF4_CLASSIC')
# Get input size and create appropriate dimensions
# TODO: add some sort of error handling
+ nt = 2
nx, ny, nz = np.shape(grid_z['e3t'])
dataset.createDimension('x', nx)
dataset.createDimension('y', ny)
dataset.createDimension('z', nz)
+ dataset.createDimension('time', nt)
+
# Create Variables
longitude = dataset.createVariable('longitude', np.float32, ('y', 'x'))
latitude = dataset.createVariable('latitude', np.float32, ('y', 'x'))
depth = dataset.createVariable('depth', np.float32, 'z')
+ time_counter = dataset.createVariable('time', np.float32, ('time'))
longitude.units, longitude.long_name = 'km', 'X'
latitude.units, latitude.long_name = 'km', 'Y'
depth.units, depth.long_name = 'm', 'Z'
+ time_counter.units = 'hours since 1950-01-01 00:00:00'
+ time_counter.long_name = 'Time (hours since 1950-01-01)'
+ time_counter.axis = 'T'
+ time_counter._CoordinateAxisType = "Time"
+ time_counter.calendar = 'gregorian'
+ time_counter.standard_name = 'time'
# Populate file with input data
longitude[:, :] = grid_h['lont'].T
latitude[:, :] = grid_h['latt'].T
depth[:] = grid_z['dept_1d']
-
- parameter = dataset.createVariable(str(params['name']), np.float64, ('z', 'y', 'x'))
+ time_counter[:] = (587340.00,588060.00)
+ parameter = dataset.createVariable(str(params['name']), np.float64, ('time','z', 'y', 'x'))
parameter.units, parameter.long_name = str(params['units']), str(params['longname'])
value_fill = np.ones(np.shape(grid_z['e3t']))
value_fill = value_fill*params['const_value']
@@ -492,5 +585,42 @@ def write_parameter(fileout, grid_h,grid_z,params):
# Close off pointer
dataset.close()
+ return 0
+
+def write_bathy(fileout, grid_h,grid_z):
+ '''
+ Writes out a
+
+ Args:
+ Returns:
+ '''
+ # Open pointer to netcdf file
+ dataset = Dataset(fileout, 'w', format='NETCDF4_CLASSIC')
+
+ # Get input size and create appropriate dimensions
+ # TODO: add some sort of error handling
+ nx, ny = np.shape(grid_h['e1t'])
+ dataset.createDimension('x', nx)
+ dataset.createDimension('y', ny)
+
+
+ # Create Variables
+ nav_lon = dataset.createVariable('nav_lon', np.float32, ('y', 'x'))
+
+ nav_lat = dataset.createVariable('nav_lat', np.float32, ('y', 'x'))
+ nav_lon.units, nav_lon.long_name = 'km', 'X'
+ nav_lat.units, nav_lat.long_name = 'km', 'Y'
+
+ Bathymetry = dataset.createVariable('Bathymetry', np.float64,('y','x'))
+ Bathymetry.units,Bathymetry.long_name = 'meters','Median depth by area'
+ nav_lon[:, :] = grid_h['lont'].T
+ nav_lat[:, :] = grid_h['latt'].T
+
+ Bathy = np.ones((nx,ny))
+ Bathy = Bathy * grid_z['dept_1d'][-1]
+
+ Bathymetry[:,:] = Bathy
+
+ dataset.close()
return 0
\ No newline at end of file
diff --git a/unit_tests/namelist_unit_test.bdy b/unit_tests/namelist_unit_test.bdy
index b536895..68539ce 100644
--- a/unit_tests/namelist_unit_test.bdy
+++ b/unit_tests/namelist_unit_test.bdy
@@ -30,11 +30,11 @@
!------------------------------------------------------------------------------
! grid information
!------------------------------------------------------------------------------
- sn_src_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/src_coordinates.nc'
- sn_src_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/src_coordinates.nc'
- sn_dst_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/dst_hgr_zps.nc'
- sn_dst_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/dst_zgr_zps.nc'
- sn_src_msk = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/src_bathy.nc'
+ sn_src_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_src_hgr_zps.nc'
+ sn_src_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_src_zgr_zps.nc'
+ sn_dst_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_dst_hgr_zps.nc'
+ sn_dst_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/test_dst_zgr_zps.nc'
+ sn_src_msk = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/mask.nc'
sn_bathy = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/dst_bathy.nc'
!------------------------------------------------------------------------------
diff --git a/unit_tests/src_data_unit_tests.ncml b/unit_tests/src_data_unit_tests.ncml
index 2b91153..3016c52 100644
--- a/unit_tests/src_data_unit_tests.ncml
+++ b/unit_tests/src_data_unit_tests.ncml
@@ -5,12 +5,6 @@
<ns0:scan location="file://Users/thopri/Projects/PyNEMO/unit_tests/test_data/" regExp=".*T\.nc$" />
</ns0:aggregation>
</ns0:netcdf>
- <ns0:netcdf>
- <ns0:aggregation dimName="time" name="sea_surface_height" type="joinExisting">
- <ns0:scan location="file://Users/thopri/Projects/PyNEMO/unit_tests/test_data/" regExp=".*T\.nc$" />
- </ns0:aggregation>
- </ns0:netcdf>
</ns0:aggregation>
<ns0:variable name="thetao" orgName="thetao" />
- <ns0:variable name="zos" orgName="zos" />
</ns0:netcdf>
diff --git a/unit_tests/test_gen.py b/unit_tests/test_gen.py
index 95c2984..15ec96e 100644
--- a/unit_tests/test_gen.py
+++ b/unit_tests/test_gen.py
@@ -90,22 +90,28 @@ def _main():
if write_coord_H + write_coord_Z == 0:
print("Success!")
-
# plot orginal, rotatated and source lat and lon
gt.plot_grids(grid_h2['latt'],grid_h2['lont'],grid_rot['latt'],grid_rot['lont'],grid_h3['latt'], \
grid_h3['lont'],grid_h1['latt'],grid_h1['lont'])
- # save new lat lon to dataset
- #ds.nav_lat.values[:], ds.nav_lon.values[:] = new_lat,new_lon
- # save dataset to netcdf
- #ds.to_netcdf('unit_tests/test_data/rot_'+str(rot)+'_dst_hgr_zps.nc')
- # close data files
- #ds.close()
- #src.close()
- out_fname = 'unit_tests/test_data/output_boundary.nc'
+
+ # write bathy files (constant bathy)
+ bathy_fname = 'unit_tests/test_data/dst_bathy.nc'
+ bathy = gt.write_bathy(bathy_fname,grid_h2,grid_z2)
+ if bathy == 0:
+ print('Success!')
+
+ # write boundary files (constant parameters)
+ out_fname = 'unit_tests/test_data/output_boundary_T.nc'
params = {'name':'thetao','const_value':15.0,'longname':'temperature','units':'degreesC'}
boundary = gt.write_parameter(out_fname,grid_h1,grid_z1,params)
if boundary == 0:
print('Success!')
+ #write_mask
+ mask_fname = 'unit_tests/test_data/mask.nc'
+ mask = gt.write_mask(mask_fname,grid_h1,grid_z1)
+ if mask == 0:
+ print('Success!')
+
if __name__ == '__main__':
_main()
\ No newline at end of file
--
GitLab