Skip to content

Commit

Permalink
add Raymond's update
Browse files Browse the repository at this point in the history
  • Loading branch information
@qiangshu
qiangshu committed Oct 31, 2024
1 parent bd79a1e commit 482c6cb
Showing 1 changed file with 69 additions and 55 deletions.
124 changes: 69 additions & 55 deletions bdschism/bdschism/port_boundary.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,43 +8,60 @@
Script to convert various data formats (from calsim, csv files) into
SCHISM flux, salt and temp time history (.th) files
"""
import os
import yaml
from argparse import ArgumentParser
import pandas as pd
import matplotlib.pylab as plt
from vtools.data.vtime import minutes
from vtools.functions.filter import ts_gaussian_filter
from vtools.functions.interpolate import rhistinterp
from schimpy.unit_conversions import CFS2CMS, ec_psu_25c
from pyhecdss import get_ts
from schimpy import model_time
from pyhecdss import get_ts
from schimpy.unit_conversions import CFS2CMS, ec_psu_25c
from vtools.functions.interpolate import rhistinterp
from vtools.functions.filter import ts_gaussian_filter
from vtools.data.vtime import minutes
import matplotlib.pylab as plt
import pandas as pd
from argparse import ArgumentParser
import yaml
import os

dir = os.path.dirname(__file__)

source_map_file = '../examples/port_boundary_examples/port_boundary_map.csv'
schism_flux_file = '../data/time_history/flux.th'
schism_salt_file = '../data/time_history/salt.th'
schism_temp_file = '../data/time_history/temp.th'
out_file_flux = 'flux.th.ported'
out_file_salt = 'salt.th.ported'
out_file_temp = 'temp.th.ported'
boundary_kinds = ['flux']
sd = [2014,1,1]
ed = [2014,12,31]

with open(os.path.join(dir, 'config.yaml'), 'r') as f:
config = yaml.safe_load(f)


# Read in parameters from the YAML file
source_map_file = os.path.join(dir, config['file']['source_map_file'])
schism_flux_file = os.path.join(dir, config['file']['schism_flux_file'])
schism_salt_file = os.path.join(dir, config['file']['schism_salt_file'])
schism_temp_file = os.path.join(dir, config['file']['schism_temp_file'])
out_file_flux = os.path.join(config['file']['out_file_flux'])
out_file_salt = os.path.join(config['file']['out_file_salt'])
out_file_temp = os.path.join(config['file']['out_file_temp'])
boundary_kinds = config['param']['boundary_kinds']
sd = config['param']['start_date']
ed = config['param']['end_date']

dt = minutes(15)
start_date = pd.Timestamp(sd[0], sd[1], sd[2])
end_date = pd.Timestamp(ed[0], ed[1], ed[2])
df_rng = pd.date_range(start_date, end_date, freq=dt)
source_map = pd.read_csv(source_map_file, header=0)

def read_csv(file, var, name,p=2.):
# Read in the reference SCHISM flux, salt and temperature files
# to be used as a starting point and to substitute timeseries not
# available from other data sources.

flux = pd.read_csv(schism_flux_file, index_col=0, parse_dates=[0],
sep="\\s+", header=0)
salt = pd.read_csv(schism_salt_file, header=0, parse_dates=True,
index_col=0, sep="\\s+")
temp = pd.read_csv(schism_temp_file, header=0, parse_dates=True,
index_col=0, sep="\\s+")


def read_csv(file, var, name, p=2.):
"""
Reads in a csv file of monthly boundary conditions and interpolates
Outputs an interpolated DataFrame of that variable
"""
forecast_df = pd.read_csv(os.path.join(dir,file), index_col=0, header=0,
forecast_df = pd.read_csv(os.path.join(dir, file), index_col=0, header=0,
parse_dates=True)
forecast_df.index = forecast_df.index.to_period('M')
interp_series = rhistinterp(forecast_df[var].astype('float'),
Expand All @@ -53,49 +70,40 @@ def read_csv(file, var, name,p=2.):
interp_df[[name]] = pd.DataFrame({var: interp_series})
return interp_df

def read_dss(file,pathname,sch_name=None,p=2.):

def read_dss(file, pathname, sch_name=None, p=2.):
"""
Reads in a DSM2 dss file and interpolates
Outputs an interpolated DataFrame of that variable
"""
ts15min = pd.DataFrame()
ts=get_ts(os.path.join(dir,file), pathname = pathname)
print(pathname)
ts = get_ts(os.path.join(dir, file), pathname)
for tsi in ts:
b=(tsi[0].columns.values[0]).split("/")[2]
c=(tsi[0].columns.values[0]).split("/")[3]
f=(tsi[0].columns.values[0]).split("/")[6]
b = (tsi[0].columns.values[0]).split("/")[2]
c = (tsi[0].columns.values[0]).split("/")[3]
f = (tsi[0].columns.values[0]).split("/")[6]
if p != 0:
ts15min[[sch_name]]=rhistinterp(tsi[0],dt,p=p).reindex(df_rng)
ts15min[[sch_name]] = rhistinterp(tsi[0], dt, p=p).reindex(df_rng)
else:
ts15min[[sch_name]]= tsi[0]
ts15min[[sch_name]] = tsi[0]
print("Reading " + b + " " + f)
if ts15min.empty:
raise ValueError(f'Warning: DSS data not found for {b}')
return ts15min


for boundary_kind in boundary_kinds:

source_map = source_map.loc[source_map['boundary_kind'] == boundary_kind]

"""
Read in the reference SCHISM flux, salt and temperature files
to be used as a starting point and to substitute timeseries not
available from other data sources.
"""
if boundary_kind == 'flux':
flux = pd.read_csv(schism_flux_file,index_col=0,parse_dates=[0],
sep="\\s+",header=0)

if boundary_kind == 'flow':
dd = flux.copy().reindex(df_rng)
out_file = out_file_flux
elif boundary_kind == 'ec':
salt = pd.read_csv(schism_salt_file,header=0,parse_dates=True,
index_col=0,sep="\s+")
dd = salt.copy().reindex(df_rng)
out_file = out_file_salt
elif boundary_kind == 'temp':
temp = pd.read_csv(schism_temp_file,header=0,parse_dates=True,
index_col=0,sep="\s+")
dd = temp.copy().reindex(df_rng)
out_file = out_file_temp

Expand All @@ -104,7 +112,7 @@ def read_dss(file,pathname,sch_name=None,p=2.):
name = row['schism_boundary']
source_kind = row['source_kind']
source_file = str(row['source_file'])
derived = str(row['derived']).capitalize()=='True'
derived = str(row['derived']).capitalize() == 'True'
var = row['var']
sign = row['sign']
convert = row['convert']
Expand All @@ -124,13 +132,13 @@ def read_dss(file,pathname,sch_name=None,p=2.):
csv = pd.DataFrame()
vars = var.split(';')
for v in vars:
csv[[v]] = read_csv(source_file, v, name,p = p)
csv[[v]] = read_csv(source_file, v, name, p=p)
dts = eval(formula).to_frame(name).reindex(df_rng)
dfi = ts_gaussian_filter(dts, sigma=100)
else:
print(f"Updating SCHISM {name} with interpolated monthly\
forecast {var}")
dfi = read_csv(source_file, var, name,p = p)
dfi = read_csv(source_file, var, name, p=p)

elif source_kind == 'DSS':
# Substitute in CalSim value.
Expand All @@ -141,13 +149,13 @@ def read_dss(file,pathname,sch_name=None,p=2.):
dss = pd.DataFrame()
for pn in vars:
b = pn.split("/")[2]
dss[[b]] = read_dss(source_file,pathname = pn,
sch_name = name,p = p)
dss[[b]] = read_dss(source_file, pathname=pn,
sch_name=name, p=p)
dts = eval(formula).to_frame(name).reindex(df_rng)
dfi = ts_gaussian_filter(dts, sigma=100)
else:
else:
print(f"Updating SCHISM {name} with DSS variable {var}")
dfi = read_dss(source_file,pathname = var,sch_name = name,p = p)
dfi = read_dss(source_file, pathname=var, sch_name=name, p=p)

elif source_kind == 'CONSTANT':
# Simply fill with a constant specified.
Expand All @@ -156,15 +164,21 @@ def read_dss(file,pathname,sch_name=None,p=2.):

# Do conversions.
if convert == 'CFS_CMS':
dfi = dfi*CFS2CMS*sign
dfi = dfi * CFS2CMS * sign
elif convert == 'EC_PSU':
dfi = ec_psu_25c(dfi)*sign
dfi = ec_psu_25c(dfi) * sign
else:
dfi = dfi


# Trim dfi so that it starts where flux ends, so that dfi doesn't
# overwrite any historical data
dfi = dfi[dfi.index >= flux.index[-1]]

# Update the dataframe.
dd.update(dfi, overwrite=True)

print(dfi)

print(dd)

# Format the outputs.
Expand All @@ -179,7 +193,7 @@ def read_dss(file,pathname,sch_name=None,p=2.):
sep=" ")

dd.plot()


print('Done')
plt.show()

0 comments on commit 482c6cb

Please sign in to comment.