fixed phase issue on FES tide boundary generation

pynemo/tests/nemo_tide_test.py

@@ -46,7 +46,7 @@ logging.basicConfig(filename='nrct.log', level=logging.INFO)
 
 # TODO: add TPXO read and subset functionality currently only uses FES as "truth"
 
-def main(bdy_file='inputs/namelist_cmems.bdy',amplitude_threshold = 0.25,phase_threshold=10.00,model_res=1/16,model='fes'):
+def main(bdy_file='inputs/namelist_cmems.bdy',amplitude_threshold = 0.25,phase_threshold=15.00,model_res=1/16,model='fes'):
     logger.info('============================================')
     logger.info('Start Tide Test Logging: ' + time.asctime())
     logger.info('============================================')

pynemo/tide/fes_extract_HC.py

@@ -21,16 +21,16 @@ class HcExtract(object):
         tide_model = 'fes'
         if tide_model == 'fes':  # Define stuff to generalise Tide model
 
-            hRe_name = 'amp_stack'
-            hIm_name = 'pha_stack'
+            hRe_name = 'hRe'
+            hIm_name = 'hIm'
             lon_z_name = 'lon_z'
             lat_z_name = 'lat_z'
-            URe_name = 'amp_stack'
-            UIm_name = 'pha_stack'
+            URe_name = 'URe'
+            UIm_name = 'UIm'
             lon_u_name = 'lon_u'
             lat_u_name = 'lat_u'
-            VRe_name = 'amp_stack'
-            VIm_name = 'pha_stack'
+            VRe_name = 'VRe'
+            VIm_name = 'VIm'
             lon_v_name = 'lon_v'
             lat_v_name = 'lat_v'
             mz_name = 'mask_z'
@@ -60,18 +60,18 @@ class HcExtract(object):
             self.mask_dataset[mz_name] = mask_z
 
             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] = 1
+            mask_z[mask_z != 18446744073709551616.00000] = 1
             mask_u[mask_u == 18446744073709551616.00000] = 0
             self.mask_dataset[mu_name] = mask_u
 
             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] = 1
+            mask_z[mask_z != 18446744073709551616.00000] = 1
             mask_v[mask_v == 18446744073709551616.00000] = 0
             self.mask_dataset[mv_name] = mask_v
 
             #read and convert the height_dataset file to complex and store in dicts
-            amp_stack = []
-            pha_stack = []
+            hRe = []
+            hIm = []
             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)):
@@ -82,17 +82,15 @@ class HcExtract(object):
                 amp = amp/100.00
                 phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_Z.nc').variables['phase'][:])))
                 phase[phase == 18446744073709551616.00000] = 0
-                # convert from -180-180 to 0-360
-                phase[phase < 0.0] = phase[phase < 0.0] + 360.0
-                amp_stack.append(amp)
-                pha_stack.append(phase)
-            amp_stack = np.stack(amp_stack)
-            pha_stack = np.stack(pha_stack)
-            self.height_dataset = [lon_z,lat_z,amp_stack,pha_stack]
+                hRe.append(amp*np.cos(phase*(np.pi/180)))
+                hIm.append(-amp*np.sin(phase*(np.pi/180)))
+            hRe = np.stack(hRe)
+            hIm = np.stack(hIm)
+            self.height_dataset = [lon_z,lat_z,hRe,hIm]
 
             #read and convert the velocity_dataset files to complex
-            amp_stack = []
-            pha_stack = []
+            URe = []
+            UIm = []
             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)):
@@ -103,15 +101,14 @@ class HcExtract(object):
                 amp = amp/100.00
                 phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ug'][:])))
                 phase[phase == 18446744073709551616.00000] = 0
-                phase[phase < 0.0] = phase[phase < 0.0] + 360.0
-                amp_stack.append(amp)
-                pha_stack.append(phase)
-            amp_stack = np.stack(amp_stack)
-            pha_stack = np.stack(pha_stack)
-            self.Uvelocity_dataset = [lon_u,lat_u,amp_stack,pha_stack]
-
-            amp_stack = []
-            pha_stack = []
+                URe.append(amp*np.cos(phase*(np.pi/180)))
+                UIm.append(-amp*np.sin(phase*(np.pi/180)))
+            URe = np.stack(URe)
+            UIm = np.stack(UIm)
+            self.Uvelocity_dataset = [lon_u,lat_u,URe,UIm]
+
+            VRe = []
+            VIm = []
             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)):
@@ -122,12 +119,11 @@ class HcExtract(object):
                 amp = amp/100.00
                 phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Vg'][:])))
                 phase[phase == 18446744073709551616.00000] = 0
-                phase[phase < 0.0] = phase[phase < 0.0] + 360.0
-                amp_stack.append(amp)
-                pha_stack.append(phase)
-            amp_stack = np.stack(amp_stack)
-            pha_stack = np.stack(pha_stack)
-            self.Vvelocity_dataset = [lon_v,lat_v,amp_stack,pha_stack]
+                VRe.append(amp*np.cos(phase*(np.pi/180)))
+                VIm.append(-amp*np.sin(phase*(np.pi/180)))
+            VRe = np.stack(VRe)
+            VIm = np.stack(VIm)
+            self.Vvelocity_dataset = [lon_v,lat_v,VRe,VIm]
 
             # open grid variables these are resampled TPXO parameters so may not work correctly.
             self.grid = Dataset(settings['tide_fes']+'grid_fes.nc')
@@ -209,10 +205,9 @@ class HcExtract(object):
         amp = np.zeros((len(nc_dataset[2]), lon.shape[0]))
         gph = np.zeros((len(nc_dataset[2]), lon.shape[0]))
 
-        data_amp = np.array(np.ravel(nc_dataset[2]))
-        data_pha = np.array(np.ravel(nc_dataset[3]))
-        data_amp = data_amp.reshape(nc_dataset[2].shape)
-        data_pha = data_pha.reshape(nc_dataset[2].shape)
+        data = np.array(np.ravel(nc_dataset[2]), dtype=complex)
+        data.imag = np.array(np.ravel(nc_dataset[3]))
+        data = data.reshape(nc_dataset[2].shape)
         #data = data.reshape(1,nc_dataset['M2'][real_var_name].shape)
         # data[data==0] = np.NaN
 
@@ -243,26 +238,29 @@ class HcExtract(object):
         #interpolate the mask values
         maskedpoints = interpolate.interpn((x_values, y_values), mask, lonlat)
 
-        data_temp_amp = np.zeros((data_amp.shape[0], lon.shape[0], 1, ))
-        data_temp_pha = np.zeros((data_pha.shape[0], lon.shape[0], 1,))
-        for cons_index in range(data_amp.shape[0]):
+        data_temp = np.zeros((data.shape[0], lon.shape[0], 2, ))
+        for cons_index in range(data.shape[0]):
             #interpolate real values
-            data_temp_amp[cons_index, :, 0] = interpolate_data(x_values, y_values,
-                                                                data_amp[cons_index, :, :],
+            data_temp[cons_index, :, 0] = interpolate_data(x_values, y_values,
+                                                                data[cons_index, :, :].real,
                                                                 maskedpoints, lonlat)
             #interpolate imag values
-            data_temp_pha[cons_index, :, 0] = interpolate_data(x_values, y_values,
-                                                                data_pha[cons_index, :, :],
+            data_temp[cons_index, :, 1] = interpolate_data(x_values, y_values,
+                                                                data[cons_index, :, :].imag,
                                                                 maskedpoints, lonlat)
 
             #for velocity_dataset values
             if height_data is not None:
-                data_temp_amp[cons_index, :, 0] = data_temp_amp[cons_index, :, 0]/height_data*100
-                data_temp_pha[cons_index, :, 0] = data_temp_pha[cons_index, :, 0]/height_data*100
+                data_temp[cons_index, :, 0] = data_temp[cons_index, :, 0]/height_data*100
+                data_temp[cons_index, :, 1] = data_temp[cons_index, :, 1]/height_data*100
 
-            amp[cons_index, :] = data_temp_amp[cons_index, :, 0]
-            gph[cons_index, :] = data_temp_pha[cons_index, :, 0]
+            zcomplex = np.array(data_temp[cons_index, :, 0], dtype=complex)
+            zcomplex.imag = data_temp[cons_index, :, 1]
 
+            amp[cons_index, :] = np.absolute(zcomplex)
+            gph[cons_index, :] = np.arctan2(-1*zcomplex.imag, zcomplex.real)
+        gph = gph*180.0/np.pi
+        gph[gph < 0] = gph[gph < 0]+360.0
         return amp, gph
 
 def interpolate_data(lon, lat, data, mask, lonlat):