ufs-coastal DATM+SCHISM DuckNC #127
Comments
|
@mansurjisan next steps discussed today, per my notes (plz edit/comment as needed)
|
|
I completed the following simulations and uploaded the water elevation and wave height plots to Google Drive.
|
|
@mansurjisan I used your scripts to generate input files for the RT. That are working fine but I did not integrate them to the workflow or try to run the case yet. I have also couple of suggestion and question about the scripts,
Is there any way to get those files from public space. Do we need to register?
Anyway, I think overall it is great to have these script. Thanks for your help. I am plaining to move these scripts to app level and integrate with the workflow. |
@mansurjisan I saw you tested replacing seawater with gsw in your pyschism fork right? |
|
@mansurjisan I have one minor question. In |
Hi @uturuncoglu , Regarding your questions about gen_bnd.py:
For the Duck, NC RT case, I am using a different approach to generate required elev.th file. Pyschism currently doesn't have the function to generate this file. I've used OSU OTPS for this, but I'm working on generating it directly with pyschism. |
@janahaddad, yes, you're correct! I've addressed this issue in my fork of pyschism. Here's what I did:
If you're interested in implementing these changes, I'd be happy to provide more details or assist with any questions you might have about the modification process. |
|
@mansurjisan @janahaddad JFYI, I integrated https://github.com/oceanmodeling/ufs-coastal-app/tree/feature/workflow I think that |
|
@mansurjisan BTW, why you are calling |
|
@mansurjisan |
Hi @uturuncoglu , thanks for letting me know. I'll look into this as soon as the Hercules cluster is back online. |
For most of our RT cases, we will likely use either TPXO or HYCOM to generate boundary conditions. Regarding TPXO, since it requires user registration to download the files, we may not be able to distribute those files via GitHub. One potential solution is to use HYCOM for all RT cases, which would avoid this issue. We can discuss this further and explore what steps are needed during our Monday meeting. |
Hi @uturuncoglu , I’ve uploaded the missing hgrid.ll file to the /data directory, which should resolve the issue with gen_bctides.py. You can find it here: GitHub link. Please give it a try and let me know if you encounter any further issues. |
|
@mansurjisan Thanks but I think that the issue with |
Hi @uturuncoglu , I see, thanks for the clarification. You can access the TPXO dataset files directly from my directory on Hercules: Unfortunately, due to TPXO registration restrictions, I’m unable to upload the files to GitHub. Let me know if this resolves the issue or if you encounter any other problems. |
|
@mansurjisan Yes. I am using same location but when I run the script. It could not find the required files and I think it is not looking your directory. So, I wonder if you tested that script recently. |
Hi @uturuncoglu , Yes, I just ran the script and it was able to access the required files without any issues. I also placed these files in |
|
@mansurjisan Yes, I think that it is the trick. So, it is getting them from |
Hi @uturuncoglu , I think you're right, it seems like the script is defaulting to I'll work on this and keep you updated! |
|
@mansurjisan Thank you. |
|
@uturuncoglu , here is an issue I am currently investigating related to the use of elev.th file in SCHISM within the UFS-Coastal. We can discuss it more during today's meeting. |
Updated simulations with elev2D.th.nc file instead of using elev.th file due to possible bug
SCHISM Required Files Table (Duck, NC RT 1994 Case)
Note: This table includes the core files typically required for each configuration. Additional files might be necessary depending on specific model setups and options chosen in the parameter files.
Next steps
|
|
related to Implementing Flexible Build Steps in UFS-Coastal Using ecFlow's Event-based Triggers |
Updated table with DATM+SCHISM run:
|
|
@uturuncoglu , I’ve attached a PDF outlining the steps I followed to create the ESMF mesh file for wind forcing using the ERA5 wind data. Comparing it with the actual data, it seems reasonable, but I’d appreciate some clarification on which variable (wind components/mask/topography) we should use when setting up the ESMF mesh file. Let’s go over it in today’s meeting if possible! |
|
Here are some updates on the Duck, NC test case using observed wind speed and water level data from Ali Abdolali's Hurricane Sandy Wave Watch 3 setup.
Water Level
Wind Speed
Wind Direction
|
Interpolating Observed Wind Data into ERA5 GridScript path: Wind_Interp Prerequisitesimport xarray as xr
import pandas as pd
import numpy as np
from metpy.units import units # I used metpy functions to calculate wind components from the provided wind speed and directions. I am not sure whether we can use wind speed and direction directly as input parameters in CDEPS.
from metpy.calc import wind_componentsRequired input files:
Step 1: Reading and Validating Wind ObservationsThe script starts by reading and validating the observed wind data: def read_wind_data(filename):
"""
Read wind data from observed wind speed and direction text file.
"""
# space-separated file with columns: date, time, speed, direction
df = pd.read_csv(filename, sep='\s+', names=['date', 'time', 'speed', 'direction'])
# Convert to numeric, replacing invalid values with NaN
df['speed'] = pd.to_numeric(df['speed'], errors='coerce')
df['direction'] = pd.to_numeric(df['direction'], errors='coerce')
# Remove invalid entries
df = df.dropna()
# Sanity check
df = df[(df['speed'] >= 0) & (df['speed'] < 100) & # reasonable wind speed range
(df['direction'] >= 0) & (df['direction'] <= 360)] # valid direction range
# datetime index
df['datetime'] = pd.to_datetime(df['date'] + ' ' + df['time'])
df = df.set_index('datetime')
return dfStep 2: Converting Wind ComponentsConvert observed wind speed and direction to U and V components: def calculate_wind_components(speed, direction):
"""
Calculate U and V components using MetPy with error handling.
"""
speed = np.clip(speed, 0, 100) * units('m/s') # clip to reasonable range
direction = np.clip(direction, 0, 360) * units('degrees')
u, v = wind_components(speed, direction)
return float(u.magnitude), float(v.magnitude)Step 3: Processing Wind ObservationsCreate time series of U and V components from observations: def process_wind_observations(wind_df):
"""
Process wind observations to create time series of U and V components
"""
u_series = []
v_series = []
for idx, row in wind_df.iterrows():
u, v = calculate_wind_components(row['speed'], row['direction'])
u_series.append(u)
v_series.append(v)
wind_df['u10'] = u_series
wind_df['v10'] = v_series
return wind_dfStep 4: Main Interpolation ProcessThe main interpolation function handles multiple aspects of data processing: def interpolate_era5_with_obs_wind(ds, wind_df, n_timesteps=None):
"""
Interpolate observed wind data in to ERA5 grid with 30-minute intervals.
"""
# 4.1: Initial Setup
if n_timesteps is None:
n_timesteps = len(ds.valid_time)
# Limit to specified number of timesteps
ds = ds.isel(valid_time=slice(0, n_timesteps))
# 4.2: Time Handling
# Create new time array with 30-min intervals
time_orig = pd.to_datetime(ds.valid_time.values, unit='s')
time_new = pd.date_range(start=time_orig[0],
end=time_orig[-1],
freq='30min')
# Convert to unix timestamp for xarray
time_new_unix = (time_new - pd.Timestamp("1970-01-01")) // pd.Timedelta("1s")
# 4.3: Wind Observation Processing
# Process wind observations to get U/V components
wind_df = process_wind_observations(wind_df)
# Ensure wind observations cover the required time period
wind_df = wind_df.reindex(pd.date_range(start=time_orig[0],
end=time_orig[-1],
freq='30min'))
# Interpolate any gaps in wind data
wind_df = wind_df.interpolate(method='linear')
# 4.4: Initialize New Dataset
ds_new = xr.Dataset(
coords={
'valid_time': time_new_unix,
'latitude': ds.latitude,
'longitude': ds.longitude
}
)
# 4.5: MSL Pressure Interpolation
print("Interpolating msl...")
msl_data = ds['msl'].values
new_msl = np.zeros((len(time_new),) + msl_data.shape[1:])
# Manual interpolation for MSL
for i in range(len(msl_data)-1):
idx = i * 2
# Copy original hour value
new_msl[idx] = msl_data[i]
# Calculate interpolated value for 30-min mark
if idx + 1 < len(time_new):
new_msl[idx + 1] = (msl_data[i] + msl_data[i + 1]) / 2
# Set final timestep
if len(time_new) > 0:
new_msl[-1] = msl_data[-1]
# Add MSL to new dataset
ds_new['msl'] = (('valid_time', 'latitude', 'longitude'), new_msl)
ds_new['msl'].attrs = ds['msl'].attrs
# 4.6: Wind Component Processing
# Initialize wind component arrays
u_data = np.zeros((len(time_new),) + ds['u10'].shape[1:])
v_data = np.zeros((len(time_new),) + ds['v10'].shape[1:])
# Fill wind component arrays with interpolated observations
for t in range(len(time_new)):
# Assign uniform wind field from observations
u_data[t,:,:] = wind_df['u10'].iloc[t]
v_data[t,:,:] = wind_df['v10'].iloc[t]
if t % 10 == 0:
print(f"Timestep {t+1}/{len(time_new)}")
print(f"U10={wind_df['u10'].iloc[t]:.2f} m/s, V10={wind_df['v10'].iloc[t]:.2f} m/s")
# 4.7: Add Wind Components to Dataset
ds_new['u10'] = (('valid_time', 'latitude', 'longitude'), u_data)
ds_new['v10'] = (('valid_time', 'latitude', 'longitude'), v_data)
# Copy variable attributes
ds_new['u10'].attrs = ds['u10'].attrs
ds_new['v10'].attrs = ds['v10'].attrs
# 4.8: Copy Remaining Variables and Attributes
# Copy any remaining variables
for var in ds.variables:
if var not in ['u10', 'v10', 'msl', 'valid_time', 'latitude', 'longitude']:
ds_new[var] = ds[var]
# Copy coordinate attributes
for coord in ['latitude', 'longitude']:
ds_new[coord].attrs = ds[coord].attrs
# Set time attributes
ds_new.valid_time.attrs = {
'long_name': 'time',
'standard_name': 'time',
'units': 'seconds since 1970-01-01',
'calendar': 'proleptic_gregorian'
}
# Copy global attributes
ds_new.attrs = ds.attrs
return ds_newStep 5: Final Processing and Output
# Rename valid_time to time
ds_30min = ds_30min.rename({'valid_time': 'time'})
# Invert latitudes
ds_30min = ds_30min.reindex(latitude=ds_30min.latitude[::-1])
# Set encoding for output
encoding = {
'time': {'dtype': 'int64', '_FillValue': None},
'u10': {'dtype': 'float32', '_FillValue': -9999.0},
'v10': {'dtype': 'float32', '_FillValue': -9999.0},
'msl': {'dtype': 'float32', '_FillValue': -9999.0}
}
# Save to netCDF
ds_30min.to_netcdf('era5_data_30min_obs_wind_rot_fix.nc', encoding=encoding) |
Input files used in the simulation are: Configuration Files
Grid Files
Boundary and Initial Conditions
Data Atmosphere Files
Runtime Files
Output Control Files
System Directories
|
|
@mansurjisan Re: pressure data: some scripts to pull & concatenate these here if you want to take a look at it |
Thank you, Jana! I will look in to this and try to interpolate in to the atmospheric forcing file. |
(9/10/2024 moving @mansurjisan updates from SSM team repo)
Note that this Issue also includes steps to be integrated at the app-level workflow oceanmodeling/ufs-coastal-app#8 & oceanmodeling/ufs-coastal-app#9, and testing for other NOAA-OCS workflows outside ufs-coastal (noaa-ocs-modeling/SurgeTeamCoordination/issues/563)
Repository Link
UFS-Coastal-Inputs on GitHub
Objective
The primary goal of this project is to develop scripts that will be integrated into the UFS coastal app workflow to generate input files for regression tests (RT). These scripts ensure consistency and allow for easier validation of the input file generation process across different test cases.
Current Status
Key Features
PySCHISMlibrary for SCHISM model input file generationRecent Activities
gen_gr3_input.pygen_bctides.pygen_bnd.pyScript Details
1. gen_gr3_input.py
/scripts/coastal_ike_shinnecock_atm2sch/gen_gr3_input.pypython gen_gr3_input.py2. gen_bctides.py
/scripts/coastal_ike_shinnecock_atm2sch/gen_bctides.py3. gen_bnd.py
/scripts/coastal_ike_shinnecock_atm2sch/gen_bnd.pypython gen_bnd.pyProject Structure
Next Steps
Recent Activities
Visualizations
Model Grid File
Depth-Averaged Velocity Map
https://drive.google.com/file/d/1r8XpCgZp7EGT_WnVgao2I9VC8m58_rha/view?usp=sharing
Wind Input Visualization
https://drive.google.com/file/d/1dp5VYbAoaY8BG97WDd_qkVMUcfCPIX-_/view?usp=sharing
https://github.com/user-attachments/assets/88284baa-7466-46cf-b70e-3315ad23e5b7
Snapshot of elevation
https://drive.google.com/drive/folders/1aqsMBu_0c4fNmoTMmcEp7HUWygky0vJZ?usp=sharing
The text was updated successfully, but these errors were encountered: