Skip to content

GitLab

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

Merge pull request #47 from jdha/PyNEMO3 

Py nemo3
parents 083ed633 298e9b48
Hide whitespace changes
Inline Side-by-side
Showing with 122 additions and 321 deletions
pynemo/pynemo_exe.py
View file @ ecde918e
......@@ -5,7 +5,7 @@ Entry for the project
'''
import sys, getopt
import profile
from . import profile
import logging
# Logging set to info
......@@ -20,14 +20,14 @@ def main():
try:
opts, dummy_args = getopt.getopt(sys.argv[1:], "hs:g", ["help","setup=","mask_gui"])
except getopt.GetoptError:
print "usage: pynemo -g -s <namelist.bdy> "
print("usage: pynemo -g -s <namelist.bdy> ")
sys.exit(2)
for opt, arg in opts:
if opt == "-h":
print "usage: pynemo [-g] -s <namelist.bdy> "
print " -g (optional) will open settings editor before extracting the data"
print " -s <bdy filename> file to use"
print("usage: pynemo [-g] -s <namelist.bdy> ")
print(" -g (optional) will open settings editor before extracting the data")
print(" -s <bdy filename> file to use")
sys.exit()
elif opt in ("-s", "--setup"):
setup_file = arg
......@@ -35,7 +35,7 @@ def main():
mask_gui = True
if setup_file == "":
print "usage: pynemo [-g] -s <namelist.bdy> "
print("usage: pynemo [-g] -s <namelist.bdy> ")
sys.exit(2)
#Logger
......@@ -43,7 +43,7 @@ def main():
t0 = time.time()
profile.process_bdy(setup_file, mask_gui)
t1 = time.time()
print "Execution Time: %s" % (t1-t0)
print("Execution Time: %s" % (t1-t0))
if __name__ == "__main__":
main()
pynemo/pynemo_ncml_generator.py
View file @ ecde918e
......@@ -7,8 +7,8 @@ Used for development purposes to display the ncml editor dialog.
@author: Shirley Crompton, UK Science and Technology Facilities Council
'''
import sys
from PyQt4.QtGui import *
from gui import nemo_ncml_generator as ncml_generator
from PyQt5.QtWidgets import *
from .gui import nemo_ncml_generator as ncml_generator
import logging
# Logging set to info
logging.basicConfig(level=logging.INFO)
......@@ -21,4 +21,4 @@ def main():
sys.exit(app.exec_())
if __name__ == '__main__':
main()
\ No newline at end of file
main()
pynemo/pynemo_settings_editor.py
View file @ ecde918e
......@@ -5,10 +5,10 @@ Created on 7 Jan 2015
'''
# pylint: disable=E1103
# pylint: disable=no-name-in-module
from PyQt4 import QtGui
from PyQt5 import QtGui, QtWidgets
from gui.nemo_bdy_input_window import InputWindow
import nemo_bdy_setup
from .gui.nemo_bdy_input_window import InputWindow
from . import nemo_bdy_setup
import sys, getopt
......@@ -16,9 +16,9 @@ def open_settings_window(fname):
""" Main method which starts a Qt application and gives user
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"""
app = QtGui.QApplication(sys.argv)
app = QtWidgets.QApplication(sys.argv)
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')
ex = InputWindow(setup)
......@@ -29,7 +29,7 @@ def open_settings_dialog(setup):
""" 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
but carries on with the execution"""
app = QtGui.QApplication(sys.argv)
app = QtWidgets.QApplication(sys.argv)
ex = InputWindow(setup)
ex.nl_editor.btn_cancel.clicked.connect(app.quit)
return app.exec_(), ex.mpl_widget.mask
......@@ -41,12 +41,12 @@ def main():
try:
opts, dummy_args = getopt.getopt(sys.argv[1:], "hs:", ["help", "setup="])
except getopt.GetoptError:
print "usage: pynemo_settings_editor -s <namelist.bdy> "
print("usage: pynemo_settings_editor -s <namelist.bdy> ")
sys.exit(2)
for opt, arg in opts:
if opt == "-h":
print "usage: pynemo_settings_editor -s <namelist.bdy> "
print("usage: pynemo_settings_editor -s <namelist.bdy> ")
sys.exit()
elif opt in ("-s", "--setup"):
setup_file = arg
......
pynemo/reader/directory.py
View file @ ecde918e
......@@ -5,7 +5,7 @@ This is an abstraction for the data repository
from os import listdir
import numpy as np
from netCDF4 import Dataset
from netCDF4 import netcdftime
from cftime import utime
import copy
import logging
class Reader(object):
......@@ -50,7 +50,7 @@ class Reader(object):
dir_list[i] = self.directory + dir_list[i]
dir_list.sort()
return filter(None, dir_list)
return [_f for _f in dir_list if _f]
def _delta_time_interval(self, time1, time2):
""" Get the difference between the two times in days and hours"""
......@@ -74,7 +74,7 @@ class Reader(object):
x = [filename, index]
group.data_list.append(x)
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))
group.units = varid.units
group.calendar = varid.calendar
......@@ -94,7 +94,7 @@ class Reader(object):
error"""
days = 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],
self.grid_source_data[grid_type].time_counter[1])
days.add(day)
......
pynemo/reader/factory.py
View file @ ecde918e
......@@ -15,13 +15,13 @@ from netCDF4 import Dataset
def GetReader(uri, t_adjust, reader_type=None):
if reader_type is None:
print uri
print(uri)
if uri.endswith(".ncml"):
reader_type = "NcML"
elif os.path.isdir(uri):
reader_type = "Directory"
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
if reader_type == "NcML":
return NcMLReader(uri,t_adjust)
......
pynemo/reader/ncml.py
View file @ ecde918e
......@@ -10,7 +10,7 @@ import string
import logging
import numpy as np
import jnius_config
from netCDF4 import netcdftime
from cftime import utime
ncmlpath, file_name = os.path.split(__file__)
ncmlpath = os.path.join(ncmlpath, "jars", "netcdfAll-4.6.jar")
jnius_config.set_classpath('.',ncmlpath)
......@@ -22,7 +22,7 @@ try:
proxy_port = proxylist[1]
jnius_config.add_options('-Dhttp.proxyHost='+proxy_host,'-Dhttp.proxyPort='+proxy_port)
except:
print "Didn't find a proxy environment variable"
print("Didn't find a proxy environment variable")
NetcdfDataset = None
NcMLReader = None
Section = None
......@@ -37,7 +37,7 @@ try:
Section = autoclass('ucar.ma2.Section')
init_jnius()
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"
class Reader(object):
......@@ -78,7 +78,7 @@ class Reader(object):
grid.time_counter = timevar[:]+t_adjust
grid.date_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]))
def close(self):
......
pynemo/tests/gcoms_break_depth_test.py
View file @ ecde918e
......@@ -23,8 +23,8 @@ class Test(unittest.TestCase):
self.lon,self.lat = np.meshgrid(self.lon, self.lat)
#gcoms_break_depth.gcoms_boundary_masks(self.bathy, -1, 0)
print self.bathy.shape
print self.bathy
print(self.bathy.shape)
print(self.bathy)
pass
......@@ -39,10 +39,10 @@ class Test(unittest.TestCase):
self.bathy[self.bathy>=0] = 0
self.bathy = self.bathy*-1
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.assertEquals(tmp[40,0],1,"Set the break select correctly 40")
self.assertEquals(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[32,0],1,"Set the break select correctly")
self.assertEqual(tmp[40,0],1,"Set the break select correctly 40")
self.assertEqual(tmp[50,0],1,"Set the break select correctly 50")
self.assertEqual(tmp[60,0],1,"Set the break select correctly 60")
def testGcomsBreakDepth(self):
r = 18
......
pynemo/tests/nemo_bdy_msk_c_test.py
View file @ ecde918e
......@@ -11,8 +11,8 @@ class Test(unittest.TestCase):
def testMaskData(self):
mask = Mask('data/grid_C/NNA_R12_bathy_meter_bench.nc')
print mask.data.shape
self.assertEqual(mask.data.shape,(401L,351L),'Mask reading failed')
print(mask.data.shape)
self.assertEqual(mask.data.shape,(401,351),'Mask reading failed')
if __name__ == "__main__":
......
pynemo/tests/nemo_bdy_setup_test.py
View file @ ecde918e
......@@ -45,7 +45,7 @@ class Test(unittest.TestCase):
#test
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
os.remove('test.bdy')
......@@ -64,12 +64,12 @@ class Test(unittest.TestCase):
#test
setup = Setup('test.bdy')
print setup.settings
self.assertEquals(setup.settings['nonambigious'],False,"Didn't recognize logical setting")
self.assertEquals(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.assertEquals(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")
print(setup.settings)
self.assertEqual(setup.settings['nonambigious'],False,"Didn't recognize logical setting")
self.assertEqual(setup.settings['floatval'],10.9,"Didn't recognize rn value in setting")
self.assertEqual(setup.settings['floatval2'],20.9,"Didn't recongnize nn value in setting")
self.assertEqual(setup.settings['stringval'],'coordinates.nc',"Didn't recognize cn string value in setting")
self.assertEqual(setup.settings['stringval2'],'gregorian',"Didn't recognize sn string value in setting")
#delete file
os.remove('test.bdy')
......
pynemo/tests/reader_ncml_test.py
View file @ ecde918e
......@@ -34,7 +34,7 @@ class Test(unittest.TestCase):
testfile = os.path.join(testpath, "testremote.ncml")
sd = Reader(testfile,0).grid
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):
init_jnius()
......@@ -80,11 +80,11 @@ class Test(unittest.TestCase):
testpath, file_name = os.path.split(__file__)
testfile = os.path.join(testpath, "testremote.ncml")
repo = Reader(testfile,0)
self.assertEquals(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(len(repo['t'].time_counter), 8, "Time counter should be 8")
self.assertEqual(repo['t'].time_counter[0], 691416000, "The first time value doesn't match")
repo = Reader(testfile,100)
self.assertEquals(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(len(repo['t'].time_counter), 8, "Time counter should be 8")
self.assertEqual(repo['t'].time_counter[0], 691416100, "The first time value doesn't match")
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testName']
......
pynemo/tide/nemo_bdy_tide.py
View file @ ecde918e
......@@ -77,8 +77,8 @@ class Extract:
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)
source_tree = sp.cKDTree(list(zip(lon, 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,
distance_upper_bound=0.5)
......@@ -106,7 +106,7 @@ class Extract:
# Need to identify missing points and throw a warning and set values to zero
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):
......@@ -150,7 +150,7 @@ class Extract:
else:
# 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_Re[con] = np.zeros(self.nn_id[:,0].size)
......
pynemo/tide/nemo_bdy_tide2.py
View file @ ecde918e
'''
'''
# 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)
Non
\ No newline at end of file
pynemo/tide/nemo_bdy_tide3.py
View file @ ecde918e
......@@ -7,7 +7,7 @@ Module to extract constituents for the input grid mapped onto output grid
# pylint: disable=E1103
# pylint: disable=no-name-in-module
import copy
import tpxo_extract_HC
from . import tpxo_extract_HC
import numpy as np
from netCDF4 import Dataset
from pynemo import nemo_bdy_grid_angle
......@@ -263,7 +263,7 @@ def constituents_index(constituents, inputcons):
"""
retindx = np.zeros(len(inputcons))
count = 0
for value in inputcons.values():
for value in list(inputcons.values()):
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
count = count+1
......
pynemo/tide/tpxo_extract_HC.py
View file @ ecde918e
......@@ -18,21 +18,21 @@ class TpxoExtract(object):
"""initialises the Extract of tide information from the netcdf
Tidal files"""
# Set tide model
tide_model = 'TPXO'
tide_model = 'TPXO'
if tide_model == 'TPXO': # Define stuff to generalise Tide model
hRe_name = 'hRe'
hIm_name = 'hIm'
lon_z_name = 'lon_z'
lat_z_name = 'lat_z'
if tide_model == 'TPXO': # Define stuff to generalise Tide model
hRe_name = 'hRe'
hIm_name = 'hIm'
lon_z_name = 'lon_z'
lat_z_name = 'lat_z'
URe_name = 'URe'
UIm_name = 'UIm'
lon_u_name = 'lon_u'
lat_u_name = 'lat_u'
UIm_name = 'UIm'
lon_u_name = 'lon_u'
lat_u_name = 'lat_u'
VRe_name = 'VRe'
VIm_name = 'VIm'
lon_v_name = 'lon_v'
lat_v_name = 'lat_v'
VIm_name = 'VIm'
lon_v_name = 'lon_v'
lat_v_name = 'lat_v'
mz_name = 'mz'
mu_name = 'mu'
mv_name = 'mv'
......@@ -52,18 +52,11 @@ class TpxoExtract(object):
self.cons = []
for ncon in range(self.height_dataset.variables['con'].shape[0]):
self.cons.append(self.height_dataset.variables['con'][ncon, :].tostring().strip())
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']
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:
print 'Don''t know that tide model'
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']
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:
print('Don''t know that tide model')
# Wrap coordinates in longitude if the domain is global
glob = 0
......@@ -130,7 +123,7 @@ class TpxoExtract(object):
VRe_name, VIm_name, lon_v_name, lat_v_name,
lon, lat, depth, maskname=mv_name)
else:
print 'Unknown grid_type'
print('Unknown grid_type')
return
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):
# n = lon.size
# m = lat.size
if lon.size != data.shape[0] or lat.size != data.shape[1]:
print 'Check Dimensions'
print('Check Dimensions')
return np.NaN
if glob == 1:
lon = np.concatenate(([lon[0] - 2 * lon_resolution, lon[0] - lon_resolution, ],
......
pynemo/utils/gcoms_break_depth.py
View file @ ecde918e
......@@ -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')]
print lat_ob, lon_ob
print(lat_ob, lon_ob)
len_lat = len(lat[:,0])
len_lon = len(lon[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):
i_e = i_0 + 2
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))
print 'method2', lon_slice, lat_slice
print('method2', lon_slice, lat_slice)
lat_slice = slice(j_0,j_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]
lon_pts = lon[lat_slice, lon_slice]
print lat_pts, lon_pts
print lat_lon_index[0], lat_lon_index[1]
print len_lon, len_lat, lat_lon_index[0], lat_lon_index[1]
print(lat_pts, lon_pts)
print(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])
dx,px=seawater.dist(lat_pts[0,:], lon_pts[0,:])
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))
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]
......@@ -194,8 +194,8 @@ def polcoms_select_domain(bathy, lat, lon, roi, dr):
lon_pts = np.ravel(lon_pts)
# 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
lat_pts = np.insert(lat_pts,range(0,len(lat_pts)), lat_ob[idx])
lon_pts = np.insert(lon_pts,range(0,len(lon_pts)), lon_ob[idx])
lat_pts = np.insert(lat_pts,list(range(0,len(lat_pts))), lat_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)
#distances repeat themselves so only pick every alternative distance
distance_pts = distance_pts[0][::2]
......@@ -230,19 +230,19 @@ def polcoms_select_domain(bathy, lat, lon, roi, dr):
i_e = i_0 + 2
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))
print 'method2', lon_slice, lat_slice
print('method2', lon_slice, lat_slice)
lat_slice = slice(j_0,j_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]
lon_pts = lon[lat_slice, lon_slice]
print lat_pts, lon_pts
print lat_lon_index[0], lat_lon_index[1]
print len_lon, len_lat, lat_lon_index[0], lat_lon_index[1]
print(lat_pts, lon_pts)
print(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])
dx,px=seawater.dist(lat_pts[0,:], lon_pts[0,:])
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))
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]
......@@ -252,8 +252,8 @@ def polcoms_select_domain(bathy, lat, lon, roi, dr):
lon_pts = np.ravel(lon_pts)
# 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
lat_pts = np.insert(lat_pts,range(0,len(lat_pts)), lat_lb[idx])
lon_pts = np.insert(lon_pts,range(0,len(lon_pts)), lon_lb[idx])
lat_pts = np.insert(lat_pts,list(range(0,len(lat_pts))), lat_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)
#distances repeat themselves so only pick every alternative distance
distance_pts = distance_pts[0][::2]
......
pynemo/utils/nemo_bdy_lib.py
View file @ ecde918e
......@@ -74,7 +74,7 @@ def bdy_sections(nbidta,nbjdta,nbrdta,rw):
count = 0
flag = 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[count] = 0 # use index 0 as the starting point
count += 1
......@@ -85,8 +85,8 @@ def bdy_sections(nbidta,nbjdta,nbrdta,rw):
while count <= nbdy:
lcl_pt = zip([outer_rim_i[id_order[count-1]]],
[outer_rim_j[id_order[count-1]]])
lcl_pt = list(zip([outer_rim_i[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)
if an_id[0,1] in id_order:
......@@ -127,4 +127,29 @@ def bdy_transport():
Keyword arguments:
"""
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)
pynemo/variable.info
View file @ ecde918e
......@@ -28,6 +28,7 @@ ln_tra = boundary conditions for T and S
ln_ice = ice boundary condition
nn_rimwidth = width of the relaxation zone
ln_tide = if .true. : produce bdy tidal conditions
sn_tide_model = 'FES'
clname(0) = constituent name
ln_trans =
nn_year_000 = year start
......@@ -47,4 +48,4 @@ nn_alpha = Euler rotation angle
nn_beta = Euler rotation angle
nn_gamma = Euler rotation angle
rn_mask_max_depth = Maximum depth to be ignored for the mask
rn_mask_shelfbreak_dist = Distance from the shelf break
\ No newline at end of file
rn_mask_shelfbreak_dist = Distance from the shelf break
setup.cfg
View file @ ecde918e
[bdist_wheel]
python-tag = py27
\ No newline at end of file
python-tag = py37
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