oceanspy icon indicating copy to clipboard operation
oceanspy copied to clipboard

Grid coords for LLC data

Open asiddi24 opened this issue 2 years ago • 1 comments

  • OceanSpy version: 0.2.0

Description

While investigating budget closures in the ECCOv4r4 dataset on SciServer, I realized that the grid object that is imported for it is defined with the Xp1 and Yp1 defined on the right side of the grid center. This creates inaccuracies while using the grid object for computations with xgcm. Here is what the grid object looks like right now.

od.grid

<xgcm.Grid>
Z Axis (not periodic, boundary=None):
  * center   Z --> left
  * left     Zl --> center
  * outer    Zp1 --> center
  * right    Zu --> center
X Axis (not periodic, boundary=None):
  * center   X --> right
  * right    Xp1 --> center
time Axis (not periodic, boundary=None):
  * center   time_midp --> outer
  * outer    time --> center
Y Axis (not periodic, boundary=None):
  * center   Y --> right
  * right    Yp1 --> center

Here is how the dataset looks like.

<xarray.Dataset>
Dimensions:     (time: 312, Zl: 50, face: 13, Y: 90, X: 90, Z: 50, Xp1: 90,
                 Yp1: 90, time_midp: 311, Zp1: 51, Zu: 50, nv: 2)
Coordinates: (12/42)
    CS          (face, Y, X) float32 dask.array<chunksize=(13, 90, 90), meta=np.ndarray>
    Depth       (face, Y, X) float32 dask.array<chunksize=(13, 90, 90), meta=np.ndarray>
    PHrefC      (Z) float32 dask.array<chunksize=(50,), meta=np.ndarray>
    PHrefF      (Zp1) float32 dask.array<chunksize=(51,), meta=np.ndarray>
    SN          (face, Y, X) float32 dask.array<chunksize=(13, 90, 90), meta=np.ndarray>
    XC          (face, Y, X) float32 dask.array<chunksize=(13, 90, 90), meta=np.ndarray>
    ...          ...
    rAw         (face, Y, Xp1) float32 dask.array<chunksize=(13, 90, 90), meta=np.ndarray>
    rAz         (face, Yp1, Xp1) float32 dask.array<chunksize=(13, 90, 90), meta=np.ndarray>
  * time        (time) datetime64[ns] 1992-01-16T12:00:00 ... 2017-12-16
    time_bnds   (time, nv) datetime64[ns] dask.array<chunksize=(1, 2), meta=np.ndarray>
  * time_midp   (time_midp) datetime64[ns] 1992-01-31T12:00:00 ... 2017-12-01
    timestep    (time) int64 dask.array<chunksize=(1,), meta=np.ndarray>
Dimensions without coordinates: nv
Data variables: (12/30)
    ADVr_SLT    (time, Zl, face, Y, X) float32 dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray>
    ADVr_TH     (time, Zl, face, Y, X) float32 dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray>
    ADVx_SLT    (time, Z, face, Y, Xp1) float32 dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray>
    ADVx_TH     (time, Z, face, Y, Xp1) float32 dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray>
    ADVy_SLT    (time, Z, face, Yp1, X) float32 dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray>
    ADVy_TH     (time, Z, face, Yp1, X) float32 dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray>
    ...          ...
    UVELMASS    (time, Z, face, Y, Xp1) float32 dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray>
    VVELMASS    (time, Z, face, Yp1, X) float32 dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray>
    WVELMASS    (time, Zl, face, Y, X) float32 dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray>
    oceFWflx    (time, face, Y, X) float32 dask.array<chunksize=(1, 13, 90, 90), meta=np.ndarray>
    oceQsw      (time, face, Y, X) float32 dask.array<chunksize=(1, 13, 90, 90), meta=np.ndarray>
    oceSPtnd    (time, Z, face, Y, X) float32 dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray>

Unlike the other regional datasets available on SciServer, here the len of X and Xp1 is the same. So, using the same grid definition in the intake catalog generates the wrong grid object. This is how the correct grid object should look like with the Xp1 and Yp1 defined on the left side of center.

<xgcm.Grid>
X Axis (not periodic, boundary=None):
  * center   X --> left
  * left     Xp1 --> center
Y Axis (not periodic, boundary=None):
  * center   Y --> left
  * left     Yp1 --> center
Z Axis (not periodic, boundary=None):
  * center   Z --> left
  * left     Zl --> center
  * outer    Zp1 --> center
  * right    Zu --> center
time Axis (not periodic, boundary=None):
  * center   time --> inner
  * inner    time_midp --> center

I will create a PR about this and simply change the catalog entry for the ECCOv4r4 dataset. I'm guessing we'll have to do the same for the LLC4320 dataset too ? @Mikejmnez

@MaceKuailv this should affect the way you're working with the grid object as well, in case you're using the same one.

asiddi24 avatar Aug 30 '22 23:08 asiddi24

This is true! Nice catch @asiddi24 . The question now is if this issue is propagated to after transforming the dataset (after the cutout). Have you taken a look? This is because the data as originally stored (and read) cannot be used to compute budgets until after getting rid of the face dimension.

Mikejmnez avatar Aug 31 '22 13:08 Mikejmnez