0

I've merged a 3D surface pressure field (ERA5, converted from Pa to hPa, function of lat,lon and time) with a 4D variable which is also a function of pressure levels (lat,lon,time,level).

So, my netcdf file has two fields, Temperature which is 4D:

float t(time, level, latitude, longitude)

surface pressure, which is 3d:

float sp(time, latitude, longitude)

The pressure dimension "level" is of course a vector:

int level(level)

What I want to do is make a mask for temperature for all locations where the pressure exceeds the surface pressure.

I know how to use nco to make a mask using a simple threshold:

 ncap2 -s 'mask=(level>800)' t_ps.nc mask.nc

But of course when I try to use the surface pressure

ncap2 -s 'mask=(level>sp)' t_ps.nc mask.nc

I get the error

ncap2: ERROR level and template sp share no dimensions

I think what I need to do is make a new variable like "level3d" which duplicates the pressure "level" to be a function of lat and lon, which I can then use to efficiently make the mask, yes? But I'm not sure how to do this with a dimension (I thought about cdo enlarge but couldn't get it to work).

By the way, instead of posting the data, this is the python api script I used to retrieve it

import cdsapi
c = cdsapi.Client()
c.retrieve(
    'reanalysis-era5-single-levels-monthly-means',
    {
        'format': 'netcdf',
        'product_type': 'monthly_averaged_reanalysis',
        'variable': 'surface_pressure',
        'year': '2020',
        'month': '03',
        'time': '00:00',
    },
    'ps.nc')

c.retrieve(
    'reanalysis-era5-pressure-levels-monthly-means',
    {
        'format': 'netcdf',
        'product_type': 'monthly_averaged_reanalysis',
        'variable': 'temperature',
        'pressure_level': [
            '1', '2', '3',
            '5', '7', '10',
            '20', '30', '50',
            '70', '100', '125',
            '150', '175', '200',
            '225', '250', '300',
            '350', '400', '450',
            '500', '550', '600',
            '650', '700', '750',
            '775', '800', '825',
            '850', '875', '900',
            '925', '950', '975',
            '1000',
        ],
        'year': '2020',
        'month': '03',
        'time': '00:00',
    },
    't.nc')
ClimateUnboxed
  • 7,106
  • 3
  • 41
  • 86

3 Answers3

1

If I understand the question properly, you should be able to do this using CDO, using either aexpr or expr along with clev (which should give you access to the pressure level). If one of the variables is a surface variable then CDO should just use it for calculations for all vertical levels when you mix and matched 3D/4D variables. NCO, I guess, is a lot stricter and requires a 100% identical data structure.

The following is likely to work:

cdo -L -aexpr,'mask=(clev(Temperature)>sp)' infile outfile
ClimateUnboxed
  • 7,106
  • 3
  • 41
  • 86
Robert Wilson
  • 3,192
  • 11
  • 19
  • Thanks Robert, I tried that but for some reason that is giving me a zero mask which is not a function of "level" : float mask(time, latitude, longitude). I'm not sure why it doesn't produce a 4D result? – ClimateUnboxed Jun 01 '23 at 12:23
  • That is strange. I'm not sure what could produce that behaviour. ERA5 should have *very* CF-compliant files, so CDO should behave as expected with them. Unless it was the merging that is causing the confusion. Was that done with CDO or NCO? – Robert Wilson Jun 01 '23 at 12:36
  • 1
    Ok, THANKS so much for your answer, super useful, combining it with enlarge, I managed to do what I wanted, a bit clunky but it works :) thanks again. ps: merging was done with cdo... – ClimateUnboxed Jun 01 '23 at 12:45
  • Hmm. Seems aexpr behaves slightly inconsistently with other CDO methods. Maybe one of those rare cases where using ``reduce_dims`` on the file will help, if the file with surface pressure has a level in it? – Robert Wilson Jun 02 '23 at 07:44
1

Your diagnosis of the NCO behavior is essentially correct. The "broadcast"

ncap2 -s 'mask=(level>sp)' t_ps.nc mask.nc

fails because level and sp are arrays (not scalars) that share no dimensions. The fix would be to create and use a temporary 3D version of level with something like

ncap2 -s 'level_3D[level,latitude,longitude]=level;mask=(level_3D>sp)' t_ps.nc mask.nc
Charlie Zender
  • 5,929
  • 14
  • 19
0

THANKS so much to Robert for pointing me to clev, with this I found the following solution using enlarge too. It is a bit clunky, maybe someone has a more streamlined solution? This works directly now on the two fields downloaded in my example api script.

# make 3d level field
cdo enlarge,t.nc -expr,"l=clev(t)" t.nc level.nc

# merge together
cdo -O  merge t.nc ps.nc level.nc t_ps_p.nc

# make mask, ###note here I have a 100 factor to change Pa->hPa
cdo setctomiss,0 -expr,'mask=(l*100)<sp' t_ps_p.nc mask.nc

# mask the data
cdo mul t.nc mask.nc masked_t.nc
ClimateUnboxed
  • 7,106
  • 3
  • 41
  • 86