How to extract info from multilayer NETCDF

Esteban82 · August 12, 2021, 2:24pm

Intro:

i am working with the GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc file which has 192 levels (z_levels) each for monthly solution (indicated as days since 01-01-2002). Important: NOT all month all included.

My problem:

I managed to get a map. On it I indicate the layer number but I would like to plot the date instead (month-year). How can I do it? Thanks in advance.

What I try:

1. gdalinfo

I run

gdalinfo NETCDF:"GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc":lwe_thickness

An I get the info of all the bands (below an extract of band 1).

Band 1 Block=720x1 Type=Float64, ColorInterp=Undefined
  NoData Value=-99999
  Unit Type: cm
  Metadata:
    comment=Coastline Resolution Improvement (CRI) filter is applied
    coordinates=time lat lon
    grid_mapping=WGS84
    long_name=Liquid_Water_Equivalent_Thickness
    NETCDF_DIM_time=106.5
    NETCDF_VARNAME=lwe_thickness
    standard_name=Liquid_Water_Equivalent_Thickness
    units=cm
    valid_max=778.7073956181606
    valid_min=-1719.85936245136
    _FillValue=-99999

There is a field with the information of time that I need (NETCDF_DIM_time). How can I get that value (106.5)?

2. grdinfo:
I try using grdinfo to extract the info. If I run
gmt grdinfo GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc?lwe_thickness -Q

I get all the z-levels (but I can find a way to get only the value corresponding to the layer that I want to plot).

GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: Title: JPL GRACE and GRACE-FO MASCON RL06Mv2 CRI
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: Command: GRACE and GRACE-FO JPL RL06Mv2-CRI
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: Remark: 
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: Gridline node registration used [Geographic grid]
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: Grid file format: nd = GMT netCDF format (64-bit float), CF-1.7
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: x_min: 0.25 x_max: 359.75 x_inc: 0.5 (30 min) name: longitude n_columns: 720
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: y_min: -89.75 y_max: 89.75 y_inc: 0.5 (30 min) name: latitude n_rows: 360
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: z_min: 106.5 z_max: 6924.5 z_inc: (variable) name: time [days since 2002-01-01T00:00:00Z] n_levels: 192
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: z_levels: 106.5, 129.5, 227.5, 258, 288.5, 319, 349.5, 380.5, 410, 439.5, 470, 495.5, 561.5, 592.5, 623, 653, 684, 714.5, 736.5, 777, 805.5, 836, 866.5, 897, 927.5, 958.5, 989, 1019.5, 1050, 1080.5, 1111.5, 1141, 1170.5, 1201, 1231.5, 1262, 1292.5, 1323.5, 1354, 1384.5, 1415, 1445.5, 1476.5, 1506, 1535.5, 1566, 1596.5, 1627, 1657.5, 1688.5, 1719, 1749.5, 1780, 1810.5, 1841.5, 1871, 1900.5, 1931, 1961.5, 1992, 2022.5, 2053.5, 2084, 2114.5, 2145, 2175.5, 2206.5, 2236.5, 2266.5, 2297, 2327.5, 2358, 2388.5, 2419.5, 2450, 2480.5, 2511, 2541.5, 2572.5, 2602, 2631.5, 2662, 2692.5, 2723, 2753.5, 2784.5, 2815, 2845.5, 2876, 2906.5, 2937.5, 2967, 2996.5, 3027, 3057.5, 3088, 3118.5, 3149.5, 3180, 3210.5, 3241, 3269.5, 3335.5, 3361.5, 3392, 3422.5, 3485.5, 3514.5, 3545, 3575.5, 3591.5, 3652, 3667.5, 3697.5, 3727.5, 3746.5, 3819, 3849.5, 3880.5, 3908, 3974.5, 4002.5, 4033.5, 4062, 4128, 4153.5, 4184, 4214.5, 4306.5, 4337, 4367.5, 4391.5, 4457.5, 4488, 4518.5, 4546, 4610.5, 4641, 4671.5, 4703, 4769.5, 4793, 4822.5, 4853, 4864, 4943.5, 4975.5, 5004.5, 5104.5, 5128.5, 5157, 5188.5, 5253, 5280, 5309.5, 5346.5, 5444.5, 5471.5, 5499, 5568.5, 5592.5, 5610.5, 5640, 6010, 6034, 6147.5, 6163, 6193.5, 6224.5, 6253, 6283.5, 6314, 6344.5, 6375, 6405.5, 6436.5, 6467, 6497.5, 6528, 6558.5, 6589.5, 6619.5, 6649.5, 6680, 6710.5, 6741, 6771.5, 6802.5, 6833, 6863.5, 6894, 6924.5
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: v_min: -89.75 v_max: 89.75 name: Liquid_Water_Equivalent_Thickness [cm]
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: scale_factor: 1 add_offset: 0
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: format: classic

3. grdinfo on layer
However. f I run grdinfo on a paritcular layer there is NO information about the date:
gmt grdinfo GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc?lwe_thickness[1]

GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: Title: JPL GRACE and GRACE-FO MASCON RL06Mv2 CRI
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: Command: GRACE and GRACE-FO JPL RL06Mv2-CRI
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: Remark: 
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: Gridline node registration used [Geographic grid]
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: Grid file format: nd = GMT netCDF format (64-bit float), CF-1.7
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: x_min: 0.25 x_max: 359.75 x_inc: 0.5 (30 min) name: longitude n_columns: 720
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: y_min: -89.75 y_max: 89.75 y_inc: 0.5 (30 min) name: latitude n_rows: 360
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: v_min: -89.75 v_max: 89.75 name: Liquid_Water_Equivalent_Thickness [cm]
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: scale_factor: 1 add_offset: 0
GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc: format: classic

PlanetGus · August 12, 2021, 2:36pm

So it is a discontinuous dataset ?

I guess you could use your solution 1 in combination with grep and awk ?

Esteban82 · August 12, 2021, 3:45pm

Yes, is discontinuous.

NC_GLOBAL#months_missing=2002-06;2002-07;2003-06;2011-01;2011-06;2012-05;2012-10;2013-03;2013-08;2013-09;2014-02;2014-07;2014-12;2015-06;2015-10;2015-11;2016-04;2016-09;2016-10;2017-02;2017-07;2017-08;2017-09;2017-10;2017-11;2017-12;2018-01;2018-02;2018-03;2018-04;2018-05;2018-08-2018-09

Thanks, I think grep and awk should work. I was thinking only in use gmt.

PlanetGus · August 12, 2021, 4:15pm

Well, if you know the list, you can probably pull something out using gmt math but at this point it’s almost easier to do it by hand.

Esteban82 · August 12, 2021, 4:17pm

Yes, I think on doing it by hand but an issue is that new data can be added every month. (I think right now there are 196 levels).

PlanetGus · August 12, 2021, 4:30pm

Maybe far fetched, but what about extracting the desired layer entirely. Like this the only info you’d get would be from this layer only (thus you read only band N)

Esteban82 · August 12, 2021, 5:08pm

Yes, that would be great but I can’t see/found the info of the date when I extract only one layer.

PlanetGus · August 12, 2021, 5:33pm

Even with grdcut ?

Esteban82 · August 12, 2021, 8:32pm

I am not sure if it is possible because is information of the header, not the data.

Esteban82 · August 12, 2021, 8:40pm

I found a workaround. The information can be found in the time_bounds subdataset.

Here a small sample (2nd column).

0.25    6924.5  6909
0.75    6924.5  6940
0.25    6894    6879
0.75    6894    6909
0.25    6863.5  6848
0.75    6863.5  6879

So I extract the data (grd2xyz) and the with gmt convert reverse the order (-I) and eliminate duplicates (-Td). With -qo(layer) I can get only the value that I need.

gmt grd2xyz GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc?time_bounds -o1 | gmt convert -I -Td -qo191

PlanetGus · August 12, 2021, 8:49pm

Is also have a workaround (maybe easier) :

gmt math NCFILE?"time" =

This extract all the time values.

You can tweak that with -fT to have absolute time using --TIME_UNIT and --TIME_EPOCH

Esteban82 · August 12, 2021, 9:59pm

Many thanks @PlanetGus. Your idea is better.

With:
$ gmt math GRCTellus.JPL.200204_202012.GLO.RL06M.MSCNv02CRI.nc?time = -fT --FORMAT_CLOCK_OUT=- --FORMAT_DATE_OUT=mm-yyyy

I get all the 192 dates exactly how I need them: