janahaddad commented 2 months ago

(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

Scripts for the Coastal Ike Shinnecock RT case have been successfully developed.
Successfully executed scripts using an ecFlow-based workflow [noaa-ocs-modeling/SurgeTeamCoordination/issues/563] on the Hercules Cluster.
Next step would be to generate input files for the new Duck, NC RT case.

Key Features

Scripts designed for integration into the UFS coastal app workflow
Uses the PySCHISM library for SCHISM model input file generation
Produces input files identical to those used in existing regression tests
Expandable to new test cases (e.g., Duck, NC)

Recent Activities

Completed development of scripts for the Coastal Ike Shinnecock RT case:
- gen_gr3_input.py
- gen_bctides.py
- gen_bnd.py
Verified that the generated input files are identical to those currently used in the Coastal Ike Shinnecock RT
Prepared for expansion to the Duck, NC test case

Script Details

1. gen_gr3_input.py

Location: /scripts/coastal_ike_shinnecock_atm2sch/gen_gr3_input.py
Purpose: Generates .gr3 format files (albedo, diffmax, diffmin, watertype, windrot_geo2proj, Manning's roughness) to match existing regression test inputs
Current state: Sets constant values for all parameters
Input: Requires hgrid.gr3
Usage: python gen_gr3_input.py

2. gen_bctides.py

Location: /scripts/coastal_ike_shinnecock_atm2sch/gen_bctides.py
Purpose: Generates bctides.in using TPXO, matching the file used in the current regression test

Usage:

python gen_bctides.py ../../data/hgrid.ll 2008-08-23 20.0 '[[3,0,0,0]]' Q1,O1,P1,K1,N2,M2,S2,K2,Mm,Mf,M4,MN4,MS4,2N2,S1 tpxo --earth_tidal_potential Y --cutoff_depth 40

3. gen_bnd.py

Location: /scripts/coastal_ike_shinnecock_atm2sch/gen_bnd.py
Purpose: Generates time series of water surface elevations (elev2D.th.nc file) based on Hycom
Note: This script is developed but not currently used in the existing regression test
Usage: python gen_bnd.py

Project Structure

UFS-Coastal-Inputs/
│
├── data/
│   ├── hgrid.gr3
│   └── vgrid.in
│
├── outputs/
│   ├── albedo.gr3
│   ├── bctides.in
│   ├── diffmax.gr3
│   ├── diffmin.gr3
│   ├── manning.gr3
│   ├── watertype.gr3
│   └── windrot_geo2proj.gr3
│
├── scripts/
│   ├── coastal_ike_shinnecock_atm2sch/
│   │   ├── gen_bctides.py
│   │   ├── gen_bnd.py
│   │   ├── gen_gr3_input.py
│   │   └── README
│   │
│   └── Test_Duck_NC/
│       └── README
│
├── LICENSE
└── README.md

Next Steps

Generate atmospheric forcing files for the RT tests.
Generate input files for the Duck, NC regression test case.

Recent Activities

Made progress on the Duck, NC test case:
- Performed simulation using the standalone SCHISM model setup for Duck Island, NC test case.
- Plotted some output variables; elevation and depth-averaged velocity maps from the model simulation.
- Created an atmospheric forcing file using the ERA5 dataset for the Duck, NC test case.
- Conducted a successful model run with the new forcing file.
- Encountered a challenge when changing the grid from Cartesian to latitude/longitude coordinate system.

Visualizations

Model Grid File

schism_grid_plot_zoomedin

janahaddad commented 2 months ago

@mansurjisan next steps discussed today, per my notes (plz edit/comment as needed)

compare your results with Dan's 1994 obs
Standalone SCHISM with atm for Dorian case
with UFS Coastal for Dorian

mansurjisan commented 2 months ago

I completed the following simulations and uploaded the water elevation and wave height plots to Google Drive.

Model Configuration	Water Elevation	Wave Height	Directory in Hercules
Exp 1: SCHISM+WWM (reference)	Link	Link	`/work/noaa/nosofs/mjisan/schism/schism_verification_tests/Test_WWM_Duck`
Exp 2: SCHISM+WWM+ATM	-	-	`/work/noaa/nosofs/mjisan/schism/schism_verification_tests/Test_WWM_Duck_ATM`
Exp 3: Standalone SCHISM	Link	-	`/work/noaa/nosofs/mjisan/schism/schism_verification_tests/Test_DUCK_SCH`
Exp 4: Standalone SCHISM in UFS	Link	-	`/work/noaa/nosofs/mjisan/ufs-weather-model/tests/stmp/mjisan/FV3_RT/RT_DUCK_NC_SCHISM_STD`

uturuncoglu commented 2 months ago

@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,

The scripts are looking to the ../../data/hgrid.gr3 directory. It would be nice to follow same approach with gen_bctides.py and make inputs as command line argument. This will make easy to integrate them in the app level and also enable flexibility to use with other mesh files.
BTW, gen_bctides.py script is looking for hgrid.ll not hgrid.gr3. Same for vgrid file. So it would be nice to be consistent with the steps and correct README file. I am not sure SCHSIM expects hgrid.ll or hgrid.gr3.

I had also issue with same script (also not sure if LOCAL_TPXO_DATABASE is used or not) as following,

FileNotFoundError: No TPXO file found at "/home/tufuk/.local/share/tpxo/h_tpxo9.v1.nc".
New users will need to register and request a copy of the TPXO9 NetCDF file (specifically `h_tpxo9.v1.nc`) from the authors at https://www.tpxo.net.
Once you obtain `h_tpxo9.v1.nc`, you can follow one of the following options: 
1) copy or symlink the file to "/home/tufuk/.local/share/tpxo/h_tpxo9.v1.nc"
2) set the environment variable `h_tpxo9.v1.nc` to point to the file

Is there any way to get those files from public space. Do we need to register?

gen_bnd.py also looking to data directory. It would be nice to make it command line argument for flexibility.
I am also getting pyschism/forcing/hycom/gofs.py:8: UserWarning: The seawater library is deprecated! Please use gsw instead. import seawater as sw warning. It would be nice to use that package. I know that this is psychic side and maybe we need to contact with them to make that changes.

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.

janahaddad commented 2 months ago

I am also getting pyschism/forcing/hycom/gofs.py:8: UserWarning: The seawater library is deprecated! Please use gsw instead. import seawater as sw warning. It would be nice to use that package.

@mansurjisan I saw you tested replacing seawater with gsw in your pyschism fork right?

uturuncoglu commented 2 months ago

@mansurjisan I have one minor question. In gen_bnd.py script, it seems that ocean_bnd_ids=[0] is fixed to 0. Does it mean we are only creating boundary file for single boundary? To be more generic, I think this part needs to be configurable. Right? I am not sure how this works but maybe it is ordered like North, South, etc. If you could give more information that would be great. BTW, did you test this with multiple boundary?

mansurjisan commented 2 months ago

@mansurjisan I have one minor question. In gen_bnd.py script, it seems that ocean_bnd_ids=[0] is fixed to 0. Does it mean we are only creating boundary file for single boundary? To be more generic, I think this part needs to be configurable. Right? I am not sure how this works but maybe it is ordered like North, South, etc. If you could give more information that would be great. BTW, did you test this with multiple boundary?

Hi @uturuncoglu ,

Regarding your questions about gen_bnd.py:

This script isn't currently used in the Coastal Ike Shinnecock RT test case. It's an unmodified example script from the PySchism tutorial.
You're correct that ocean_bnd_ids=[0] is currently fixed. The number and order of boundaries are actually defined in the .gr3 file, not as cardinal directions.
I used the gen_bctides.py for generating the bctides.in file for Ike Shinnecock test case, but we'll likely need gen_bnd.py for Hycom data in the future.
I plan to modify gen_bnd.py to read the number of boundaries directly from the .gr3 file, making it more flexible for multiple boundaries.
I haven't tested it with multiple boundaries yet, but that's on the to-do list.

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.

mansurjisan commented 2 months ago

I am also getting pyschism/forcing/hycom/gofs.py:8: UserWarning: The seawater library is deprecated! Please use gsw instead. import seawater as sw warning. It would be nice to use that package.

@mansurjisan I saw you tested replacing seawater with gsw in your pyschism fork right?

@janahaddad, yes, you're correct! I've addressed this issue in my fork of pyschism. Here's what I did:

First, I installed the gsw package to replace the deprecated seawater library.
Then, I modified the relevant files in the pyschism codebase, changing instances of import seawater as sw to import gsw
You can see all these changes in my fork of pyschism, which is available here: https://github.com/mansurjisan/pyschism/

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.

uturuncoglu commented 2 months ago

@mansurjisan @janahaddad JFYI, I integrated gen_bnd.py and gen_gr3.py with the workflow (see updated scripts in the ush/schism folder) and I'll continue with the gen_bctides.py. The current implementation is in following link,

https://github.com/oceanmodeling/ufs-coastal-app/tree/feature/workflow

I think that gen_bnd.py and gen_gr3.py are still not generic to cover all the cases. For example, we have values = [2.000000e-1, 1.0, 1e-6, 4, 0.00000000, 2.5000000e-02] gen_gr3.py. I could get those from the config.yaml but I am not sure they could only get single solar value or could be varying spatially. Any idea? Anyway, we could improve them once we could reproduce the DATM+SCHSIM RT with the workflow. Once I have everything, I'll also add documentation. So, you could also test in your end to be sure it is working fine.

uturuncoglu commented 2 months ago

@mansurjisan BTW, why you are calling Bctides twice in the following script? Is there any specific reason? https://github.com/mansurjisan/UFS-Coastal-Inputs/blob/3c918cf78aa7204aafe6cf2dbc741cfcf4cd15c7/scripts/coastal_ike_shinnecock_atm2sch/gen_bctides.py#L120

uturuncoglu commented 2 months ago

@mansurjisan gen_bctides.py script could not find the directory that has input files in my case. There might be an issue with that script. I wonder if you try the version that is available in the repository?

mansurjisan commented 2 months ago

@mansurjisan gen_bctides.py script could not find the directory that has input files in my case. There might be an issue with that script. I wonder if you try the version that is available in the repository?

Hi @uturuncoglu , thanks for letting me know. I'll look into this as soon as the Hercules cluster is back online.

mansurjisan commented 2 months ago

@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,

The scripts are looking to the ../../data/hgrid.gr3 directory. It would be nice to follow same approach with gen_bctides.py and make inputs as command line argument. This will make easy to integrate them in the app level and also enable flexibility to use with other mesh files.

BTW, gen_bctides.py script is looking for hgrid.ll not hgrid.gr3. Same for vgrid file. So it would be nice to be consistent with the steps and correct README file. I am not sure SCHSIM expects hgrid.ll or hgrid.gr3.

I had also issue with same script (also not sure if LOCAL_TPXO_DATABASE is used or not) as following,
FileNotFoundError: No TPXO file found at "/home/tufuk/.local/share/tpxo/h_tpxo9.v1.nc".
New users will need to register and request a copy of the TPXO9 NetCDF file (specifically `h_tpxo9.v1.nc`) from the authors at https://www.tpxo.net.
Once you obtain `h_tpxo9.v1.nc`, you can follow one of the following options: 
1) copy or symlink the file to "/home/tufuk/.local/share/tpxo/h_tpxo9.v1.nc"
2) set the environment variable `h_tpxo9.v1.nc` to point to the file
Is there any way to get those files from public space. Do we need to register?

gen_bnd.py also looking to data directory. It would be nice to make it command line argument for flexibility.

I am also getting pyschism/forcing/hycom/gofs.py:8: UserWarning: The seawater library is deprecated! Please use gsw instead. import seawater as sw warning. It would be nice to use that package. I know that this is psychic side and maybe we need to contact with them to make that changes.

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.

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.

mansurjisan commented 2 months ago

@mansurjisan gen_bctides.py script could not find the directory that has input files in my case. There might be an issue with that script. I wonder if you try the version that is available in the repository?

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.

uturuncoglu commented 2 months ago

@mansurjisan Thanks but I think that the issue with gen_bctides.py is not related with the grid file. It could not find the input files related with the tpxo dataset even if I try to set the environment variable.

mansurjisan commented 2 months ago

@mansurjisan Thanks but I think that the issue with gen_bctides.py is not related with the grid file. It could not find the input files related with the tpxo dataset even if I try to set the environment variable.

Hi @uturuncoglu , I see, thanks for the clarification. You can access the TPXO dataset files directly from my directory on Hercules: /work/noaa/nosofs/mjisan/pyschism-main/PySCHISM_tutorial/data/TPXO

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.

uturuncoglu commented 2 months ago

@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.

mansurjisan commented 2 months ago

@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 /home/mjisan/.local/share/tpxo , so it's possible the script might be reading from that shared directory.

uturuncoglu commented 2 months ago

@mansurjisan Yes, I think that it is the trick. So, it is getting them from /home/mjisan/.local not the original location even if you set environment variable in the top of the script. Is it possible to fix that part. So, it could pick from the location set by environment variable. I try different things but I could not have success. It would be nice to have that capability. BTW, we could place the data somewhere on the platform for our usage. Of course, external users could get and place it to their .local directory.

mansurjisan commented 2 months ago

@mansurjisan Yes, I think that it is the trick. So, it is getting them from /home/mjisan/.local not the original location even if you set environment variable in the top of the script. Is it possible to fix that part. So, it could pick from the location set by environment variable. I try different things but I could not have success. It would be nice to have that capability. BTW, we could place the data somewhere on the platform for our usage. Of course, external users could get and place it to their .local directory.

Hi @uturuncoglu , I think you're right, it seems like the script is defaulting to /home/mjisan/.local/share even when the environment variable is set. I'll take a closer look at that part of the script to make sure it picks up the location set by the environment variable instead. Also, I agree, it would be a good idea to place the data somewhere that is easily accessible.

I'll work on this and keep you updated!

uturuncoglu commented 2 months ago

@mansurjisan Thank you.

mansurjisan commented 1 month ago

@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.

elev.th_error_handling.pdf

mansurjisan commented 1 month ago

Updated simulations with elev2D.th.nc file instead of using elev.th file due to possible bug

Model Configuration	Water Elevation	Wave Height	Directory in Hercules
Exp 1: SCHISM+WWM (reference)	Link	Link	`/work/noaa/nosofs/mjisan/schism/schism_verification_tests/Test_WWM_Duck_Elev2D_SCH_WWM`
Exp 2: SCHISM+WWM+ATM	Link	Link	`/work/noaa/nosofs/mjisan/schism/schism_verification_tests/Test_WWM_Duck_Elev2D_SCH_WWM_ATM`
Exp 3: SCHISM+ATM	Link	-	`/work/noaa/nosofs/mjisan/schism/schism_verification_tests/Test_WWM_Duck_Elev2D_SCH_ATM`
Exp 4: Standalone SCHISM	Link	-	`/work/noaa/nosofs/mjisan/schism/schism_verification_tests/Test_WWM_Duck_Elev2D_SCH_ST`
Exp 5: Standalone SCHISM in UFS	Link	-	`/work/noaa/nosofs/mjisan/ufs-weather-model/tests/stmp/mjisan/FV3_RT/RT_DUCK_NC_SCHISM_STD`

SCHISM Required Files Table (Duck, NC RT 1994 Case)

File Name	SCHISM Standalone	SCHISM + WWM	SCHISM + WWM + ATM
hgrid.gr3	✓	✓	✓
vgrid.in	✓	✓	✓
param.nml	✓	✓	✓
bctides.in	✓	✓	✓
manning.gr3	✓	✓	✓
elev.ic	✓	✓	✓
rough.gr3	✓	✓	✓
elev2D.th.nc	✓	✓	✓
wwminput.nml		✓	✓
wwmbnd.gr3		✓	✓
wwmbnd.XY		✓	✓
windrot_geo2proj.gr3		✓	✓
sflux_air_1.0001.nc			✓
sflux_prc_1.0001.nc			✓
sflux_inputs.txt			✓

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.

Wrote a script that converts the type 1 boundary condition to type 4: https://github.com/mansurjisan/SCHISM-DUCK-NC-RT/blob/main/write_elev2dnc.py (required input elev.th)

Next steps

run these cases within UFS-Coastal framework
update the existing ecFlow workflow code to preprocess, execute and post-process output files within UFS-Coastal.

mansurjisan commented 1 month ago

related to Implementing Flexible Build Steps in UFS-Coastal Using ecFlow's Event-based Triggers https://github.com/noaa-ocs-modeling/SurgeTeamCoordination/issues/627

mansurjisan commented 3 weeks ago

Updated table with DATM+SCHISM run:

Model Configuration	Water Elevation	Wave Height	Directory in Hercules
Exp 1: SCHISM+WWM (reference)	Link	Link	`/work/noaa/nosofs/mjisan/schism/schism_verification_tests/Test_WWM_Duck_Elev2D_SCH_WWM`
Exp 2: SCHISM+WWM+ATM	Link	Link	`/work/noaa/nosofs/mjisan/schism/schism_verification_tests/Test_WWM_Duck_Elev2D_SCH_WWM_ATM`
Exp 3: SCHISM+ATM	Link	-	`/work/noaa/nosofs/mjisan/schism/schism_verification_tests/Test_WWM_Duck_Elev2D_SCH_ATM`
Exp 4: Standalone SCHISM	Link	-	`/work/noaa/nosofs/mjisan/schism/schism_verification_tests/Test_WWM_Duck_Elev2D_SCH_ST`
Exp 5: Standalone SCHISM in UFS	Link	-	`/work/noaa/nosofs/mjisan/ufs-weather-model/tests/stmp/mjisan/FV3_RT/RT_DUCK_NC_SCHISM_STD`
Exp 6: DATM+SCHISM in UFS	[]()	-	`/work/noaa/nosofs/mjisan/ufs-weather-model/tests/stmp/mjisan/FV3_RT/RT_DUCK_DATM_SCH/`

mansurjisan commented 3 weeks ago

@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!

ESMF_ATM_Mesh_for_DuckNC.pdf

mansurjisan commented 2 weeks ago

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.

I extracted the water level and wind speed data from ww3_shel.nml file. Attached are the time series plots based on the data.

Water Level

water_level_analysis2

Wind Speed

wind_speed_analysis

Wind Direction

wind_direction_analysis

mansurjisan commented 2 weeks ago

I used the timeseries of water elevation data as boundary condition in the simulation. I used the type 4 boundary condition (elev2D.th.nc) since type 1 boundary condition (elev.th) can't be used due to a known bug. I used the write_elev2dnc.py to convert the data from elev.th to elev2D.th.nc file.
For wind forcing, I interpolated observed wind forcing onto ERA5 mesh at every 30-minute interval (that is the interval at which the observations are provided). So, wind speed remains constant over the entire domain and updates with every time step. Below are the steps I followed for interpolation:

Interpolating Observed Wind Data into ERA5 Grid

Script path: Wind_Interp

Prerequisites

import 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_components

Required input files:

ERA5 netCDF file (e.g., 'era5_data_201210_sandy.nc') . This is the original ERA5 data for Hurricane Sandy over which I interpolated the observed data.
Observed wind data file (space-separated text file with columns: date, time, speed, direction)

Step 1: Reading and Validating Wind Observations

The 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 df

Step 2: Converting Wind Components

Convert 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 Observations

Create 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_df

Step 4: Main Interpolation Process

The 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_new

Step 5: Final Processing and Output

Renaming valid_time to time
Inverting latitudes
Setting proper encoding for netCDF output
Explicit fill value specification in encoding


# 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)

mansurjisan commented 2 weeks ago

I completed the model run for the Hurricane Sandy case for the Duck, NC domain for the CDEPS and SCHISM coupled configuration. Model run is located in this directory in Hercules: /work/noaa/nosofs/mjisan/ufs-weather-model/tests/stmp/mjisan/FV3_RT/RT_DUCK_DATM_SCH_V2

Input files used in the simulation are:

Configuration Files

File	Description
`model_configure`	UFS model configuration settings
`ufs.configure`	UFS system configuration file
`param.nml`	SCHISM namelist file
`datm_in`	Data atmosphere input configuration
`fd_ufs.yaml`	UFS field dictionary configuration

Grid Files

File	Description
`hgrid.gr3`	Horizontal grid file
`hgrid.ll`	Horizontal grid in lat/lon format
`vgrid.in`	Vertical grid configuration
`windrot_geo2proj.gr3`	Wind rotation grid file
`manning.gr3`	Manning's roughness coefficient grid
`diffmax.gr3`	Maximum diffusivity grid
`diffmin.gr3`	Minimum diffusivity grid
`rough.gr3`	Surface roughness grid

Boundary and Initial Conditions

File	Description
`bctides.in`	Tidal boundary conditions
`elev.ic`	Initial water elevation conditions
`elev2D.th.nc`	2D elevation time history

Data Atmosphere Files

File	Description
`datm.streams`	Data atmosphere stream configuration
`era5/`	Directory containing ERA5 atmospheric forcing data

Runtime Files

File	Description
`fv3.exe`	executable
`job_card`	Job submission configuration

Output Control Files

File	Description
`station.in`	Station locations for point output
`rpointer.atm`	Atmosphere restart pointer
`rpointer.cpl`	Coupler restart pointer

System Directories

Directory	Description
`outputs/`	Model output files
`RESTART/`	Restart files directory
`era5/`	Wind input files directory
`modulefiles/`	Module files for environment setup

mansurjisan commented 2 weeks ago

Animations of the model result

a. Water level

b. Surface wind speed

c. Comparisons of water level between without atmospheric forcing and with atmospheric forcing simulation

janahaddad commented 1 week ago

@mansurjisan Re: pressure data: Not sure if this is worth exploring/doing, but there is pressure data at FRF here, I'd guess that Ali's wind observations are from the same lat/lon as these (thre's also FRF wind speed and direction obs here)

some scripts to pull & concatenate these here if you want to take a look at it

mansurjisan commented 1 week ago

@mansurjisan Re: pressure data: Not sure if this is worth exploring/doing, but there is pressure data at FRF here, I'd guess that Ali's wind observations are from the same lat/lon as these (thre's also FRF wind speed and direction obs here)

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.

oceanmodeling / ufs-weather-model

ufs-coastal DATM+SCHISM DuckNC #127

Repository Link

Objective

Current Status

Key Features

Script Details

1. gen_gr3_input.py

2. gen_bctides.py

3. gen_bnd.py

Project Structure

Next Steps

Recent Activities

Visualizations

Model Grid File

Depth-Averaged Velocity Map

Wind Input Visualization

Snapshot of elevation

Updated simulations with elev2D.th.nc file instead of using elev.th file due to possible bug

Next steps

Updated table with DATM+SCHISM run:

Water Level

Wind Speed

Wind Direction

Interpolating Observed Wind Data into ERA5 Grid

Prerequisites

Step 1: Reading and Validating Wind Observations

Step 2: Converting Wind Components

Step 3: Processing Wind Observations

Step 4: Main Interpolation Process

Step 5: Final Processing and Output

Configuration Files

Grid Files

Boundary and Initial Conditions

Data Atmosphere Files

Runtime Files

Output Control Files

System Directories

Animations of the model result