How to use the ds.vertical() method?

[JairoVS]

Hello Robert! I have data on air pollutants concentrations through the atmosphere, divided into 47 vertical levels. I need to sum up the concentrations from 0 to 30 to work with only tropospheric data. I tried using the ds.vertical_sum(), but the sum is calculated for all 47 values. Is it possible to use any of the vertical methods to specify which levels I want to sum up and avoid getting the sum of all levels?

I highly appreciate any guidance you can provide on this!

[robertjwilson]

You would need to subset the levels first. Something like this would work:

levs = [ds.levels[0], ds.levels[29]]
ds.subset(levels = levs)
ds.run()
# check you've got 30 levels
ds.vertical_sum()

[JairoVS]

Thank you for your quick response! I will do so then.

[JairoVS]

Hello @robertjwilson! I tried subsetting the levels as suggested, but I failed. The dataset has 47 vertical levels associated with atmospheric pressure. Level 47 has the lowest value (3.8e-5, at the top of the atmosphere), and level 1 has the highest value (0.99, at the ground level).

I tried subsetting the levels using the index as you suggested "levs = [ds.levels[0], ds.levels[29]]" and also tried using the actual value of the levels "levs = [ds.levs[0.99], ds.levels[0.19]]" but it did not work.

And also tried switching the order of the values, but it did not work.
levs = [ds.levels[29], ds.levels[0]]
levs=[ds.levels[0.19], ds.levels[0.99]]

Again, thank you very much for any guidance you can provide me on how to make it work, and I apologize for bothering you with my questions.

P.S. I paste @here the list with the values of the levels, just in case it helps to solve the issue.

[3.8000005e-05,
0.000138674505,
0.00041406405,
0.00113378455,
0.0028636805,
0.0066847565,
0.014542271,
0.024079705,
0.03427161,
0.04828166,
0.06745011,
0.085439015,
0.10051436,
0.11825,
0.139115,
0.163661504,
0.192541017,
0.226513533,
0.26647905,
0.31279158,
0.3562501,
0.3937501,
0.4312501,
0.4687501,
0.5062501,
0.5437501,
0.58125015,
0.6187502,
0.65625015,
0.6875001,
0.7125001,
0.737500105,
0.762499965,
0.78750002,
0.8100002,
0.82750016,
0.842500125,
0.85750006,
0.87249996,
0.88749996,
0.90249991,
0.91749991,
0.93250006,
0.947499955,
0.962499776,
0.9774999,
0.992500024]

[robertjwilson]

I would need the data to figure out what is going wrong. Can you upload a file or point to where to download it? It sounds like something might be mispecified in the file.

[JairoVS]

Thank you very much @robertjwilson!

I add a file with a subsetted region in this comment. Thank you so much for letting me know if you can correctly open it or need to make any changes. As you will see, there are concentrations for two species in 47 layers for every lat,lon in the dataset.

I very much appreciate your time and support!

mxc_jul19.zip

[robertjwilson]

I see the problem. nctoolkit uses CDO to subset this, but the CDO call only takes integer values, so it won't work very well for your data as it goes from 0-1. I'll see if I can fix that or just ask the CDO developers to adjust their method.

In the meantime you can use the following by calling CDO internally. Note: the level order displayed by ds.levels is the reverse of what's in the file raw data. This is automatically ordered by nctoolkit, because in general that's easier for users. But this is a rare case where it might be confusing. I might need a warning to flag it's been re-ordered.

ds.cdo_command(f"sellevidx,{47-29}/47")
ds.run()
# check you've got 30 levels
ds.vertical_sum()

[JairoVS]

Thank you Robert! I tried what you said, and it worked! But, when calling the method using "{47-29}/47", I got data from top layer 47 to 18 (29 layers). I had to do it "{31-1}/1". That way, I got data from layers 30 to 1 (ground level). Maybe the levels are assigned as in the raw data when using CDO?

The issue is a bit more complex because for every lat and lot I will have to extract data from ground level up to a different top level (depending on the height of the troposphere for that specific grid). I have a NETCDF with the corresponding top layer for every lat and lot, with the same resolution as the datasets, do you have any recommendation on how could I perform the extraction using a specific top value for every lat,lon?

Thank you again for your time and support!

[robertjwilson]

I think in that case would need to create a Lon/lat/level mask file so that everything outside your condition is a missing value. You could do that using assign. Then multiply your original files with that mask. Without knowing what you condition is it’s hard to provide a specific solution.

[JairoVS]

Thank you for your reply @robertjwilson. I will have to calculate the concentrations in every grid by summing up the values from the ground to the troposphere.

I have attached a file with the troposphere level (TropLev variable in the file) for every grid in the mxc_jul19 file. For example, at (-100.9375, 17.75), the TropLev is 34.857594, meaning that for that lon, lat I will have to sum the concentration from level 1 to level 35. And the same process for every grid.

I am not sure if I explained myself well. An additional question, is it possible to round variables using lambda functions? I checked the list of accepted functions but did not find one for rounding values. I appreciate any guidance you can provide me on this.

Thank you so very much for all your time and support!

TropLev_MXC.zip