Commit dd951d03 authored by Ricardo Torres's avatar Ricardo Torres 💬

merge resolution

parents 3249fce6 8ef7d57d
......@@ -2465,7 +2465,16 @@ def read_fvcom_mesh(mesh):
Parameters
----------
mesh : str
Full path to the FVCOM unstructured grid file (.dat usually).
Full path to the FVCOM unstructured grid file (.dat usually). The file
name should also contain '_grd' in the name e.g. 'my_file_grd.dat'.
The file contains information about the triangle indicies and the
x, y, z coodinates of the grid nodes.
The file header should be two lines:
Node Number = nnn
Cell Number = eee
Followed by 'element_index, tri1, tri2, tri3' for eee cells
Followed by 'node_index, x, y, z' for nnn nodes
Returns
-------
......@@ -2481,24 +2490,71 @@ def read_fvcom_mesh(mesh):
"""
fileRead = open(mesh, 'r')
# Skip the file header (two lines)
lines = fileRead.readlines()[2:]
# Read the file and header (two lines)
lines = fileRead.readlines()
fileRead.close()
header = lines[:2] # two lines of header
lines = lines[2:]
node_number = 0
cell_number = 0
triangles = []
nodes = []
x = []
y = []
z = []
for line in lines:
ttt = line.strip().split()
if len(ttt) == 5:
if 'Node' in header[0]:
node_number = int(header[0].split('=')[1])
cell_number = int(header[1].split('=')[1])
elif 'Cell' in header[0]:
node_number = int(header[1].split('=')[1])
cell_number = int(header[0].split('=')[1])
if ((node_number == 0) | (cell_number == 0)):
# if header info not present get grid using column numbers
for line in lines:
ttt = line.strip().split()
if len(ttt) == 5:
t1 = int(ttt[1]) - 1
t2 = int(ttt[2]) - 1
t3 = int(ttt[3]) - 1
triangles.append([t1, t2, t3])
elif len(ttt) == 4:
x.append(float(ttt[1]))
y.append(float(ttt[2]))
z.append(float(ttt[3]))
nodes.append(int(ttt[0]))
else:
# Check if triangles or xyz comes first
tri_first = node_number == int(lines[-1].strip().split()[0])
# Get index ranges of the two sections
if tri_first:
tri_st = 0
tri_en = cell_number
xyz_st = cell_number
xyz_en = cell_number + node_number
else:
xyz_st = 0
xyz_en = node_number
tri_st = node_number
tri_en = node_number + cell_number
# Extract data from the two sections
for l in range(tri_st, tri_en):
ttt = lines[l].strip().split()
t1 = int(ttt[1]) - 1
t2 = int(ttt[2]) - 1
t3 = int(ttt[3]) - 1
triangles.append([t1, t2, t3])
elif len(ttt) == 4:
for l in range(xyz_st, xyz_en):
ttt = lines[l].strip().split()
x.append(float(ttt[1]))
y.append(float(ttt[2]))
z.append(float(ttt[3]))
......
......@@ -1540,7 +1540,34 @@ class FileReader(Domain):
def add_river_flow(self, river_nc_file, river_nml_file):
"""
TODO: docstring.
Read in the river forcing data.
The river forcing is in two files the netcdf with time varying flux,
temperature, salinity for each river node.
The nml text file contains information about which node the river is
attached to.
Parameters
----------
river_nc_file : str
Path to the netcdf file.
river_nml_file : str
Path to the text nml file.
Returns
-------
Populates:
self.river.time_dt
self.river.river_time_sec
self.river.river_nodes
self.river.river_lat
self.river.river_lon
self.river.river_fluxes
self.river.total_flux
self.river.river_temp
self.river.river_salt
self.river.raw_fluxes
self.river.raw_temp
self.river.raw_salt
"""
......@@ -1549,7 +1576,12 @@ class FileReader(Domain):
river_nc = Dataset(river_nc_file, 'r')
time_raw = river_nc.variables['Times'][:]
self.river.time_dt = [datetime.strptime(b''.join(this_time).decode('utf-8').rstrip(), '%Y/%m/%d %H:%M:%S') for this_time in time_raw]
try:
self.river.time_dt = [datetime.strptime(b''.join(this_time).decode('utf-8').rstrip().replace('/', '').replace('-', '').replace('T', '').replace(' ', '').replace(':', '').replace('.', ''), '%Y%m%d%H%M%S%f') for this_time in time_raw]
except:
self.river.time_dt = [datetime.strptime(b''.join(this_time).decode('utf-8').rstrip().replace('/', '').replace('-', '').replace('T', '').replace(' ', '').replace(':', '').replace('.', ''), '%Y%m%d%H%M%S') for this_time in time_raw]
ref_date = datetime(1900, 1, 1)
mod_time_sec = [(this_dt - ref_date).total_seconds() for this_dt in self.time.datetime]
......@@ -1565,12 +1597,18 @@ class FileReader(Domain):
river_flux_raw = river_nc.variables['river_flux'][:, rivers_in_grid]
self.river.river_fluxes = np.asarray([np.interp(mod_time_sec, self.river.river_time_sec, this_col) for this_col in river_flux_raw.T]).T
self.river.total_flux = np.sum(self.river.river_fluxes, axis=1)
self.river.raw_fluxes = river_flux_raw
river_temp_raw = river_nc.variables['river_temp'][:, rivers_in_grid]
self.river.river_temp = np.asarray([np.interp(mod_time_sec, self.river.river_time_sec, this_col) for this_col in river_temp_raw.T]).T
self.river.raw_temp = river_temp_raw
river_salt_raw = river_nc.variables['river_salt'][:, rivers_in_grid]
self.river.river_salt = np.asarray([np.interp(mod_time_sec, self.river.river_time_sec, this_col) for this_col in river_salt_raw.T]).T
self.river.raw_sal = river_salt_raw
self.river.river_lat = self.grid.lat[self.river.river_nodes]
self.river.river_lon = self.grid.lon[self.river.river_nodes]
river_nc.close()
......
......@@ -21,7 +21,7 @@ import matplotlib.pyplot as plt
import numpy as np
from pandas import read_hdf
from PyFVCOM.grid import get_boundary_polygons, vincenty_distance
from PyFVCOM.grid import get_boundary_polygons, vincenty_distance, node_to_centre
from PyFVCOM.plot import Time, Plotter
from PyFVCOM.read import FileReader
from PyFVCOM.stats import calculate_coefficient, rmse
......@@ -489,9 +489,6 @@ class ValidationComparison():
if self.horizontal_match == 'nearest':
self.chosen_mod_nodes = np.squeeze(np.asarray([self.fvcom_data.closest_element(this_ll) for this_ll in self.chosen_obs_ll]))[:, np.newaxis]
self.chosen_mod_nodes_weights = np.ones(len(self.chosen_mod_nodes))[:, np.newaxis]
# self.chosen_mod_nodes = self.fvcom_data.closest_element(self.chosen_obs_ll)
# self.chosen_mod_nodes_weights = np.ones(len(self.chosen_mod_nodes))
elif self.horizontal_match == 'interp':
chosen_mod_nodes = []
chosen_mod_nodes_weights = []
......@@ -541,12 +538,11 @@ class ValidationComparison():
self.mod_depths = -self.mod_h[:,np.newaxis, :]*self.fvcom_data.grid.siglay[np.newaxis, :,:]
elif self.mode == 'elements':
setattr(self.fvcom_data.data, 'zeta_centre', pf.grid.node_to_centre(self.fvcom_data.zeta, self.fvcom_data))
self.mod_h = self.fvcom_data.data.zeta_centre + self.grid.h_center
setattr(self.fvcom_data.data, 'zeta_centre', node_to_centre(self.fvcom_data.data.zeta, self.fvcom_data))
self.mod_h = self.fvcom_data.data.zeta_centre + self.fvcom_data.grid.h_center
self.mod_depths = self.mod_h[:,np.newaxis, :]*self.fvcom_data.grid.siglay_center[np.newaxis, :,:]
self.mod_obs_depths_all_t = np.sum(self.mod_depths[:,:,self.chosen_mod_nodes]*self.chosen_mod_nodes_weights[np.newaxis, np.newaxis,:], axis=-1)
self.mod_obs_depths = np.diagonal(np.squeeze(self.mod_obs_depths_all_t[self.chosen_mod_times,:]))
self.mod_obs_depths = np.asarray([np.sum(np.sum(self.mod_depths[:,:,self.chosen_mod_nodes[i,...]]*self.chosen_mod_nodes_weights[i,...], axis=-1)[self.chosen_mod_times[i,...],:]*self.chosen_mod_times_weights[i,...],axis=0) for i in np.arange(0,len(self.chosen_mod_times))])
if self.vertical_match == 'nearest':
self.chosen_mod_depths = np.asarray([np.argmin(np.abs(this_mod_obs_dep - this_dep)) for this_mod_obs_dep, this_dep in zip(self.mod_obs_depths, self.obs_data.grid.depth[self.chosen_obs])])[:,np.newaxis]
......
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