Apply for help about vertical interpolation in MPAS-ocean

[shipf16]

Dear Scientist,

Hello, everyone, I meet a technical problem with the vertical interpolation and apply for your guidance. For interpolating in-depth, I would like to compare and compute the difference in temperature and salinity between the ocean outputs from E3SM v2 and ocean reanalysis data (EN4.2.2). However, the vertical layers of MPAS are 60 layers and EN4.2.2 data has 42 layers. I am thinking whether I could convert the 60 layers of MPAS-ocean to 42 layers of EN4.2.2.

However, in the MPAS-ocean, I find that the depth coordinate varies with space and time. For example, in the layer 10, the depth varies from one grid cell to another. In this case, I am thinking about whether I could interpolate the MPAS-ocean outputs to the same depth from EN4 data. Are there any similar NCO commands like this (ncremap -P mwf -s grd_ocn -g grd_out) to perform the vertical interpolation prior to the horizontal interpolation? I want to compare the difference between MPAS-ocean outputs and EN4 data on each same depth. Could you please provide me with more guidance and help?

Thank you very much!

[xylar]

@shipf16, thanks for your question! I also received your email.

We don't have a great tool for doing vertical interpolation like you are requesting. What we do have is a python function in a conda package called mpas-analysis. Here is an example script showing how it can be used:

#!/usr/bin/env python

from mpas_analysis.shared.interpolation import interp_1d
import xarray as xr
import numpy as np


sim = \
    '20221129.CMPASO-JRA1p4_TL319_EC30to60E2r2.fix-time-series-stats.chrysalis'
# you could also iterate over years and months
dates = ['0001-01-01']

# a list of 3D variables you want to interpolate
variables = ['timeMonthly_avg_activeTracers_temperature',
             'timeMonthly_avg_activeTracers_salinity',
             'timeMonthly_avg_zMid']

# we often use a restart file as an MPAS mesh file, needed here for masking
mesh_filename = f'{sim}.mpaso.rst.0001-02-01_00000.nc'

# you will have to figure out how to read in the EN4 vertical grid as an
# xarray dataset and extract the (positive up) z coordinate

# note that the coordinate and the dimension need to be named differently
# (i.e. not following CF-conventions) because of a quirk in the interp_1d function,
# so z and nz in this example.
z = xr.DataArray(data=np.linspace(0, -5000, 42), dims=['nz'])

# now do the interpolation

# first, make a cell mask

ds_mesh = xr.open_dataset(mesh_filename)
# the index of the bottom cell, converted from Fortran to Python indexing
max_level_cell = ds_mesh.maxLevelCell - 1
n_vert_levels = ds_mesh.sizes['nVertLevels']

vert_index = xr.DataArray.from_dict({'dims': ('nVertLevels',),
                                         'data': np.arange(n_vert_levels)})

# make a mask for where cells are above the bathymetry
cell_mask = vert_index <= max_level_cell

# then, interpolate the requested variable for each date
for date in dates:
    in_monthly_filename = \
        f'{sim}.mpaso.hist.am.timeSeriesStatsMonthly.{date}.nc'
    out_monthly_filename = f'mpaso.T_S_to_EN4.{date}.nc'
    ds = xr.open_dataset(in_monthly_filename)

    # just keep the variables we want
    ds = ds[variables]
    
    # drop the Time dimension, since this interferes with the interpolation
    ds = ds.isel(Time=0)

    for variable in variables:
        # mask out invalid cells
        ds[variable] = ds[variable].where(cell_mask)

    # add the destination vertical coordinate
    ds['z'] = z

    # do the interpolation
    ds = interp_1d(ds, inInterpDim='nVertLevels',
                   inInterpCoord='timeMonthly_avg_zMid',
                   outInterpDim='nz', outInterpCoord='z')

    # the coordinate and dimension need to have different names for interp_1d
    # but can be renamed after interpolation to be the same if you like:
    ds = ds.rename({'nz': 'z'})

    ds.to_netcdf(out_monthly_filename)

It is probably a good idea to do this vertical interpolation before you use ncremap to do horizontal interpolation.

If you paste an ncdump -h on one of the EN4 files you're using, that could help me give you instructions on how to read in the vertical coordinate as needed in this example.

I should also note that this performs bilinear interpolation using the zMid coordinate in MPAS. It doesn't extrapolate above the middle of the top layer or below the middle of the bottom layer. Ideally, one would use conservative interpolation instead, accounting for the sum of overlaps between the bounds of each level on the MPAS vertical grid with the corresponding bounds of the EN4 vertical grid. But we don't have a tool to support that at this time.

[shipf16]

Thank you very much for your great help! I follow your nice guidance and keep the rest of codes unchanged including the for loop over dates. We could get the outputs from vertical interpolation. The script is shown as below and shared with others. Could you please help me check the script?

In the above results, I still couldn't understand why most regions do not have interpolated values in the first layer using Method 1. Could you please help me better understand this? Could you please tell me which Method you will recommend more for the vertical and horizontal interpolation based on the results of Method 1 and Method 2?

Again, many thanks for your great support and help!

#!/usr/bin/env python

from mpas_analysis.shared.interpolation import interp_1d
import xarray as xr
import numpy as np
import argparse

parser = argparse.ArgumentParser(
    description='Vertically interpolate MPAS-Ocean files',
    formatter_class=argparse.RawTextHelpFormatter)
parser.add_argument("-m", "--mpas_mesh", dest="mpas_mesh",
                    help="An MPAS-Ocean mesh file such as a restart file",
                    required=True)
parser.add_argument("-i", "--mpas_monthly", dest="mpas_monthly",
                    help="An MPAS-Ocean monthly average file to interpolate",
                    required=True)
parser.add_argument("-e", "--en4_file", dest="en4_file",
                    help="An file with the EN4 vertical coordinate",
                    required=True)
parser.add_argument("-o", "--out_file", dest="out_file",
                    help="The output file after vertical interpolation",
                    required=True)
args = parser.parse_args()

in_monthly_filename = args.mpas_monthly
mesh_filename = args.mpas_mesh
en4_filename = args.en4_file
out_monthly_filename = args.out_file

# you could also iterate over years and months
dates = ['0400-01-01']

# a list of 3D variables you want to interpolate
variables = ['timeMonthly_avg_activeTracers_temperature',
             'timeMonthly_avg_activeTracers_salinity',
             'timeMonthly_avg_zMid']
# you will have to figure out how to read in the EN4 vertical grid as an
# xarray dataset and extract the (positive up) z coordinate

# note that the coordinate and the dimension need to be named differently
# (i.e. not following CF-conventions) because of a quirk in the interp_1d function,

ds_en4 = xr.open_dataset(en4_filename)
z = -ds_en4.depth
z = z.drop('depth')

# now do the interpolation

# first, make a cell mask

ds_mesh = xr.open_dataset(mesh_filename)
# the index of the bottom cell, converted from Fortran to Python indexing
max_level_cell = ds_mesh.maxLevelCell - 1
n_vert_levels = ds_mesh.sizes['nVertLevels']

vert_index = xr.DataArray.from_dict({'dims': ('nVertLevels',),
                                         'data': np.arange(n_vert_levels)})

# make a mask for where cells are above the bathymetry
cell_mask = vert_index <= max_level_cell

# then, interpolate the requested variable for each date
for date in dates:
    ds = xr.open_dataset(in_monthly_filename)

    # just keep the variables we want
    ds = ds[variables]

    # drop the Time dimension, since this interferes with the interpolation
    ds = ds.isel(Time=0)

    for variable in variables:
        # mask out invalid cells
        ds[variable] = ds[variable].where(cell_mask)

    # add the destination vertical coordinate
    ds['z'] = z

    # do the interpolation
    ds = interp_1d(ds, inInterpDim='nVertLevels',
                   inInterpCoord='timeMonthly_avg_zMid',
                   outInterpDim='depth', outInterpCoord='z')

    # the coordinate and dimension need to have different names for interp_1d
    # but can be renamed after interpolation to be the same if you like:

    ds = ds.rename({'depth': 'z'})

    ds = ds.transpose('z', 'nCells')

    ds.to_netcdf(out_monthly_filename)

[xylar]

Thank you very much for your great help! I follow your nice guidance and keep the rest of codes unchanged including the for loop over dates.

I would take out dates and the for loop at this point because it isn't helpful. You aren't using date for anything. Other than that, the script looks good!

#!/usr/bin/env python

from mpas_analysis.shared.interpolation import interp_1d
import xarray as xr
import numpy as np
import argparse

parser = argparse.ArgumentParser(
    description='Vertically interpolate MPAS-Ocean files',
    formatter_class=argparse.RawTextHelpFormatter)
parser.add_argument("-m", "--mpas_mesh", dest="mpas_mesh",
                    help="An MPAS-Ocean mesh file such as a restart file",
                    required=True)
parser.add_argument("-i", "--mpas_monthly", dest="mpas_monthly",
                    help="An MPAS-Ocean monthly average file to interpolate",
                    required=True)
parser.add_argument("-e", "--en4_file", dest="en4_file",
                    help="An file with the EN4 vertical coordinate",
                    required=True)
parser.add_argument("-o", "--out_file", dest="out_file",
                    help="The output file after vertical interpolation",
                    required=True)
args = parser.parse_args()

in_monthly_filename = args.mpas_monthly
mesh_filename = args.mpas_mesh
en4_filename = args.en4_file
out_monthly_filename = args.out_file

# a list of 3D variables you want to interpolate
variables = ['timeMonthly_avg_activeTracers_temperature',
             'timeMonthly_avg_activeTracers_salinity',
             'timeMonthly_avg_zMid']
# you will have to figure out how to read in the EN4 vertical grid as an
# xarray dataset and extract the (positive up) z coordinate

# note that the coordinate and the dimension need to be named differently
# (i.e. not following CF-conventions) because of a quirk in the interp_1d function,

ds_en4 = xr.open_dataset(en4_filename)
z = -ds_en4.depth
z = z.drop('depth')

# now do the interpolation

# first, make a cell mask

ds_mesh = xr.open_dataset(mesh_filename)
# the index of the bottom cell, converted from Fortran to Python indexing
max_level_cell = ds_mesh.maxLevelCell - 1
n_vert_levels = ds_mesh.sizes['nVertLevels']

vert_index = xr.DataArray.from_dict({'dims': ('nVertLevels',),
                                         'data': np.arange(n_vert_levels)})

# make a mask for where cells are above the bathymetry
cell_mask = vert_index <= max_level_cell

# then, interpolate the requested variable for each date
ds = xr.open_dataset(in_monthly_filename)

# just keep the variables we want
ds = ds[variables]

# drop the Time dimension, since this interferes with the interpolation
ds = ds.isel(Time=0)

for variable in variables:
    # mask out invalid cells
    ds[variable] = ds[variable].where(cell_mask)

# add the destination vertical coordinate
ds['z'] = z

# do the interpolation
ds = interp_1d(ds, inInterpDim='nVertLevels',
               inInterpCoord='timeMonthly_avg_zMid',
               outInterpDim='depth', outInterpCoord='z')

# the coordinate and dimension need to have different names for interp_1d
# but can be renamed after interpolation to be the same if you like:

ds = ds.rename({'depth': 'z'})

ds = ds.transpose('z', 'nCells')

ds.to_netcdf(out_monthly_filename)

[xylar]

In the above results, I still couldn't understand why most regions do not have interpolated values in the first layer using Method 1. Could you please help me better understand this?

The depth of the top EN4 layer is at z = -5.02158976 m according to the data you printed out for me earlier. The interpolation approach used in method 1 will give NaN values if zMid for the top layer is below this top EN4 layer. It doesn't do any extrapolation. So depending on the sea-surface height, the middle of the top MPAS-Ocean layer will be above that top EN4 layer and at other locations, it will be below.

We could make a modified version of method 1 that extrapolates to the top of the ocean rather than cutting off at the middle of the top layer (and similarly it could extrapolate to the bottom of the bottom layer, rather than cutting off at its middle). We don't have this capability right now.

Could you please tell me which Method you will recommend more for the vertical and horizontal interpolation based on the results of Method 1 and Method 2?

This really depends on your needs. Method 1 is more accurate. You didn't happen to provide the details of method 2 so I can't speak to it in detail but it seems like that might be more useful to you because it seems like you are able to handle the top layer in the expected way. In my own work, I nearly always use method 1 because we also have vertical grids that are not anything like z-levels (e.g. when we include the ocean cavities under ice shelves around Antarctica) so it would not be accurate at all to use method 2. But you are using a simulation where zMid and the 1D depth coordinate are going to be rather close, as you have seen.

[shipf16]

Thank you very much for your detailed guidance. Many thanks for your great help and support! I hope this answer will also be helpful to other users to analyze the MPAS-ocean outputs in E3SM. Again, thank you very much!

[shipf16]

Hello, I am so sorry to bother you again. Today I check the ocean EN4.2.2 reanalysis data, and find that the depth in this data is not the exact depth. For example, for the first layer at 5.02 m, we could see that the first layer ranges from 0.00000 to 10.0475 meters using the software Panoply.

I am still worried that using Method 1, most of Atlantic ocean shows missing data in the top first layer in the above comparison. I plan to assimilate EN4.2.2 reanalysis data into the MPAS-ocean model in E3SM. I am concerned that the unused observed EN4.2.2 data in the first layer (5 m) may affect the performance of sea surface temperature.

I would like to provide the details of Method 2 and hope for your further guidance. For method 2, I plan to use this add_depth.py python script written by you to add the depth coordinates: python add_depth.py -c $depthfile -i $inputfile -o $outputfile. The $depthfile is mpas_refBottomDepth_60lyr.nc, which has information of the vertical coordinates at each point in space. After doing this, I could get the variable depth named reference depth of the center of each vertical level. I show the detailed information for the depth variable for your review.

I plan to do horizontal interpolation first to get the Grid data in latitude and longitude, and then do vertical interpolation. Could you please tell me whether I could do vertical 1D interpolation based on the new variable depth and how to write codes to interpolate the depth from MPAS-ocean to EN4.2.2 reanalysis data? I know the depth coordinate varies with space and time, but maybe changes in time could be neglected? Could you please provide me with more guidance and help?

Many thanks for your great help and support!

double depth(depth) ;
long_name = "reference depth of the center of each vertical level" ;
standard_name = "depth" ;
units = "meters" ;
axis = "Z" ;
positive = "down" ;
valid_min = 0. ;
valid_max = 5499.9902734375 ;
bounds = "depth_bnds" ;

depth = 5, 15, 25, 35, 45, 55, 65, 75, 85, 95, 105, 115, 125, 135, 145, 155,
165.0983984375, 175.47904296875, 186.29125, 197.66025390625,
209.71134765625, 222.57826171875, 236.40880859375, 251.37013671875,
267.65416015625, 285.48361328125, 305.1191796875, 326.86794921875,
351.0934375, 378.2275390625, 408.784609375, 443.37765625, 482.7366796875,
527.72796875, 579.3728515625, 638.862578125, 707.56328125, 787.0025,
878.825234375, 984.705859375, 1106.20421875, 1244.566953125,
1400.4971875, 1573.94640625, 1764.00328125, 1968.94421875, 2186.4565625,
2413.9715625, 2649.00125, 2889.3846875, 3133.4046875, 3379.79375,
3627.670625, 3876.4521875, 4125.7684375, 4375.3928125, 4625.190625,
4875.08375, 5125.0284375, 5375 ;

[xylar]

I plan to assimilate EN4.2.2 reanalysis data into the MPAS-ocean model in E3SM.

I see, I thought the purpose was for analysis. If the purposes is assimilating EN4 in MPAS-Ocean, you probably want to be interpolating in the other direction. You should be interpolating the EN4 data to the MPAS-Ocean depth. Since EN4 is on a regular grid, you just need to do interpolation first to the MPAS horizontal location (e.g. using ncremap from EN4 to MPAS, rather than the other way around) and then do simple 1D vertical interpolation from your z coordinate to MPAS-O's zMId.

I think the way you're doing it is the wrong way around for assimilating data into MPAS-Ocean. What you want to do is much easier than what I originally thought -- it's easy to get data from a regular grid in 3D to the MPAS-Ocean points. It's much harder to go the other direction.

Please let me know if what I just said is accurate. If so, do you need specific guidance in this direction? I may have time to write an example script tomrrow.

[shipf16]

Hello, thanks a lot for your quick reply. I am not familiar with MPAS-Ocean and apply for your further guidance. Could you please help me check if my thinking is right? I would like to assimilate EN4.2.2 reanalysis data into the MPAS-ocean model in E3SM using the 4DEnVar method. EN4 is the regular grid with latitude, longitude and depth. We know that we could first use ncremap to do horizontal interpolation for MPAS-Ocean outputs (ncremap -P mpas -m map_EC30to60E2r2_to_cmip6_180x360_aave.20220301.nc inputfile outputfile), and then we also get the regular MPAS-Ocean outputs (depth=60, lat=180, lon=360) with latitude and longitude grid. I have already finished the horizontal interpolation of EN4 to the same horizontal resolution (depth=42, lat=180, lon=360).

However, the problem is that MPAS-Ocean outputs and EN4 have different vertical layers, and the depth coordinate of MPAS-Ocean varies with space and time. For example, in the layer 10, the depth varies from one grid cell to another. So I am thinking to add the approximate depth using the add_depth.py command (python add_depth.py -c $depthfile -i $inputfile -o $outputfile). The $depthfile is mpas_refBottomDepth_60lyr.nc, which has information of the vertical coordinates at each point in space. After doing this, I could use the refBottomDepth variable in the current data file. So I plan to do 1D vertical interpolation of EN4 reanalysis data (depth=42, lat=180, lon=360) to the refBottomDepth of MPAS-Ocean, and I could compute the difference between EN4 and MPAS-ocean in different horizontal grids and vertical layers for data assimilation. Maybe I could also choose to vertically interpolate MPAS-O to the vertical levels of EN4 but Method 1 provides the missing data in the first top layer over most of the Atlantic ocean and makes me sad. For data assimilation, it is good to convert model data from model space to observation space in order not to lose the observation information. So it's better and I would like to keep the EN4 data in the regular grid (lat=180, lon=360).

Could you please tell me if the above plan is right? I know this is complicated and difficult work. I sincerely apply for your help and need specific guidance in this direction. If possible, it's better to help me write an example script. The interpolation is important for data assimilation, and I look forward to your further guidance. Many thanks for your kind guidance and help!!

[xylar]

@shipf16, it does seem like using a 1D depth coordinate is your best plan, then. I would have thought that bigger inaccuracies would come from interpolating model data to the EN4 grid and back but I don't know a ton about data assimilation.

the add_depth.py script should add a depth coordinate with depth_bnds. You should not use refBottomDepth as the depth coordinate, as that corresponds to depth_bnds[:, 1], the bottom of each layer rather than the middle. If you use depth instead, you should be okay. Other than that, your plan sounds okay to me. Certainly, some accuracy will be lost by using the 1D coordinate but I doubt this will be the primary source of error in your process.

Let me know how it goes. I will also keep you posted on whether we create a version of the more accurate vertical interpolation that extrapolates beyond the middle of the top and bottom layers.

[shipf16]

Hello, thank you very much for your reply. I notice that refBottomDepth is not the depth coordinate, and I should use the depth instead. Thank you very much for your reminder, and I avoid this mistake.

I apply for your further guidance about the add_depth.py script. The add_depth.py script adds the depth variable (depth=60) named reference depth of the center of each vertical level. Could you please tell me if this script vertically interpolates the variables into the same vertical level, or this script just provides me the approximate depth for each vertical layer? If this script just provides the approximate depth, could you please tell me if I could do 1D vertical interpolation of EN4 reanalysis data at this approximate depth, and then compute the difference between MPAS-ocean outputs and EN4 reanalysis in each vertical layer at this approximate depth, and then input this difference to data assimilation system?

I am still worried about this vertical interpolation. Method 1 provides the accurate depth for MPAS-O outputs, but I am so sorry to notice the missing data of the Atlantic ocean in this first layer after 1D vertical interpolation. The ocean's first layer has strong influence on the sea surface temperature. Could you please tell me if we could do some changes in the Method 1 script to extrapolate the first layer, especially where zMid for the top layer is below this top EN4 layer (maybe this way is not accurate, but really want to try)? If we could fill in the missing values, Method 1 is probably better the Method 2 in the accurate vertical layer.

I see, I thought the purpose was for analysis. If the purpose is assimilating EN4 in MPAS-Ocean, you probably want to be interpolating in the other direction. You should be interpolating the EN4 data to the MPAS-Ocean depth. Since EN4 is on a regular grid, you just need to do interpolation first to the MPAS horizontal location (e.g. using ncremap from EN4 to MPAS, rather than the other way around) and then do simple 1D vertical interpolation from your z coordinate to MPAS-O's zMId.

Could you please provide your opinion for interpolating the EN4 data to the MPAS-Ocean depth? The horizontal interpolation of converting MPAS-ocean outputs to the regular grid with latitude and longitude is easy to achieve. However, the vertical interpolation still bothers me and seems the big difficulty due to the time-varying, 3D depth zMid coordinate. Could you please give me more specific guidance in this direction or help me write an example script? Please accept my apologies for any inconvenience caused and thank you very much for your great help!

[xylar]

Could you please tell me if this script vertically interpolates the variables into the same vertical level, or this script just provides me the approximate depth for each vertical layer?

The script just provides the approximate vertical coordinate. All other data remains the same as before.

If this script just provides the approximate depth, could you please tell me if I could do 1D vertical interpolation of EN4 reanalysis data at this approximate depth, and then compute the difference between MPAS-ocean outputs and EN4 reanalysis in each vertical layer at this approximate depth, and then input this difference to data assimilation system?

Yes, you should be able to do that. There will be small errors because zMid differs slightly from depth but as I suggested above, those errors are likely small compared with other errors in the data assimilation process.

Method 1 provides the accurate depth for MPAS-O outputs, but I am so sorry to notice the missing data of the Atlantic ocean in this first layer after 1D vertical interpolation. The ocean's first layer has strong influence on the sea surface temperature. Could you please tell me if we could do some changes in the Method 1 script to extrapolate the first layer, especially where zMid for the top layer is below this top EN4 layer (maybe this way is not accurate, but really want to try)? If we could fill in the missing values, Method 1 is probably better the Method 2 in the accurate vertical layer.

I agree with you that Method 1 as currently available is not going to serve your needs. I do not currently have time to fix Method 1 in the way that you need or to instruct you on how to do it yourself. I do have an issue in MPAS-Tools, one of our repositories for this kind of work:
MPAS-Dev/MPAS-Tools#479
Please keep an eye on that issue for now and maybe I will be able to write a function to serve your needs in the next month or so.

In the meantime, I think Method 2 is your best bet.

Could you please provide your opinion for interpolating the EN4 data to the MPAS-Ocean depth? The horizontal interpolation of converting MPAS-ocean outputs to the regular grid with latitude and longitude is easy to achieve. However, the vertical interpolation still bothers me and seems the big difficulty due to the time-varying, 3D depth zMid coordinate. Could you please give me more specific guidance in this direction or help me write an example script? Please accept my apologies for any inconvenience caused and thank you very much for your great help!

If you were interested in interpolating from the regularly spaced EN4 grid to the MPAS coordinates, it is easy to first use ncremap to remap EN4 to MPAS cell centers. Then, you would iterate over cells, interpolating (e.g. with numpy.interp from the EN4 depth coordinate to the MPAS zMid coordinate at that cell). I would be willing to provide a script to do this but it takes a non-trivial amount of time to write such a thing and at this point I feel a little bit like we're running around in circles.

For now, can you see how far you can get with Method 2, and we can return to other possibilities if it becomes clear that that approach isn't good enough?

[shipf16]

Hello, thank you very much for your detailed guidance. I have kept an eye on that issue of MPAS-Tools in your link, and I look forward to your new development for vertical interpolation of MPAS-Ocean. Please don't hesitate to contact me if the new function is completed for this accurate vertical interpolation. I will update the new vertical interpolation with your help.

The vertical interpolation is very important for users to analyze the difference between MPAS-Ocean outputs and other ocean reanalysis data. I will follow your guidance to use Method 2 to add the approximate vertical coordinate, and compute the difference between MPAS-ocean and EN4 reanalysis on each vertical level according to our discussion above.

Thank you very much! Please accept my heartfelt gratitude for your kind help and guidance. I am so grateful to work with the most supportive associates. Many thanks for your great support and help!!!