Skip to content

GitLab

Commit 13915375 authored by thopri's avatar thopri
Browse files

fixed large spikes in FES cosine and sin output

parent a1d80702
Hide whitespace changes
Inline Side-by-side
Showing with 38 additions and 135 deletions
inputs/namelist_cmems.bdy
View file @ 13915375
......@@ -52,7 +52,7 @@
! CMEMS Data Source Configuration
!------------------------------------------------------------------------------
ln_use_cmems = .true.
ln_download_cmems = .true.
ln_download_cmems = .false.
sn_cmems_dir = '/Users/thopri/Projects/PyNEMO/inputs/' ! where to download CMEMS input files (static and variable)
ln_download_static = .false.
ln_subset_static = .false.
......@@ -103,7 +103,7 @@
!------------------------------------------------------------------------------
! unstructured open boundaries tidal parameters
!------------------------------------------------------------------------------
ln_tide = .false. ! =T : produce bdy tidal conditions
ln_tide = .true. ! =T : produce bdy tidal conditions
sn_tide_model = 'fes' ! Name of tidal model (fes|tpxo)
clname(1) = 'M2' ! constituent name
clname(2) = 'S2'
......
pynemo/tide/fes_extract_HC.py
View file @ 13915375
......@@ -53,26 +53,20 @@ class HcExtract(object):
lon_resolution = lon_z[1] - lon_z[0]
data_in_km = 0 # added to maintain the reference to matlab tmd code
# extract example amplitude grid for Z, U and V and change NaNs to 0 (for land) and other values to 1 (for water)
mask_z = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_Z.nc').variables['amplitude'][:]))
where_are_NaNs = np.isnan(mask_z)
where_no_NaNs = np.invert(np.isnan(mask_z))
mask_z[where_no_NaNs] = 1
mask_z[where_are_NaNs] = 0
# extract example amplitude grid for Z, U and V and change NaNs to 1 (for land) and other values to 0 (for water)
mask_z = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_Z.nc').variables['amplitude'][:])))
mask_z[mask_z == 18446744073709551616.00000] = 0
mask_z[mask_z != 0] = 1
self.mask_dataset[mz_name] = mask_z
mask_u = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_U.nc').variables['Ua'][:]))
where_are_NaNs = np.isnan(mask_u)
where_no_NaNs = np.invert(np.isnan(mask_u))
mask_u[where_no_NaNs] = 1
mask_u[where_are_NaNs] = 0
mask_u = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_U.nc').variables['Ua'][:])))
mask_u[mask_u == 18446744073709551616.00000] = 0
mask_u[mask_u != 0] = 1
self.mask_dataset[mu_name] = mask_u
mask_v = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_V.nc').variables['Va'][:]))
where_are_NaNs = np.isnan(mask_v)
where_no_NaNs = np.invert(np.isnan(mask_v))
mask_v[where_no_NaNs] = 1
mask_v[where_are_NaNs] = 0
mask_v = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_V.nc').variables['Va'][:])))
mask_v[mask_v == 18446744073709551616.00000] = 0
mask_v[mask_v != 0] = 1
self.mask_dataset[mv_name] = mask_v
......@@ -82,8 +76,10 @@ class HcExtract(object):
lat_z = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_Z.nc').variables['lat'][:])
lon_z = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_Z.nc').variables['lon'][:])
for ncon in range(len(constituents)):
amp = np.array(np.rot90(Dataset(settings['tide_fes']+str(constituents[ncon])+'_Z.nc').variables['amplitude'][:]))
phase = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_Z.nc').variables['phase'][:]))
amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+str(constituents[ncon])+'_Z.nc').variables['amplitude'][:])))
amp[amp == 18446744073709551616.00000] = 0
amp = amp/100.00
phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_Z.nc').variables['phase'][:])))
hRe.append(amp*np.sin(phase))
hIm.append(amp*np.cos(phase))
hRe = np.stack(hRe)
......@@ -96,8 +92,10 @@ class HcExtract(object):
lat_u = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_U.nc').variables['lat'][:])
lon_u = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_U.nc').variables['lon'][:])
for ncon in range(len(constituents)):
amp = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ua'][:]))
phase = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ug'][:]))
amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ua'][:])))
amp[amp == 18446744073709551616.00000] = 0
amp = amp/100.00
phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ug'][:])))
URe.append(amp*np.sin(phase))
UIm.append(amp*np.cos(phase))
URe = np.stack(URe)
......@@ -109,8 +107,10 @@ class HcExtract(object):
lat_v = np.array(Dataset(settings['tide_fes'] + constituents[ncon] + '_V.nc').variables['lat'][:])
lon_v = np.array(Dataset(settings['tide_fes'] + constituents[ncon] + '_V.nc').variables['lon'][:])
for ncon in range(len(constituents)):
amp = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Va'][:]))
phase = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Vg'][:]))
amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Va'][:])))
amp[amp == 18446744073709551616.00000] = 0
amp = amp/100.00
phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Vg'][:])))
VRe.append(amp*np.sin(phase))
VIm.append(amp*np.cos(phase))
VRe = np.stack(VRe)
......
test_scripts/bdy_var_plot.py
View file @ 13915375
......@@ -61,13 +61,13 @@ def nemo_bdy_order(fname):
id_order[0,] = 0
flag = False
mark = 0
source_tree = sp.cKDTree(zip(nbi_tmp, nbj_tmp), balanced_tree=False, compact_nodes=False)
source_tree = sp.cKDTree(list(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]]),
nn_dist, nn_id = source_tree.query(list(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?
......@@ -207,19 +207,19 @@ def plot_bdy(fname, bdy_ind, bdy_dst, bdy_brk, varnam, t, rw):
# create a pseudo bathymetry from the depth data
bathy = np.zeros_like(coords)
mbath = np.sum(dta[n].mask == 0, axis=0)
#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]
#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)
#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].add_patch(bathy_patch)
ax[n].set_title('BDY points along section: ' + str(n))
patches = []
colors = []
......@@ -248,17 +248,17 @@ def plot_bdy(fname, bdy_ind, bdy_dst, bdy_brk, varnam, t, rw):
# 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))
#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].set_ylim((0, np.max(bathy)))
ax[n].invert_yaxis()
return f
fname = '/Users/thopri/Projects/PyNEMO/outputs/NNA_R12_bdyT_y2017m11.nc'
fname = '/Users/thopri/Projects/PyNEMO/outputs/NNA_R12_bdyT_y2017m01.nc'
print(fname)
ind, dst, brk = nemo_bdy_order(fname)
f = plot_bdy(fname, ind, dst, brk, 'thetao', 0, 0)
......
test_scripts/test_cmems.py deleted 100644 → 0
View file @ a1d80702
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Thu Oct 31 15:49:52 2019
Test script to test download CMEMS data function in PyNEMO. To be incorprated
into profile.py
@author: thopri
"""
from pynemo import nemo_bdy_dl_cmems as dl_cmems
from pynemo import nemo_bdy_setup as setup
from calendar import monthrange
import sys
import time
import threading
Setup = setup.Setup('/Users/thopri/Projects/PyNEMO/inputs/namelist_cmems.bdy') # default settings file
settings = Setup.settings
if settings['use_cmems'] == True:
if settings['year_end'] - settings['year_000'] > 0:
date_min = settings['year_000']+'-01-01'
date_max = settings['year_end']+'-12-31'
days_mth = monthrange(settings['year_end'],settings['month_end'])
date_min = str(settings['year_000'])+'-'+str(settings['month_000'])+'-01'
date_max = str(settings['year_end'])+'-'+str(settings['month_end'])+'-'+str(days_mth[1])
cmes_config= {
'ini_config_template' : settings['cmes_config_template'],
'user' : settings['cmes_usr'],
'pwd' : settings['cmes_pwd'],
'motu_server' : settings['motu_server'],
'service_id' : settings['cmes_model'],
'product_id' : settings['cmes_product'],
'date_min' : date_min,
'date_max' : date_max,
'latitude_min' : settings['latitude_min'],
'latitude_max' : settings['latitude_max'],
'longitude_min' : settings['longitude_min'],
'longitude_max' : settings['longitude_max'],
'depth_min' : settings['depth_min'],
'depth_max' : settings['depth_max'],
'variable' : settings['cmes_variable'],
'src_dir' : settings['src_dir'],
'out_name' : settings['cmes_output'],
'config_out' : settings['cmes_config']
}
chk = dl_cmems.chk_motu()
dl = dl_cmems.request_CMEMS(cmes_config)
class progress_bar_loading(threading.Thread):
def run(self):
global stop
global kill
print 'Loading.... ',
sys.stdout.flush()
i = 0
while stop != True:
if (i%4) == 0:
sys.stdout.write('\b/')
elif (i%4) == 1:
sys.stdout.write('\b-')
elif (i%4) == 2:
sys.stdout.write('\b\\')
elif (i%4) == 3:
sys.stdout.write('\b|')
sys.stdout.flush()
time.sleep(0.2)
i+=1
if kill == True:
print '\b\b\b\b ABORT!',
else:
print '\b\b done!',
kill = False
stop = False
p = progress_bar_loading()
p.start()
try:
#anything you want to run.
time.sleep(1)
stop = True
except KeyboardInterrupt or EOFError:
kill = True
stop = True
\ 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