search

NOAA-ORR-ERD / gridded

A single API for accessing / working with gridded model results on multiple grid types
https://noaa-orr-erd.github.io/gridded/index.html
The Unlicense
64 stars 14 forks source link

Loading and plotting DELFT-3D Flow NetCDF output #47

Closed JulesBlm closed 4 years ago

JulesBlm commented 4 years ago

Hi all, I'm trying to plot DELFT-3D FLOW NetCDF output. Plotting top-down plots of variables is fairly straightforward with xarray and matplotlib, but I'm having a hard time plotting vertical sections on an ocean sigma-grid. The output follow CF-1.6 and SGRID-0.3 conventions. I'm loading the NetCDF file using xarray which gives me the following info.

<xarray.Dataset>
Dimensions:        (KMAXOUT: 81, KMAXOUT_RESTR: 80, LSED: 4, LSEDTOT: 4, LSTSCI: 4, LTUR: 2, M: 42, MC: 42, N: 62, NC: 62, NSRC: 21, SIG_INTF: 81, SIG_LYR: 80, avgtime: 1, length_7: 7, nlyr: 77, nlyrp1: 78, time: 49)
Coordinates:
    XCOR           (MC, NC) float32 ...
    YCOR           (MC, NC) float32 ...
    XZ             (M, N) float32 ...
    YZ             (M, N) float32 ...
  * SIG_LYR        (SIG_LYR) float32 -0.05 -0.15 -0.225 ... -0.99925 -0.99975
  * SIG_INTF       (SIG_INTF) float32 0.0 -0.1 -0.2 ... -0.999 -0.9995 -1.0
  * KMAXOUT        (KMAXOUT) int32 0 1 2 3 4 5 6 7 8 ... 73 74 75 76 77 78 79 80
  * KMAXOUT_RESTR  (KMAXOUT_RESTR) int32 0 1 2 3 4 5 6 ... 73 74 75 76 77 78 79
  * time           (time) datetime64[ns] 2019-05-02 ... 2019-05-02T04:00:00
Dimensions without coordinates: LSED, LSEDTOT, LSTSCI, LTUR, M, MC, N, NC, NSRC, avgtime, length_7, nlyr, nlyrp1
Data variables:
    ALFAS          (M, N) float32 ...
    KCU            (MC, N) int32 ...
    KCV            (M, NC) int32 ...
    KCS            (M, N) int32 ...
    DP0            (MC, NC) float32 ...
    DPS0           (M, N) float32 ...
    DPU0           (MC, N) float32 ...
    DPV0           (M, NC) float32 ...
    ... (more data variables) ...
    SBUUA          (avgtime, LSEDTOT, MC, N) float32 ...
    SBVVA          (avgtime, LSEDTOT, M, NC) float32 ...
    SSUUA          (avgtime, LSED, MC, N) float32 ...
    SSVVA          (avgtime, LSED, M, NC) float32 ...
Attributes:
    Conventions:  CF-1.6 SGRID-0.3
    institution:  Deltares
    references:   www.deltares.nl
    source:       Deltares, FLOW2D3D Version 6.02.08.6738, Nov  9 2016, 11:40:05
    history:      This file is created on 2019-10-28T13:27:00+0100, Delft3D
    LAYER_MODEL:  SIGMA-MODEL

Grid description

The NetCDF's grid variable looks like this

<xarray.DataArray 'grid' ()>
array(0, dtype=int32)
Attributes:
    units:                1
    cf_role:              grid_topology
    topology_dimension:   2
    node_dimensions:      MC NC
    face_dimensions:      M:MC (padding: low) N:NC (padding: low)
    face_coordinates:     XZ YZ
    node_coordinates:     XCOR YCOR
    vertical_dimensions:  SIG_LYR:SIG_INTF (padding: none)

MC is an array from 1 to 42. NC is an array from 1 to 61. XZ is a grid of X-coordinates of cell centres in m.

<xarray.Variable (M: 42, N: 62)>
array([[   0.,    0.,    0., ...,    0.,    0.,    0.],
       [   0.,    0.,    0., ...,   50.,   50.,   50.],
       [   0.,    0.,    0., ...,  150.,  150.,  150.],
       ...,
       [   0.,    0.,    0., ..., 3850., 3850., 3850.],
       [   0.,    0.,    0., ..., 3950., 3950., 3950.],
       [   0.,    0.,    0., ...,    0.,    0.,    0.]], dtype=float32)
Attributes:
    standard_name:  projection_x_coordinate
    long_name:      X-coordinate of cell centres
    units:          m

YZ is a grid of Y-coordinates of cell centres in m.

<xarray.Variable (M: 42, N: 62)>
array([[   0.,    0.,    0., ...,    0.,    0.,    0.],
       [   0.,    0.,    0., ..., 5850., 5950., 6050.],
       [   0.,    0.,    0., ..., 5850., 5950., 6050.],
       ...,
       [   0.,    0.,    0., ..., 5850., 5950., 6050.],
       [   0.,    0.,    0., ..., 5850., 5950., 6050.],
       [   0.,    0.,    0., ...,    0.,    0.,    0.]], dtype=float32)
Attributes:
    standard_name:  projection_y_coordinate
    long_name:      Y-coordinate of cell centres
    units:          m

SIG_LYR is an array of the sigma-layer thicknesses, which become a lot more detailed towards the bottom.

<xarray.DataArray 'SIG_LYR' (SIG_LYR: 80)>
array([-0.05   , -0.15   , -0.225  , -0.275  , -0.32   , -0.3515 , -0.3745 ,
       -0.3975 , -0.4205 , -0.4435 , -0.4665 , -0.4895 , -0.5125 , -0.5355 ,
       -0.5585 , -0.5815 , -0.6045 , -0.6275 , -0.6505 , -0.6735 , -0.693  ,
       -0.709  , -0.725  , -0.741  , -0.757  , -0.773  , -0.789  , -0.805  ,
       -0.821  , -0.837  , -0.849  , -0.857  , -0.865  , -0.873  , -0.881  ,
       -0.889  , -0.897  , -0.905  , -0.913  , -0.921  , -0.927  , -0.931  ,
       -0.935  , -0.939  , -0.943  , -0.947  , -0.951  , -0.955  , -0.959  ,
       -0.963  , -0.966  , -0.968  , -0.97   , -0.972  , -0.974  , -0.976  ,
       -0.978  , -0.98   , -0.982  , -0.984  , -0.9855 , -0.9865 , -0.9875 ,
       -0.9885 , -0.9895 , -0.9905 , -0.9915 , -0.9925 , -0.9935 , -0.9945 ,
       -0.99525, -0.99575, -0.99625, -0.99675, -0.99725, -0.99775, -0.99825,
       -0.99875, -0.99925, -0.99975], dtype=float32)
Coordinates:
  * SIG_LYR  (SIG_LYR) float32 -0.05 -0.15 -0.225 ... -0.99875 -0.99925 -0.99975
Attributes:
    standard_name:  ocean_sigma_coordinate
    long_name:      Sigma coordinates of layer centres
    units:          1
    formula_terms:  sigma: SIG_LYR eta: S1 depth: DPS

SIG_INTF is similar only with Sigma coordinates of layer interfaces.

DPS contains the bottom depth in m at each timestep, which I guess is needed for a vertical cross-section plot in m.

<xarray.DataArray 'DPS' (time: 49, M: 42, N: 62)>
[127596 values with dtype=float32]
Coordinates:
    XZ       (M, N) float32 ...
    YZ       (M, N) float32 ...
  * time     (time) datetime64[ns] 2019-05-02 ... 2019-05-02T04:00:00
Dimensions without coordinates: M, N
Attributes:
    long_name:  Bottom depth (zeta point)
    units:      m
    grid:       grid
    location:   face

Example of DataArray

This is an example of a variable in the dataset I would like to plot vertically on a sigma-grid

<xarray.DataArray 'RHO' (time: 49, KMAXOUT_RESTR: 80, M: 42, N: 62)>
[10207680 values with dtype=float32]
Coordinates:
    XZ             (M, N) float32 ...
    YZ             (M, N) float32 ...
  * KMAXOUT_RESTR  (KMAXOUT_RESTR) int32 0 1 2 3 4 5 6 ... 73 74 75 76 77 78 79
  * time           (time) datetime64[ns] 2019-05-02 ... 2019-05-02T04:00:00
Dimensions without coordinates: M, N
Attributes:
    long_name:  Density per layer in zeta point
    units:      kg/m3
    grid:       grid
    location:   face

Things I've tried

Load dataset directly with 'gridded'

When I try to load the dataset with gridded as follow 

ncfilename = '/Results/trim-GUIM.nc'
dataset = gridded.Dataset(ncfile=ncfilename)

I get this error 

ValueError: The netCDF file appears to conform to SGRID conventions, but padding values cannot be found.

Load only one variable into gridded

If I only load one variable into gridded like this

density = gridded.Variable.from_netCDF(filename=ncfilename, name='Density', varname='RHO')

I get this error

~.../gridded/gridded/pysgrid/sgrid.py in infer_location(self, variable)
    319         if len(variable.shape) < 2:
    320             return None
--> 321         difference = (shape[-2:] - self.node_lon.shape).tolist()
    322         if (difference == [1, 1] or difference == [-1, -1]) and self.center_lon is not None:
    323             location = 'center'

AttributeError: 'NoneType' object has no attribute 'shape'

I'm really lost and at this point I'm not even sure if gridded is the right tool to help me plot with sigma-layers. Any help is appreciated, even if anyone could point me to an example with an ocean sigma-grid I'd be really grateful.

ChrisBarker-NOAA commented 4 years ago

Sorry we weren't more helful here, this is indeed the kind of thing we want gridded to be able to do.

Did you find a solution?

JulesBlm commented 4 years ago

Hi, I tried debugging Gridded for a while and if I remember correctly the error was in the parse_padding function in read_netcdf.py, more specifically in the regex that parses the padding string. I couldn't figure it out so I ended up using xarray and writing my own plot methods to plot according to the sigma grid. These can be found here. The processing_2d module of gridded is very useful for processing Delft3D output.

ChrisBarker-NOAA commented 4 years ago

Thanks -- we really need to refactor the loading code -- it should have a cleaner distinction between parsing our the standards and loading the data -- making it easier to override the parsing part.

If you have a moderate size example, could you add it to this issue?

Thanks for your feedback.

JulesBlm commented 4 years ago

Sure thing. I couldn't add it to this issue directly because files are limited to 10MB here. So I uploaded to my repo here. I had to remove a bunch of data variables from the dataset to keep it under GitHub's 100MB file limit, but grid metadata and all that stuff is still there so it's a good example.

Thanks for the response

JulesBlm commented 4 years ago

I should mention that I'm 'abusing' Delft3D for geological modelling, so this output will not like a typical coastal model.

ChrisBarker-NOAA commented 4 years ago

great thanks! Hopefully we can look at this at some point.

And no problem with the "abusing" -- ideally gridded is supposed to be fairly general purpose.

-CHB