Reading NOAA STOFS (from ADCIRC model) data

[marceloandrioni]

Hello, I was wondering if it is possible to use xugrid to read and subset NOAA STOFS data created by the ADCIRC model.

The netcdf files are available in AWS buckets like this one

Bellow an example of one the smaller files:

netcdf stofs_2d_glo_fcst.61 {
dimensions:
	time = UNLIMITED ; // (960 currently)
	station = 1687 ;
	namelen = 50 ;
variables:
	double time(time) ;
		time:long_name = "model time" ;
		time:standard_name = "time" ;
		time:units = "seconds since 2022-11-24 00:00:00        ! NCDASE - BASE_DAT" ;
		time:base_date = "2022-11-24 00:00:00        ! NCDASE - BASE_DATE" ;
	char station_name(station, namelen) ;
		station_name:long_name = "station name" ;
	double x(station) ;
		x:long_name = "longitude" ;
		x:standard_name = "longitude" ;
		x:units = "degrees_east" ;
		x:positive = "east" ;
	double y(station) ;
		y:long_name = "latitude" ;
		y:standard_name = "latitude" ;
		y:units = "degrees_north" ;
		y:positive = "north" ;
	double zeta(time, station) ;
		zeta:long_name = "water surface elevation above geoid" ;
		zeta:standard_name = "sea_surface_height_above_geoid" ;
		zeta:units = "m" ;
		zeta:_FillValue = -99999. ;

// global attributes:
		:_FillValue = -99999. ;
		:dry_Value = -99999. ;
		:model = "ADCIRC" ;
		:version = "0d82f1a" ;
		:grid_type = "Triangular" ;
		:description = "2023062300 :-6 hr nowcast and +180 hr forecast ! 32 CHARACTER ALPHANUMERIC RUN D" ;
		:agrid = "OceanMesh2D" ;
		:rundes = "2023062300 :-6 hr nowcast and +180 hr forecast ! 32 CHARACTER ALPHANUMERIC RUN D" ;
		:runid = "STOFS 2D GLOBAL v5.3       ! 24 CHARACTER ALPHANUMERIC RUN IDENTIFICATION" ;
		:title = "STOFS_2D_GLOBAL.V1.1.0     ! NCPROJ - PROJECT TITLE" ;
		:institution = "NOS/OCS/CSDL/CMMB          ! NCINST - PROJECT INSTITUTION" ;
		:source = "Dogwood/Cactus             ! NCSOUR - PROJECT SOURCE" ;
		:history = "PRODUCTION                 ! NCHIST - PROJECT HISTORY" ;
		:references = "http://www.adcirc.org      ! NCREF  - PROJECT REFERENCES" ;
		:comments = "STOFS_2D_GLOBAL.V1.1.0     ! NCCOM  - PROJECT COMMENTS" ;
		:host = "NOS/OCS/CSDL/CMMB          ! NCHOST - PROJECT HOST" ;
		:convention = "CF-1.0                     ! NCCONV - CONVENTIONS" ;
		:Conventions = "UGRID-0.9.0" ;
		:contact = "[email protected]     ! NCCONT - CONTACT INFORMATION" ;
		:creation_date = "2023-06-23  3:54:51  00:00" ;
		:modification_date = "2023-06-23  3:54:51  00:00" ;
		:fort.15 = "==== Input File Parameters (below) ====" ;
		:dt = 6. ;
		:ihot = 568 ;
		:ics = 22 ;
		:nolibf = 1 ;
		:nolifa = 2 ;
		:nolica = 1 ;
		:nolicat = 1 ;
		:nwp = 7 ;
		:ncor = 1 ;
		:ntip = 2 ;
		:nws = 10 ;
		:nramp = 1 ;
		:tau0 = 0.053333 ;
		:statim = 0. ;
		:reftim = 0. ;
		:rnday = 211. ;
		:dramp = 6.75 ;
		:a00 = 0.8 ;
		:b00 = 0.2 ;
		:c00 = 0. ;
		:h0 = 0.1 ;
		:slam0 = 0. ;
		:sfea0 = 45. ;
		:cf = 0.0005 ;
		:eslm = -0.2 ;
		:cori = 0. ;
		:ntif = 8 ;
		:nbfr = 0 ;
}

The global files are very large (15GB), so I would like to do a subset and save only a smaller region. There is no mesh_topoplogy so I was wondering if it is possible to recreate from the "raw" latitude-longitude coordinates and perform the subset.

Thank you.

[Huite]

Hi @marceloandrioni,

Excuse the tardy reply, I was on vacation for two weeks. I should be able to download some sample data, and see what I encounter. In general, if you want to subset data, you will need the topology. It doesn't have to be stored in the dataset per se, since we can initialize the UGRID topology from a different dataset:

import xugrid as xu
import xarray as xu

topology_ds = xr.open_dataset("topology.nc")
ds = xr.open_dataset("data.nc")

grid = xu.Ugrid2d.from_dataset(topology_ds)
uds = xu.UgridDataset(ds, grid)

This should work, provided the dimension names match nicely between the two dataset. I'll try and see if I can get an example working in the coming days with some actual data!

[marceloandrioni]

Hello @Huite , thank you for the reply. I sent a email inquiring about the grid topology variable to the NOAA researcher responsible for the STOFS model. In the meantime, I was doing some tests with a much smaller ADCIRC file. The file can be read in ncWMS and QGIS (as show here), so I think it follows the UGRID specification correctly.

The file is available here: ADCSWAN_BG_20210111.zip

Do you know if it is possible to apply a latitude-longitude subet using xugrid? Something like this:

import xarray as xr
import xugrid as xu


ds = xr.open_dataset('ADCSWAN_BG_20210111.nc')
grid = xu.Ugrid2d.from_dataset(ds)
uds = xu.UgridDataset(ds, grid)

uds2 = uds.sel(latitude=slice(-23, -22.8),
               longitude=slice(-43, -42.8))

Thank you.

[marceloandrioni]

Hello again,

I just noticed that I should have used:

uds2 = uds.ugrid.sel(y=slice(-23, -22.8),
                     x=slice(-43, -42.8))

uds2.ugrid.to_netcdf('ADCSWAN_BG_20210111_subset.nc')

with this it works and the subset is applied.

Thank you.

[marceloandrioni]

Hello again @Huite.

I just wanted to report that after downloading one of the large global files with actual STOFS output from NOMADS, I noticed that the file has the grid topology variable.

With this I was able to load and subset the global file using xugrid, e.g.:

import xarray as xr
import xugrid as xu


ds = xr.open_dataset('stofs_2d_glo.t00z.fields.cwl.nc',
                     drop_variables=['nvel'])

ds['time'].encoding['units'] = (ds['time'].encoding['units']
                                .split('!')[0].strip())

# get only first time step
ds = ds.isel(time=[0])

uds = xu.UgridDataset(ds)

uds2 = uds.ugrid.sel(y=slice(-28, -18),
                     x=slice(-49, -35))

uds2.ugrid.to_netcdf('stofs_subset.nc')

The subset file displayed just fine using ncWMS2

Now I can use xugrid operationally to only save the subset for my region of interest. I will mark this issued as complete.

Thank you very much.