Remapping 3D FESOM2 data with cdo

[JanStreffing]

Most researchers outside of AWI will interact with FESOM2 output after remapping to regular grid, e.g. via cdo.

For AWI-ESM-1-1-MR CMIP6 data this works fine. Here is the salinity vertically interpolated at 1000m depth and remapped to a 2° regular grid via:

cdo genycon,r180x91 -selname,salt -setgrid,/pool/data/AWICM/FESOM1/MESHES/core/griddes.nc weights_unstr_2_r180x91_so.nc
cdo -remap,r180x91,weights_unstr_2_r180x91_so.nc $in $out

I do the same for FESOM2 output, and the missval being 0 messed things up:

cdo genycon,r180x91 -selname,salt -setgrid,/work/ab0246/a270092/input/fesom2/core2/core2_griddes_nodes.nc
cdo -remap,r180x91,weights_unstr_2_r180x91_so.nc $in $out

I tired the obvious -setctomiss,0 and it mostly works. But it turns out, not all of the spurious values are exactly 0:

cdo genycon,r180x91 -selname,salt -setgrid,/work/ab0246/a270092/input/fesom2/core2/core2_griddes_nodes.nc
cdo  -remap,r180x91,weights_unstr_2_r180x91_so.nc -setctomiss,0 $in $out

For some points along the coasts there is still a spuriously low salt concentration. Any advice how to remap FESOM2 output better?

[koldunovn]

Can you try to set missing points before doing interpolation?

[JanStreffing]

I had the same idea, unfortunately, no difference.

[koldunovn]

Can you list the commands you use for first setting the missing values and then interpolate?

[JanStreffing]

for var in salt;
do
        cdo -setctomiss,0 ${outdir}/${var}_${tmpstr}.nc ${outdir}/${var}_${tmpstr}_nan.nc
        cdo genycon,r180x91 -selname,${var} -setgrid,$gridfile  ${outdir}/${var}_${tmpstr}_nan.nc ${outdir}/weights_unstr_2_r180x91_${var}.nc &
done
wait

for var in salt;
do
        cdo -L -splitlevel -remap,r180x91,weights_unstr_2_r180x91_${var}.nc -selname,${var} -setgrid,$gridfile ${var}_${tmpstr}_nan.nc ${var}_${tmpstr}_ &
done
wait

for var in salt;
do
        for lvl in "000010" "000100" "001000" "004000";
        do
                mv ${var}_${tmpstr}_${lvl}.nc ${var}_${tmpstr}_${lvl}_remap.nc &
        done
done
wait

for lvl in "000010" "000100" "001000" "004000";
do
        cdo chname,salt,so salt_${tmpstr}_${lvl}_remap.nc so_${tmpstr}_${lvl}.nc &
done
wait

for var in so
do
        for lvl in 10 100 1000 4000;
        do
                #cdo -L splitseas -setctomiss,0 -expr,so=(so<20.0)?0.0:so;  -yseasmean ${var}_${tmpstr}_$(printf "%06d" $lvl).nc $outdir/${var}_${model_name}_198912-201411_${lvl}m_ &
                cdo -L -splitseas -yseasmean ${var}_${tmpstr}_$(printf "%06d" $lvl).nc $outdir/${var}_${model_name}_198912-201411_${lvl}m_ &
        done
done
wait

[koldunovn]

Looks OK to me. One thing that might happen is that interpolation in cdo ignores the missing flag, and continue to use zeros. So you can try to replace missing values by nan with something like setmissval,nan.

[chrisdane]

Jan, can you provide the ok and not-ok versions with the same colorbars? Its difficult to see the actual problem.

[JanStreffing]

@koldunovn setmissval,nan seems to have no effect. I tried it before and after -setctomiss,0. I think this is expected, as the command only changes a missval that cdo already recognizes as such to a different value. Is is possible that the low but non zero values stem from using partial cells? As in total salt content of a partial cell being divided by the whole cell area?

@chrisdane The first plot from FESOM1 is ok, the last one from FESOM2 is not. They are on the same colorbar. I've outlined areas where the problem appears and uploaded the last plot with much dpi (click for larger version).

Compared to:

[chrisdane]

How does the data look like without interpolating? Maybe what you call not ok is just the model result? :D

[JanStreffing]

Luckily that's not the case. Here is the same plot with pyfesom2:

[JanStreffing]

@helgegoessling: Here are the minimal examples on Mistral:

FESOM1 grid file: /pool/data/AWICM/FESOM1/MESHES/core/griddes.nc FESOM1 native output: /mnt/lustre01/work/ab0246/a270092/postprocessing/climatologies/cmip6_raw/AWI-ESM-1-1-LR/so_first_timestep_lvl1000.nc FESOM1 remapped: /mnt/lustre01/work/ab0246/a270092/postprocessing/climatologies/cmip6_raw/AWI-ESM-1-1-LR/so_first_timestep_lvl1000_remap.nc

FESOM2: grid file: /mnt/lustre01/work/ab0246/a270092/input/fesom2/core2/core2_griddes_nodes.nc FESOM2: native output: /mnt/lustre01/work/ab0246/a270092/download/so_first_timestep_lvl1000.nc FESOM2: remapped: /mnt/lustre01/work/ab0246/a270092/download/so_first_timestep_lvl1000_remap.nc

[pgierz]

Not relevant for current production versions of FESOM, but there is a plan to make the output CDO compliant by default. At the moment one still needs to run the spherlab program once to generate a grid des file, that will then be circumvented.

We already have UGRID compliance, and are still waiting for tests, see here: https://github.com/FESOM/fesom2/pull/300

Once I have an OK from our test person, CDO is next on the list...

[helgegoessling]

I think @chrisdane might still be right, meaning that the issue is in the data, but this leaves the question why pyfesom does not show such funny coastal values.

Using the data from the minimal example and plotting it with yet another tool - spheRlab - I get the following when I plot the corresponding polygons for each node individually:

NAs are in cyan. We see funny values scattered along "coasts" (well, the 1000m "coast"), with salinities always around 27psu.

I thought that maybe pyfesom instead only plots the triangular elements, always averaging over the corresponding three nodal values, so I did the same, getting this:

Many of the funny values are now smeared out, but they remain at around 27psu where bigger patches of them exist.

If I'm not mistaken, it's thus an actual issue with the data, leaving the primary question what is happening in FESOM (could it have to do with partial bottom cells?) and the secondary question why the pyfesom plot hides those values... now, thinking aloud, is it possible that pyfesom treats partial bottom cells as NA, although the FESOM data do provide those funny values there?

[pgierz]

One spontaneous idea: it looks as if @JanStreffing is using something that (in the background maybe?) uses contourf. This does some form of, for lack of a better word, "smearing" or interpolation to make "nice" (or "un-nice") contours. It would be interesting to see what happens instead if you directly plot the values without interpolation, e.g. with something like pcolormesh. It could be, that it is mixing an extreme negative value (probably a default for land points) with the real value in the ocean, thus producing these brown patches.

[JanStreffing]

The same values can also be seen in ncview:

[pgierz]

~Is that an LGM mesh?~

Nevermind, just saw level 1000 in the title of the screenshot

[helgegoessling]

No, but at 1000m depth!

[JanStreffing]

no, it's present day, but 1000m depth

[pgierz]

You do a remap with cdo genycon,r180x91. Does the same problem occur with other remapping methods? (Funnily enough, there is a remapping discussion also occurring in another group)

[helgegoessling]

@pgierz , my plots above suggest it's in the data after all, so not a remapping issue.

[pgierz]

...I need a whiteboard for this, but: what happens in partial polyhedrons? However, I would expect that sort of problem to then occur at any coastline, not just Greenland and New Zealand??

[chrisdane]

It does @pgierz: Zoom in this figure you see very low salinity values (blueish) along all coastlines (color is mean over 3 nodes of triangles at original nod2d locations): Maybe its the intlevel,1000? Where is the raw model output salt.fesom.2049.nc @JanStreffing?

[JanStreffing]

Could be. though it would have to behave differently for FESOM2 with intlevel than for FESOM1. The file is here:/mnt/lustre01/work/ab0246/a270092/download/salt.fesom.2049.nc

[JanStreffing]

Ooooh, yes of course it behaves differently for intlevel. Because the values that are in the ground on the lower level are 0 instead of nan. And so when interpolating vertically, that's when these strange half way between 0 and physical values come in.

I'll try that asap.

[JanStreffing]

Indeed, simply placing -setctomiss,0 before the vertical interpolation fixed the issue. Thank you all!

Btw the reason it worked with pyfesom, is that pyfesom does the vertical interpolation itself. I'm guessing its configured to set 0=nan without the need to tell it to explicitly.

[pgierz]

So it is written down in case some one goes down this particular rabbit hole again: what was the final CDO command?

[helgegoessling]

Oh gosh, I assumed the "native output" was native output ;-) ... should have looked into the cdo history :confounded:

[helgegoessling]

So it is written down in case some one goes down this particular rabbit hole again: what was the final CDO command?

If you look at a level that is a model (output) level, then everything's fine. If however you need to look at an intermediate level for which vertical interpolation is required, it would be (example, all in one command): cdo remapcon,r180x91 -intlevel,${level} -setctomiss,0 -setgrid,${gridfile} ${ifile} ${ofile} ... with the order of setctomiss and intlevel being critical.

[chrisdane]

@JanStreffing: the last figure you showed might also be not fully correct? The land sea mask does not fit to 1000m depth. See e.g. Panama Strait? Is this due to ncview?

@helgegoessling: but this command sets potentially correct values, e.g. temperature = 0°C, to nan, right? Wouldn't this be wrong as already asked by @JanStreffing in #289?

[helgegoessling]

but this command sets potentially correct values, e.g. temperature = 0°C, to nan, right? Wouldn't this be wrong as already asked by @JanStreffing in https://github.com/FESOM/fesom2/issues/289?

Yes. That's why it would be better to put NANs into the FESOM output in the first place. But also it's extremely unlikely to get exact zeros for actual temperatures, so for the time being I guess we can live with it.

[pgierz]

The messed up Panama is likely due to the coarse interpolation of 2˚x2˚, there were similar interpolation "hiccups" back in the MPIOM days. (oh man, I am starting to sound like a grandpa)

[pgierz]

That's why it would be better to put NANs into the FESOM output in the first place.

yes, but unfortunately NaN is treated differently in different post-processing languages. My recommendation would therefore be to tell FESOM to write a very large negative number (e.g. -9e99) and then leverage the _FillValue attribute of NetCDF.

[koldunovn]

@chrisdane FESOM did not have exact 0 in 20 years. There is a small chance there will be one in the next 20 :)

[helgegoessling]

That's why it would be better to put NANs into the FESOM output in the first place.

yes, but unfortunately NaN is treated differently in different post-processing languages. My recommendation would therefore be to tell FESOM to write a very large negative number (e.g. -9e99) and then leverage the _FillValue attribute of NetCDF.

Well, that's quasi what I meant :-).