Unverified Commit 6f50bc0c authored by jdha's avatar jdha Committed by GitHub
Browse files

Merge pull request #1 from jdha/PyNEMO3

Py nemo3
parents b1ee0b2f a666111c
...@@ -7,8 +7,8 @@ Used for development purposes to display the ncml editor dialog. ...@@ -7,8 +7,8 @@ Used for development purposes to display the ncml editor dialog.
@author: Shirley Crompton, UK Science and Technology Facilities Council @author: Shirley Crompton, UK Science and Technology Facilities Council
''' '''
import sys import sys
from PyQt4.QtGui import * from PyQt5.QtWidgets import *
from gui import nemo_ncml_generator as ncml_generator from .gui import nemo_ncml_generator as ncml_generator
import logging import logging
# Logging set to info # Logging set to info
logging.basicConfig(level=logging.INFO) logging.basicConfig(level=logging.INFO)
......
...@@ -5,10 +5,10 @@ Created on 7 Jan 2015 ...@@ -5,10 +5,10 @@ Created on 7 Jan 2015
''' '''
# pylint: disable=E1103 # pylint: disable=E1103
# pylint: disable=no-name-in-module # pylint: disable=no-name-in-module
from PyQt4 import QtGui from PyQt5 import QtGui, QtWidgets
from gui.nemo_bdy_input_window import InputWindow from .gui.nemo_bdy_input_window import InputWindow
import nemo_bdy_setup from . import nemo_bdy_setup
import sys, getopt import sys, getopt
...@@ -16,9 +16,9 @@ def open_settings_window(fname): ...@@ -16,9 +16,9 @@ def open_settings_window(fname):
""" Main method which starts a Qt application and gives user """ Main method which starts a Qt application and gives user
an option to pick a namelist.bdy file to edit. Once user selects it an option to pick a namelist.bdy file to edit. Once user selects it
it will open a dialog box where users can edit the parameters""" it will open a dialog box where users can edit the parameters"""
app = QtGui.QApplication(sys.argv) app = QtWidgets.QApplication(sys.argv)
if fname is None: if fname is None:
fname = QtGui.QFileDialog.getOpenFileName(None, 'Open file') fname = QtWidgets.QFileDialog.getOpenFileName(None, 'Open file')
setup = nemo_bdy_setup.Setup(fname)#'../../data/namelisttest.bdy') setup = nemo_bdy_setup.Setup(fname)#'../../data/namelisttest.bdy')
ex = InputWindow(setup) ex = InputWindow(setup)
...@@ -29,7 +29,7 @@ def open_settings_dialog(setup): ...@@ -29,7 +29,7 @@ def open_settings_dialog(setup):
""" This method is to start the settings window using the setup settings provided """ This method is to start the settings window using the setup settings provided
in the input. On clicking the cancel button it doesn't shutdown the applicaiton in the input. On clicking the cancel button it doesn't shutdown the applicaiton
but carries on with the execution""" but carries on with the execution"""
app = QtGui.QApplication(sys.argv) app = QtWidgets.QApplication(sys.argv)
ex = InputWindow(setup) ex = InputWindow(setup)
ex.nl_editor.btn_cancel.clicked.connect(app.quit) ex.nl_editor.btn_cancel.clicked.connect(app.quit)
return app.exec_(), ex.mpl_widget.mask return app.exec_(), ex.mpl_widget.mask
...@@ -41,12 +41,12 @@ def main(): ...@@ -41,12 +41,12 @@ def main():
try: try:
opts, dummy_args = getopt.getopt(sys.argv[1:], "hs:", ["help", "setup="]) opts, dummy_args = getopt.getopt(sys.argv[1:], "hs:", ["help", "setup="])
except getopt.GetoptError: except getopt.GetoptError:
print "usage: pynemo_settings_editor -s <namelist.bdy> " print("usage: pynemo_settings_editor -s <namelist.bdy> ")
sys.exit(2) sys.exit(2)
for opt, arg in opts: for opt, arg in opts:
if opt == "-h": if opt == "-h":
print "usage: pynemo_settings_editor -s <namelist.bdy> " print("usage: pynemo_settings_editor -s <namelist.bdy> ")
sys.exit() sys.exit()
elif opt in ("-s", "--setup"): elif opt in ("-s", "--setup"):
setup_file = arg setup_file = arg
......
...@@ -5,7 +5,7 @@ This is an abstraction for the data repository ...@@ -5,7 +5,7 @@ This is an abstraction for the data repository
from os import listdir from os import listdir
import numpy as np import numpy as np
from netCDF4 import Dataset from netCDF4 import Dataset
from netCDF4 import netcdftime from cftime import utime
import copy import copy
import logging import logging
class Reader(object): class Reader(object):
...@@ -50,7 +50,7 @@ class Reader(object): ...@@ -50,7 +50,7 @@ class Reader(object):
dir_list[i] = self.directory + dir_list[i] dir_list[i] = self.directory + dir_list[i]
dir_list.sort() dir_list.sort()
return filter(None, dir_list) return [_f for _f in dir_list if _f]
def _delta_time_interval(self, time1, time2): def _delta_time_interval(self, time1, time2):
""" Get the difference between the two times in days and hours""" """ Get the difference between the two times in days and hours"""
...@@ -74,7 +74,7 @@ class Reader(object): ...@@ -74,7 +74,7 @@ class Reader(object):
x = [filename, index] x = [filename, index]
group.data_list.append(x) group.data_list.append(x)
group.time_counter.append(varid[index]+t_adjust) group.time_counter.append(varid[index]+t_adjust)
group.date_counter.append(netcdftime.utime(varid.units, group.date_counter.append(utime(varid.units,
varid.calendar).num2date(varid[index]+t_adjust)) varid.calendar).num2date(varid[index]+t_adjust))
group.units = varid.units group.units = varid.units
group.calendar = varid.calendar group.calendar = varid.calendar
...@@ -94,7 +94,7 @@ class Reader(object): ...@@ -94,7 +94,7 @@ class Reader(object):
error""" error"""
days = set() days = set()
hrs = set() hrs = set()
for grid_type in self.grid_source_data.keys(): for grid_type in list(self.grid_source_data.keys()):
day, hr = self._delta_time_interval(self.grid_source_data[grid_type].time_counter[0], day, hr = self._delta_time_interval(self.grid_source_data[grid_type].time_counter[0],
self.grid_source_data[grid_type].time_counter[1]) self.grid_source_data[grid_type].time_counter[1])
days.add(day) days.add(day)
......
...@@ -15,13 +15,13 @@ from netCDF4 import Dataset ...@@ -15,13 +15,13 @@ from netCDF4 import Dataset
def GetReader(uri, t_adjust, reader_type=None): def GetReader(uri, t_adjust, reader_type=None):
if reader_type is None: if reader_type is None:
print uri print(uri)
if uri.endswith(".ncml"): if uri.endswith(".ncml"):
reader_type = "NcML" reader_type = "NcML"
elif os.path.isdir(uri): elif os.path.isdir(uri):
reader_type = "Directory" reader_type = "Directory"
else: else:
print "Error input should be a NcML file or URL or a Local directory" print("Error input should be a NcML file or URL or a Local directory")
return None return None
if reader_type == "NcML": if reader_type == "NcML":
return NcMLReader(uri,t_adjust) return NcMLReader(uri,t_adjust)
......
...@@ -10,7 +10,7 @@ import string ...@@ -10,7 +10,7 @@ import string
import logging import logging
import numpy as np import numpy as np
import jnius_config import jnius_config
from netCDF4 import netcdftime from cftime import utime
ncmlpath, file_name = os.path.split(__file__) ncmlpath, file_name = os.path.split(__file__)
ncmlpath = os.path.join(ncmlpath, "jars", "netcdfAll-4.6.jar") ncmlpath = os.path.join(ncmlpath, "jars", "netcdfAll-4.6.jar")
jnius_config.set_classpath('.',ncmlpath) jnius_config.set_classpath('.',ncmlpath)
...@@ -22,7 +22,7 @@ try: ...@@ -22,7 +22,7 @@ try:
proxy_port = proxylist[1] proxy_port = proxylist[1]
jnius_config.add_options('-Dhttp.proxyHost='+proxy_host,'-Dhttp.proxyPort='+proxy_port) jnius_config.add_options('-Dhttp.proxyHost='+proxy_host,'-Dhttp.proxyPort='+proxy_port)
except: except:
print "Didn't find a proxy environment variable" print("Didn't find a proxy environment variable")
NetcdfDataset = None NetcdfDataset = None
NcMLReader = None NcMLReader = None
Section = None Section = None
...@@ -37,7 +37,7 @@ try: ...@@ -37,7 +37,7 @@ try:
Section = autoclass('ucar.ma2.Section') Section = autoclass('ucar.ma2.Section')
init_jnius() init_jnius()
except ImportError: except ImportError:
print 'Warning: Please make sure pyjnius is installed and jvm.dll/libjvm.so/libjvm.dylib is in the path' print('Warning: Please make sure pyjnius is installed and jvm.dll/libjvm.so/libjvm.dylib is in the path')
time_counter_const = "time_counter" time_counter_const = "time_counter"
class Reader(object): class Reader(object):
...@@ -78,7 +78,7 @@ class Reader(object): ...@@ -78,7 +78,7 @@ class Reader(object):
grid.time_counter = timevar[:]+t_adjust grid.time_counter = timevar[:]+t_adjust
grid.date_counter = [] grid.date_counter = []
for index in range(0,len(grid.time_counter)): for index in range(0,len(grid.time_counter)):
grid.date_counter.append(netcdftime.utime(grid.units, grid.date_counter.append(utime(grid.units,
grid.calendar).num2date(grid.time_counter[index])) grid.calendar).num2date(grid.time_counter[index]))
def close(self): def close(self):
......
...@@ -23,8 +23,8 @@ class Test(unittest.TestCase): ...@@ -23,8 +23,8 @@ class Test(unittest.TestCase):
self.lon,self.lat = np.meshgrid(self.lon, self.lat) self.lon,self.lat = np.meshgrid(self.lon, self.lat)
#gcoms_break_depth.gcoms_boundary_masks(self.bathy, -1, 0) #gcoms_break_depth.gcoms_boundary_masks(self.bathy, -1, 0)
print self.bathy.shape print(self.bathy.shape)
print self.bathy print(self.bathy)
pass pass
...@@ -39,10 +39,10 @@ class Test(unittest.TestCase): ...@@ -39,10 +39,10 @@ class Test(unittest.TestCase):
self.bathy[self.bathy>=0] = 0 self.bathy[self.bathy>=0] = 0
self.bathy = self.bathy*-1 self.bathy = self.bathy*-1
tmp = gcoms_break_depth.polcoms_select_domain(self.bathy, self.lat, self.lon, roi, 200) tmp = gcoms_break_depth.polcoms_select_domain(self.bathy, self.lat, self.lon, roi, 200)
self.assertEquals(tmp[32,0],1,"Set the break select correctly") self.assertEqual(tmp[32,0],1,"Set the break select correctly")
self.assertEquals(tmp[40,0],1,"Set the break select correctly 40") self.assertEqual(tmp[40,0],1,"Set the break select correctly 40")
self.assertEquals(tmp[50,0],1,"Set the break select correctly 50") self.assertEqual(tmp[50,0],1,"Set the break select correctly 50")
self.assertEquals(tmp[60,0],1,"Set the break select correctly 60") self.assertEqual(tmp[60,0],1,"Set the break select correctly 60")
def testGcomsBreakDepth(self): def testGcomsBreakDepth(self):
r = 18 r = 18
......
...@@ -11,8 +11,8 @@ class Test(unittest.TestCase): ...@@ -11,8 +11,8 @@ class Test(unittest.TestCase):
def testMaskData(self): def testMaskData(self):
mask = Mask('data/grid_C/NNA_R12_bathy_meter_bench.nc') mask = Mask('data/grid_C/NNA_R12_bathy_meter_bench.nc')
print mask.data.shape print(mask.data.shape)
self.assertEqual(mask.data.shape,(401L,351L),'Mask reading failed') self.assertEqual(mask.data.shape,(401,351),'Mask reading failed')
if __name__ == "__main__": if __name__ == "__main__":
......
...@@ -45,7 +45,7 @@ class Test(unittest.TestCase): ...@@ -45,7 +45,7 @@ class Test(unittest.TestCase):
#test #test
setup = Setup('test.bdy') setup = Setup('test.bdy')
self.assertEquals(setup.settings,{'nonambigious':True},"Didn't recognize valid setting") self.assertEqual(setup.settings,{'nonambigious':True},"Didn't recognize valid setting")
#delete file #delete file
os.remove('test.bdy') os.remove('test.bdy')
...@@ -64,12 +64,12 @@ class Test(unittest.TestCase): ...@@ -64,12 +64,12 @@ class Test(unittest.TestCase):
#test #test
setup = Setup('test.bdy') setup = Setup('test.bdy')
print setup.settings print(setup.settings)
self.assertEquals(setup.settings['nonambigious'],False,"Didn't recognize logical setting") self.assertEqual(setup.settings['nonambigious'],False,"Didn't recognize logical setting")
self.assertEquals(setup.settings['floatval'],10.9,"Didn't recognize rn value in setting") self.assertEqual(setup.settings['floatval'],10.9,"Didn't recognize rn value in setting")
self.assertEquals(setup.settings['floatval2'],20.9,"Didn't recongnize nn value in setting") self.assertEqual(setup.settings['floatval2'],20.9,"Didn't recongnize nn value in setting")
self.assertEquals(setup.settings['stringval'],'coordinates.nc',"Didn't recognize cn string value in setting") self.assertEqual(setup.settings['stringval'],'coordinates.nc',"Didn't recognize cn string value in setting")
self.assertEquals(setup.settings['stringval2'],'gregorian',"Didn't recognize sn string value in setting") self.assertEqual(setup.settings['stringval2'],'gregorian',"Didn't recognize sn string value in setting")
#delete file #delete file
os.remove('test.bdy') os.remove('test.bdy')
......
...@@ -34,7 +34,7 @@ class Test(unittest.TestCase): ...@@ -34,7 +34,7 @@ class Test(unittest.TestCase):
testfile = os.path.join(testpath, "testremote.ncml") testfile = os.path.join(testpath, "testremote.ncml")
sd = Reader(testfile,0).grid sd = Reader(testfile,0).grid
dataset = Data2(sd.dataset,"time_counter") dataset = Data2(sd.dataset,"time_counter")
self.assertEquals(len(dataset), 8, "There should 8 datasets") self.assertEqual(len(dataset), 8, "There should 8 datasets")
def testVariable(self): def testVariable(self):
init_jnius() init_jnius()
...@@ -80,11 +80,11 @@ class Test(unittest.TestCase): ...@@ -80,11 +80,11 @@ class Test(unittest.TestCase):
testpath, file_name = os.path.split(__file__) testpath, file_name = os.path.split(__file__)
testfile = os.path.join(testpath, "testremote.ncml") testfile = os.path.join(testpath, "testremote.ncml")
repo = Reader(testfile,0) repo = Reader(testfile,0)
self.assertEquals(len(repo['t'].time_counter), 8, "Time counter should be 8") self.assertEqual(len(repo['t'].time_counter), 8, "Time counter should be 8")
self.assertEquals(repo['t'].time_counter[0], 691416000, "The first time value doesn't match") self.assertEqual(repo['t'].time_counter[0], 691416000, "The first time value doesn't match")
repo = Reader(testfile,100) repo = Reader(testfile,100)
self.assertEquals(len(repo['t'].time_counter), 8, "Time counter should be 8") self.assertEqual(len(repo['t'].time_counter), 8, "Time counter should be 8")
self.assertEquals(repo['t'].time_counter[0], 691416100, "The first time value doesn't match") self.assertEqual(repo['t'].time_counter[0], 691416100, "The first time value doesn't match")
if __name__ == "__main__": if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testName'] #import sys;sys.argv = ['', 'Test.testName']
......
...@@ -77,8 +77,8 @@ class Extract: ...@@ -77,8 +77,8 @@ class Extract:
nb.close() # Close Bathymetry file nb.close() # Close Bathymetry file
# Find nearest neighbours on the source grid to each dst bdy point # Find nearest neighbours on the source grid to each dst bdy point
source_tree = sp.cKDTree(zip(lon, lat)) source_tree = sp.cKDTree(list(zip(lon, lat)))
dst_pts = zip(dst_lon, dst_lat) dst_pts = list(zip(dst_lon, dst_lat))
nn_dist, self.nn_id = source_tree.query(dst_pts, k=4, eps=0, p=2, nn_dist, self.nn_id = source_tree.query(dst_pts, k=4, eps=0, p=2,
distance_upper_bound=0.5) distance_upper_bound=0.5)
...@@ -106,7 +106,7 @@ class Extract: ...@@ -106,7 +106,7 @@ class Extract:
# Need to identify missing points and throw a warning and set values to zero # Need to identify missing points and throw a warning and set values to zero
mv = np.sum(self.wei,axis=1) == 0 mv = np.sum(self.wei,axis=1) == 0
print '##WARNING## There are', np.sum(mv), 'missing values, these will be set to ZERO' print('##WARNING## There are', np.sum(mv), 'missing values, these will be set to ZERO')
def extract_con(self, con): def extract_con(self, con):
...@@ -150,7 +150,7 @@ class Extract: ...@@ -150,7 +150,7 @@ class Extract:
else: else:
# throw some warning # throw some warning
print '##WARNING## Missing constituent values will be set to ZERO' print('##WARNING## Missing constituent values will be set to ZERO')
self.harm_Im[con] = np.zeros(self.nn_id[:,0].size) self.harm_Im[con] = np.zeros(self.nn_id[:,0].size)
self.harm_Re[con] = np.zeros(self.nn_id[:,0].size) self.harm_Re[con] = np.zeros(self.nn_id[:,0].size)
......
''' Non
\ No newline at end of file
'''
# pylint: disable=E1103
# pylint: disable=no-name-in-module
import numpy as np
import scipy.spatial as sp
from netCDF4 import Dataset
import copy # DEBUG ONLY- allows multiple runs without corruption
from pynemo import nemo_bdy_grid_angle
from pynemo.nemo_bdy_lib import rot_rep
class Extract:
def __init__(self, setup, DstCoord, Grid):
self.g_type = Grid.grid_type
DC = copy.deepcopy(DstCoord)
dst_lon = DC.bdy_lonlat[self.g_type]['lon'][Grid.bdy_r == 0]
dst_lat = DC.bdy_lonlat[self.g_type]['lat'][Grid.bdy_r == 0]
self.dst_dep = DC.depths[self.g_type]['bdy_hbat'][Grid.bdy_r == 0]
self.harm_Im = {} # tidal boundary data: Imaginary
self.harm_Re = {} # tidal boundary data: Real
# Modify lon for 0-360 TODO this needs to be auto-dectected
dst_lon = np.array([x if x > 0 else x+360 for x in dst_lon])
fileIDb = '../data/tide/grid_tpxo7.2.nc' # TPX bathymetry file
nb = Dataset(fileIDb) # Open the TPX bathybetry file using the NetCDF4-Python library
# Open the TPX Datafiles using the NetCDF4-Python library
if self.g_type == 't':
self.fileID = '../data/tide/h_tpxo7.2.nc' # TPX sea surface height file
self.var_Im = 'hIm'
self.var_Re = 'hRe'
elif (self.g_type == 'u') or (self.g_type == 'v') :
self.fileID = '../data/tide/u_tpxo7.2.nc' # TPX velocity file
self.var_Im = 'UIm'
self.var_Re = 'URe'
self.var_Im2 = 'VIm'
self.var_Re2 = 'VRe'
self.key_tr = setup['trans']
# Determine the grid angle for rotating vector qunatities
maxJ = DC.lonlat['t']['lon'].shape[0]
maxI = DC.lonlat['t']['lon'].shape[1]
GridAngles = nemo_bdy_grid_angle.GridAngle(setup['dst_hgr'], 1, maxI, 1, maxJ, self.g_type)
dst_gcos = np.ones([maxJ, maxI])
dst_gsin = np.zeros([maxJ,maxI])
dst_gcos[1:,1:] = GridAngles.cosval
dst_gsin[1:,1:] = GridAngles.sinval
# Retain only boundary points rotation information
self.gcos = np.zeros(Grid.bdy_i.shape[0])
self.gsin = np.zeros(Grid.bdy_i.shape[0])
for p in range(Grid.bdy_i.shape[0]):
self.gcos[p] = dst_gcos[Grid.bdy_i[p,1], Grid.bdy_i[p,0]]
self.gsin[p] = dst_gsin[Grid.bdy_i[p,1], Grid.bdy_i[p,0]]
if self.g_type == 'u':
self.rot_dir = 'i'
elif self.g_type == 'v':
self.rot_dir = 'j'
else:
print 'You should not see this message!'
# We will average velocities onto the T grid as there is a rotation to be done
# Also need to account for east-west wrap-around
nc = Dataset('../data/tide/h_tpxo7.2.nc')
lon = np.ravel(np.concatenate([nc.variables['lon_z'][-2:,:],
nc.variables['lon_z'][:,:],
nc.variables['lon_z'][0:2,:]]))
lat = np.ravel(np.concatenate([nc.variables['lat_z'][-2:,:],
nc.variables['lat_z'][:,:],
nc.variables['lat_z'][0:2,:]]))
bat = np.ravel(np.concatenate([nb.variables['hz'][-2:,:],
nb.variables['hz'][:,:],
nb.variables['hz'][0:2,:]]))
msk = np.ravel(np.concatenate([nb.variables['mz'][-2:,:],
nb.variables['mz'][:,:],
nb.variables['mz'][0:2,:]]))
# Pull out the constituents that are avaibable
self.cons = []
for ncon in range(nc.variables['con'].shape[0]):
self.cons.append(nc.variables['con'][ncon,:].tostring().strip())
nc.close() # Close Datafile
nb.close() # Close Bathymetry file
# Find nearest neighbours on the source grid to each dst bdy point
source_tree = sp.cKDTree(zip(lon, lat))
dst_pts = zip(dst_lon, dst_lat)
# Upper bound set at 0.5 deg as the TPXO7.2 data are at 0.25 deg resolution and
# we don't want to grab points from further afield
nn_dist, self.nn_id = source_tree.query(dst_pts, k=4, eps=0, p=2,
distance_upper_bound=0.5)
# Create a weighting index for interpolation onto dst bdy point
# need to check for missing values
ind = nn_dist == np.inf
self.nn_id[ind] = 0 # better way of carrying None in the indices?
dx = (lon[self.nn_id] - np.repeat(np.reshape(dst_lon,[dst_lon.size, 1]),4,axis=1) ) * np.cos(np.repeat(np.reshape(dst_lat,[dst_lat.size, 1]),4,axis=1) * np.pi / 180.)
dy = lat[self.nn_id] - np.repeat(np.reshape(dst_lat,[dst_lat.size, 1]),4,axis=1)
dist_tot = np.power((np.power(dx, 2) + np.power(dy, 2)), 0.5)
self.msk = msk[self.nn_id]
self.bat = bat[self.nn_id]
dist_tot[ind | self.msk] = np.nan
dist_wei = 1/( np.divide(dist_tot,(np.repeat(np.reshape(np.nansum(dist_tot,axis=1),[dst_lat.size, 1]),4,axis=1)) ) )
self.nn_wei = dist_wei/np.repeat(np.reshape(np.nansum(dist_wei, axis=1),[dst_lat.size, 1]),4,axis=1)
self.nn_wei[ind | self.msk] = 0.
# Need to identify missing points and throw a warning and set values to zero
mv = np.sum(self.nn_wei,axis=1) == 0
if np.sum(mv) > 1:
print '##WARNING## There are', np.sum(mv), 'missing values, these will be set to ZERO'
else:
print '##WARNING## There is', np.sum(mv), 'missing value, this will be set to ZERO'
def extract_con(self, con):
if con in self.cons:
con_ind = self.cons.index(con)
# Extract the complex amplitude components
nc = Dataset(self.fileID) # pass variable ids to nc
if self.g_type == 't':
vIm = np.ravel(np.concatenate([nc.variables[self.var_Im][con_ind,-2:,:],
nc.variables[self.var_Im][con_ind,:,:],
nc.variables[self.var_Im][con_ind,0:2,:]]))
vRe = np.ravel(np.concatenate([nc.variables[self.var_Re][con_ind,-2:,:],
nc.variables[self.var_Re][con_ind,:,:],
nc.variables[self.var_Re][con_ind,0:2,:]]))
self.harm_Im[con] = np.sum(vIm[self.nn_id]*self.nn_wei,axis=1)
self.harm_Re[con] = np.sum(vRe[self.nn_id]*self.nn_wei,axis=1)
else:
uIm = np.concatenate([nc.variables[self.var_Im][con_ind,-2:,:],
nc.variables[self.var_Im][con_ind,:,:],
nc.variables[self.var_Im][con_ind,0:3,:]])
uRe = np.concatenate([nc.variables[self.var_Re][con_ind,-2:,:],
nc.variables[self.var_Re][con_ind,:,:],
nc.variables[self.var_Re][con_ind,0:3,:]])
vIm = np.concatenate([nc.variables[self.var_Im2][con_ind,-2:,:],
nc.variables[self.var_Im2][con_ind,:,:],
nc.variables[self.var_Im2][con_ind,0:2,:]])
vRe = np.concatenate([nc.variables[self.var_Re2][con_ind,-2:,:],
nc.variables[self.var_Re2][con_ind,:,:],
nc.variables[self.var_Re2][con_ind,0:2,:]])
# Deal with north pole. NB in TPXO7.2 data U and Z at 90N have different
# values for each Longitude value! Plus there's something odd with the
# hu and hv depths not being the min of surrounding T-grid depths
# TODO remove hardwired 722 index point and make generic
vIm = np.concatenate([vIm[:,:], np.concatenate([vIm[722:,-1],vIm[:722,-1]])[:,np.newaxis]],axis=1)
vRe = np.concatenate([vRe[:,:], np.concatenate([vRe[722:,-1],vRe[:722,-1]])[:,np.newaxis]],axis=1)
# Average U and V onto the T-grid
uIm = np.ravel((uIm[:-1,:] + uIm[1:,:])/2)
uRe = np.ravel((uRe[:-1,:] + uRe[1:,:])/2)
vIm = np.ravel((vIm[:,:-1] + vIm[:,1:])/2)
vRe = np.ravel((vRe[:,:-1] + vRe[:,1:])/2)
if self.key_tr: # We convert to velocity using tidal model bathymetry
harm_Im = np.sum(uIm[self.nn_id]*self.nn_wei,axis=1)/np.sum(self.bat*self.nn_wei,axis=1)
harm_Re = np.sum(uRe[self.nn_id]*self.nn_wei,axis=1)/np.sum(self.bat*self.nn_wei,axis=1)
harm_Im2 = np.sum(vIm[self.nn_id]*self.nn_wei,axis=1)/np.sum(self.bat*self.nn_wei,axis=1)
harm_Re2 = np.sum(vRe[self.nn_id]*self.nn_wei,axis=1)/np.sum(self.bat*self.nn_wei,axis=1)
else: # We convert to velocity using the regional model bathymetry
harm_Im = np.sum(uIm[self.nn_id]*self.nn_wei,axis=1)/self.dst_dep
harm_Re = np.sum(uRe[self.nn_id]*self.nn_wei,axis=1)/self.dst_dep
harm_Im2 = np.sum(vIm[self.nn_id]*self.nn_wei,axis=1)/self.dst_dep
harm_Re2 = np.sum(vRe[self.nn_id]*self.nn_wei,axis=1)/self.dst_dep
# Rotate vectors
self.harm_Im[con] = rot_rep(harm_Im, harm_Im2, self.g_type,
'en to %s' %self.rot_dir, self.gcos, self.gsin)
self.harm_Re[con] = rot_rep(harm_Re, harm_Re2, self.g_type,
'en to %s' %self.rot_dir, self.gcos, self.gsin)
self.harm_Im[con][self.msk]=0.
self.harm_Re[con][self.msk]=0.
nc.close()
else:
# throw some warning
print '##WARNING## Missing constituent values will be set to ZERO'
self.harm_Im[con] = np.zeros(self.nn_id[:,0].size)
self.harm_Re[con] = np.zeros(self.nn_id[:,0].size)
...@@ -7,7 +7,7 @@ Module to extract constituents for the input grid mapped onto output grid ...@@ -7,7 +7,7 @@ Module to extract constituents for the input grid mapped onto output grid
# pylint: disable=E1103 # pylint: disable=E1103
# pylint: disable=no-name-in-module # pylint: disable=no-name-in-module
import copy import copy
import tpxo_extract_HC from . import tpxo_extract_HC
import numpy as np import numpy as np
from netCDF4 import Dataset from netCDF4 import Dataset
from pynemo import nemo_bdy_grid_angle from pynemo import nemo_bdy_grid_angle
...@@ -263,7 +263,7 @@ def constituents_index(constituents, inputcons): ...@@ -263,7 +263,7 @@ def constituents_index(constituents, inputcons):
""" """
retindx = np.zeros(len(inputcons)) retindx = np.zeros(len(inputcons))
count = 0 count = 0
for value in inputcons.values(): for value in list(inputcons.values()):
const_name = value.replace("'", "").lower() # force inputcons entries to lowercase const_name = value.replace("'", "").lower() # force inputcons entries to lowercase
retindx[count] = [x.lower() for x in constituents].index(const_name) # force constituents to lowercase retindx[count] = [x.lower() for x in constituents].index(const_name) # force constituents to lowercase
count = count+1 count = count+1
......
...@@ -52,18 +52,11 @@ class TpxoExtract(object): ...@@ -52,18 +52,11 @@ class TpxoExtract(object):
self.cons = [] self.cons = []
for ncon in range(self.height_dataset.variables['con'].shape[0]): for ncon in range(self.height_dataset.variables['con'].shape[0]):
self.cons.append(self.height_dataset.variables['con'][ncon, :].tostring().strip()) self.cons.append(self.height_dataset.variables['con'][ncon, :].tostring().strip())
elif tide_model == 'FES': elif tide_model == 'FES':
constituents = ['2N2','EPS2','J1','K1','K2','L2','LA2','M2','M3','M4','M6','M8','MF','MKS2','MM','MN4','MS4','MSF','MSQM','MTM','MU2','N2','N4','NU2','O1','P1','Q1','R2','S1','S2','S4','SA','SSA','T2'] constituents = ['2N2','EPS2','J1','K1','K2','L2','LA2','M2','M3','M4','M6','M8','MF','MKS2','MM','MN4','MS4','MSF','MSQM','MTM','MU2','N2','N4','NU2','O1','P1','Q1','R2','S1','S2','S4','SA','SSA','T2']
print 'did not actually code stuff for FES in this routine. Though that would be ideal. Instead put it in fes_extract_HC.py' print('did not actually code stuff for FES in this routine. Though that would be ideal. Instead put it in fes_extract_HC.py')
else: else:
print 'Don''t know that tide model' print('Don''t know that tide model')
# Wrap coordinates in longitude if the domain is global # Wrap coordinates in longitude if the domain is global
glob = 0 glob = 0
...@@ -130,7 +123,7 @@ class TpxoExtract(object): ...@@ -130,7 +123,7 @@ class TpxoExtract(object):
VRe_name, VIm_name, lon_v_name, lat_v_name, VRe_name, VIm_name, lon_v_name, lat_v_name,
lon, lat, depth, maskname=mv_name) lon, lat, depth, maskname=mv_name)
else: else:
print 'Unknown grid_type' print('Unknown grid_type')
return return
def interpolate_constituents(self, nc_dataset, real_var_name, img_var_name, lon_var_name, def interpolate_constituents(self, nc_dataset, real_var_name, img_var_name, lon_var_name,
...@@ -221,7 +214,7 @@ def bilinear_interpolation(lon, lat, data, lon_new, lat_new): ...@@ -221,7 +214,7 @@ def bilinear_interpolation(lon, lat, data, lon_new, lat_new):
# n = lon.size # n = lon.size
# m = lat.size # m = lat.size
if lon.size != data.shape[0] or lat.size != data.shape[1]: if lon.size != data.shape[0] or lat.size != data.shape[1]:
print 'Check Dimensions' print('Check Dimensions')
return np.NaN return np.NaN
if glob == 1: if glob == 1:
lon = np.concatenate(([lon[0] - 2 * lon_resolution, lon[0] - lon_resolution, ], lon = np.concatenate(([lon[0] - 2 * lon_resolution, lon[0] - lon_resolution, ],
......
...@@ -153,7 +153,7 @@ def polcoms_select_domain(bathy, lat, lon, roi, dr): ...@@ -153,7 +153,7 @@ def polcoms_select_domain(bathy, lat, lon, roi, dr):
lon_ob = np.ravel(lon,order='F')[np.ravel(ob,order='F')] lon_ob = np.ravel(lon,order='F')[np.ravel(ob,order='F')]
print lat_ob, lon_ob print(lat_ob, lon_ob)
len_lat = len(lat[:,0]) len_lat = len(lat[:,0])
len_lon = len(lon[0,:]) len_lon = len(lon[0,:])
lat_lon_index = np.nonzero( np.logical_and(lat == lat_ob[0], lon == lon_ob[0])) lat_lon_index = np.nonzero( np.logical_and(lat == lat_ob[0], lon == lon_ob[0]))
...@@ -172,19 +172,19 @@ def polcoms_select_domain(bathy, lat, lon, roi, dr): ...@@ -172,19 +172,19 @@ def polcoms_select_domain(bathy, lat, lon, roi, dr):
i_e = i_0 + 2 i_e = i_0 + 2
lat_slice = slice(max(lat_lon_index[0],0),min(lat_lon_index[0]+1+1,len_lat)) lat_slice = slice(max(lat_lon_index[0],0),min(lat_lon_index[0]+1+1,len_lat))
lon_slice = slice(max(lat_lon_index[1],0),min(lat_lon_index[1]+1+1,len_lon)) lon_slice = slice(max(lat_lon_index[1],0),min(lat_lon_index[1]+1+1,len_lon))
print 'method2', lon_slice, lat_slice print('method2', lon_slice, lat_slice)
lat_slice = slice(j_0,j_e) lat_slice = slice(j_0,j_e)
lon_slice = slice(i_0,i_e) lon_slice = slice(i_0,i_e)
print 'method1', lon_slice, lat_slice print('method1', lon_slice, lat_slice)
lat_pts = lat[lat_slice, lon_slice] lat_pts = lat[lat_slice, lon_slice]
lon_pts = lon[lat_slice, lon_slice] lon_pts = lon[lat_slice, lon_slice]
print lat_pts, lon_pts print(lat_pts, lon_pts)
print lat_lon_index[0], lat_lon_index[1] print(lat_lon_index[0], lat_lon_index[1])
print len_lon, len_lat, lat_lon_index[0], lat_lon_index[1] print(len_lon, len_lat, lat_lon_index[0], lat_lon_index[1])
dy,py=seawater.dist(lat_pts[:,0], lon_pts[:,0]) dy,py=seawater.dist(lat_pts[:,0], lon_pts[:,0])
dx,px=seawater.dist(lat_pts[0,:], lon_pts[0,:]) dx,px=seawater.dist(lat_pts[0,:], lon_pts[0,:])
r = np.rint(np.ceil(dr/np.amax([dx,dy]))) r = np.rint(np.ceil(dr/np.amax([dx,dy])))
print dx, dy, r print(dx, dy, r)
lat_slice = slice(max(lat_lon_index[0]-r,0),min(lat_lon_index[0]+r+1,len_lat)) lat_slice = slice(max(lat_lon_index[0]-r,0),min(lat_lon_index[0]+r+1,len_lat))
lon_slice = slice(max(lat_lon_index[1]-r,0),min(lat_lon_index[1]+r+1,len_lon)) lon_slice = slice(max(lat_lon_index[1]-r,0),min(lat_lon_index[1]+r+1,len_lon))
lat_pts = lat[lat_slice, lon_slice] lat_pts = lat[lat_slice, lon_slice]
...@@ -194,8 +194,8 @@ def polcoms_select_domain(bathy, lat, lon, roi, dr): ...@@ -194,8 +194,8 @@ def polcoms_select_domain(bathy, lat, lon, roi, dr):
lon_pts = np.ravel(lon_pts) lon_pts = np.ravel(lon_pts)
# NOTE: seawater package calculates the distance from point to the next point in the array # NOTE: seawater package calculates the distance from point to the next point in the array
# that is the reason to insert reference point before every point # that is the reason to insert reference point before every point
lat_pts = np.insert(lat_pts,range(0,len(lat_pts)), lat_ob[idx]) lat_pts = np.insert(lat_pts,list(range(0,len(lat_pts))), lat_ob[idx])
lon_pts = np.insert(lon_pts,range(0,len(lon_pts)), lon_ob[idx]) lon_pts = np.insert(lon_pts,list(range(0,len(lon_pts))), lon_ob[idx])
distance_pts = seawater.dist(lat_pts, lon_pts) distance_pts = seawater.dist(lat_pts, lon_pts)
#distances repeat themselves so only pick every alternative distance #distances repeat themselves so only pick every alternative distance
distance_pts = distance_pts[0][::2] distance_pts = distance_pts[0][::2]
...@@ -230,19 +230,19 @@ def polcoms_select_domain(bathy, lat, lon, roi, dr): ...@@ -230,19 +230,19 @@ def polcoms_select_domain(bathy, lat, lon, roi, dr):
i_e = i_0 + 2 i_e = i_0 + 2
lat_slice = slice(max(lat_lon_index[0],0),min(lat_lon_index[0]+1+1,len_lat)) lat_slice = slice(max(lat_lon_index[0],0),min(lat_lon_index[0]+1+1,len_lat))
lon_slice = slice(max(lat_lon_index[1],0),min(lat_lon_index[1]+1+1,len_lon)) lon_slice = slice(max(lat_lon_index[1],0),min(lat_lon_index[1]+1+1,len_lon))
print 'method2', lon_slice, lat_slice print('method2', lon_slice, lat_slice)
lat_slice = slice(j_0,j_e) lat_slice = slice(j_0,j_e)
lon_slice = slice(i_0,i_e) lon_slice = slice(i_0,i_e)
print 'method1', lon_slice, lat_slice print('method1', lon_slice, lat_slice)
lat_pts = lat[lat_slice, lon_slice] lat_pts = lat[lat_slice, lon_slice]
lon_pts = lon[lat_slice, lon_slice] lon_pts = lon[lat_slice, lon_slice]
print lat_pts, lon_pts print(lat_pts, lon_pts)
print lat_lon_index[0], lat_lon_index[1] print(lat_lon_index[0], lat_lon_index[1])
print len_lon, len_lat, lat_lon_index[0], lat_lon_index[1] print(len_lon, len_lat, lat_lon_index[0], lat_lon_index[1])
dy,py=seawater.dist(lat_pts[:,0], lon_pts[:,0]) dy,py=seawater.dist(lat_pts[:,0], lon_pts[:,0])
dx,px=seawater.dist(lat_pts[0,:], lon_pts[0,:]) dx,px=seawater.dist(lat_pts[0,:], lon_pts[0,:])
r = np.rint(np.ceil(dr/np.amax([dx,dy]))) r = np.rint(np.ceil(dr/np.amax([dx,dy])))
print dx, dy, r print(dx, dy, r)
lat_slice = slice(max(lat_lon_index[0]-r,0),min(lat_lon_index[0]+r+1,len_lat)) lat_slice = slice(max(lat_lon_index[0]-r,0),min(lat_lon_index[0]+r+1,len_lat))
lon_slice = slice(max(lat_lon_index[1]-r,0),min(lat_lon_index[1]+r+1,len_lon)) lon_slice = slice(max(lat_lon_index[1]-r,0),min(lat_lon_index[1]+r+1,len_lon))
lat_pts = lat[lat_slice, lon_slice] lat_pts = lat[lat_slice, lon_slice]
...@@ -252,8 +252,8 @@ def polcoms_select_domain(bathy, lat, lon, roi, dr): ...@@ -252,8 +252,8 @@ def polcoms_select_domain(bathy, lat, lon, roi, dr):
lon_pts = np.ravel(lon_pts) lon_pts = np.ravel(lon_pts)
# NOTE: seawater package calculates the distance from point to the next point in the array # NOTE: seawater package calculates the distance from point to the next point in the array
# that is the reason to insert reference point before every point # that is the reason to insert reference point before every point
lat_pts = np.insert(lat_pts,range(0,len(lat_pts)), lat_lb[idx]) lat_pts = np.insert(lat_pts,list(range(0,len(lat_pts))), lat_lb[idx])
lon_pts = np.insert(lon_pts,range(0,len(lon_pts)), lon_lb[idx]) lon_pts = np.insert(lon_pts,list(range(0,len(lon_pts))), lon_lb[idx])
distance_pts = seawater.dist(lat_pts, lon_pts) distance_pts = seawater.dist(lat_pts, lon_pts)
#distances repeat themselves so only pick every alternative distance #distances repeat themselves so only pick every alternative distance
distance_pts = distance_pts[0][::2] distance_pts = distance_pts[0][::2]
......
...@@ -74,7 +74,7 @@ def bdy_sections(nbidta,nbjdta,nbrdta,rw): ...@@ -74,7 +74,7 @@ def bdy_sections(nbidta,nbjdta,nbrdta,rw):
count = 0 count = 0
flag = 0 flag = 0
mark = 0 mark = 0
source_tree = sp.cKDTree(zip(outer_rim_i, outer_rim_j)) source_tree = sp.cKDTree(list(zip(outer_rim_i, outer_rim_j)))
id_order = np.ones((nbdy,), dtype=np.int)*source_tree.n id_order = np.ones((nbdy,), dtype=np.int)*source_tree.n
id_order[count] = 0 # use index 0 as the starting point id_order[count] = 0 # use index 0 as the starting point
count += 1 count += 1
...@@ -85,8 +85,8 @@ def bdy_sections(nbidta,nbjdta,nbrdta,rw): ...@@ -85,8 +85,8 @@ def bdy_sections(nbidta,nbjdta,nbrdta,rw):
while count <= nbdy: while count <= nbdy:
lcl_pt = zip([outer_rim_i[id_order[count-1]]], lcl_pt = list(zip([outer_rim_i[id_order[count-1]]],
[outer_rim_j[id_order[count-1]]]) [outer_rim_j[id_order[count-1]]]))
junk, an_id = source_tree.query(lcl_pt, k=3, distance_upper_bound=1.1) junk, an_id = source_tree.query(lcl_pt, k=3, distance_upper_bound=1.1)
if an_id[0,1] in id_order: if an_id[0,1] in id_order:
...@@ -128,3 +128,28 @@ def bdy_transport(): ...@@ -128,3 +128,28 @@ def bdy_transport():
""" """
raise NotImplementedError raise NotImplementedError
def dist(self, x, y):
"""
Return the distance between two points.
"""
d = x-y
return np.sqrt(np.dot(d, d))
def dist_point_to_segment(p, s0, s1):
"""
Get the distance of a point to a segment.
*p*, *s0*, *s1* are *xy* sequences
This algorithm from
http://geomalgorithms.com/a02-_lines.html
"""
v = s1 - s0
w = p - s0
c1 = np.dot(w, v)
if c1 <= 0:
return dist(p, s0)
c2 = np.dot(v, v)
if c2 <= c1:
return dist(p, s1)
b = c1 / c2
pb = s0 + b * v
return dist(p, pb)
...@@ -28,6 +28,7 @@ ln_tra = boundary conditions for T and S ...@@ -28,6 +28,7 @@ ln_tra = boundary conditions for T and S
ln_ice = ice boundary condition ln_ice = ice boundary condition
nn_rimwidth = width of the relaxation zone nn_rimwidth = width of the relaxation zone
ln_tide = if .true. : produce bdy tidal conditions ln_tide = if .true. : produce bdy tidal conditions
sn_tide_model = 'FES'
clname(0) = constituent name clname(0) = constituent name
ln_trans = ln_trans =
nn_year_000 = year start nn_year_000 = year start
......
[bdist_wheel] [bdist_wheel]
python-tag = py27 python-tag = py37
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment