add ozone sonder plot based on 'develop aircraft' branch

[btang1]

I have added paired ozone sonder, and 2 plots (vertical time series & vertical boxplot) for ozone sonder.

this branch is based on 'develop aircraft branch'

[btang1]

@colin-harkins @rschwant @zmoon could you please give a review on this?

and I am continuing update some new features on this branch. for example (aggregate plot & density scatter plot)

[rschwant]

@zmoon and @btang1 is it possible to have this pull request just be a PR to the develop branch with only the changes you made Beiming to add the ozonesonde analysis & plots? When I look at the changed files to do the review, it includes all the files from the entire old develop_aircraft branch, which makes it hard to review? Also can you update the code to resolve the conflicts in the driver.py file.

[zmoon]

@rschwant seems like this line is necessary for aircraft plots to work, though it wasn't in develop

[rschwant]

I see, it got deleted with one of the satellite updates. Thanks for noticing this.

[zmoon]

@rschwant should be easier to review this now

[btang1]

Sorry for the obsolete code with aircraft, this is done at a very early age. Thanks zach for carrying me~

[rschwant]

Great job, this is looking really good! And it's great to have this capability added. I added some comments to make the code cleaner and also to use the functionality already available in the surface plots for providing the ozone min and max values and also for providing the model color types. It would be best to not have these namelist options be different for each plot as this is confusing to users. Let me know if you have questions. The part that Zach noted where the aircraft plots accidentally got deleted needs to be fixed this week because all the aircraft plots are broken now in develop. If you can address these comments this week, then it can go in with this PR. But if not no problem, I can submit a separate PR for this problem too.

[btang1]

Great job, this is looking really good! And it's great to have this capability added. I added some comments to make the code cleaner and also to use the functionality already available in the surface plots for providing the ozone min and max values and also for providing the model color types. It would be best to not have these namelist options be different for each plot as this is confusing to users. Let me know if you have questions. The part that Zach noted where the aircraft plots accidentally got deleted needs to be fixed this week because all the aircraft plots are broken now in develop. If you can address these comments this week, then it can go in with this PR. But if not no problem, I can submit a separate PR for this problem too.

thanks @rschwant I will try to address the changes (very helpful) and have a test run. I will let you know the results. if it works, I will update the PR. and then we can merge it.

[zmoon]

@rschwant we have an example now but only in YAML (runs with melodies-monet run). We could upload the figures and create a docs page for it, to display the figures. I could also upload the UFS-AQM data (unfortunately couldn't use your file, as the time didn't overlap with any of the GML ozonesondes) for the example data loading mechanism, like you did.

[rschwant]

Great yes that sounds good. Feel free to upload more UFS-AQM data too. I just wanted to limit it so that we would use less space.

[rschwant]

Great job addressing my comments so far!

I just have a question about setting up a jupyter notebook to connect your example to the ReadTheDocs completely.

See these here:
https://github.com/NOAA-CSL/MELODIES-MONET/tree/develop/docs/examples

Is there a missing commit maybe for this?

And then in your yaml file you comment out the code for an aggregate_plot_os, but since this is not available yet in the code can you just delete these lines instead just for clarity. Also the first line of the yaml file is missing a G in General.

Thanks, and we can talk about this in our upcoming meeting too!

[DWesl]

I'm trying to use this branch for stratospheric ozone and keep running into problems resulting from the model data being on pressure levels with a 1-D coordinate variable for model pressure. Is this an expected use-case, and am I missing something simple?

[btang1]

I'm trying to use this branch for stratospheric ozone and keep running into problems resulting from the model data being on pressure levels with a 1-D coordinate variable for model pressure. Is this an expected use-case, and am I missing something simple?
I am NOT exactly can solve you question. but we can look into some details of my code to make it more clear to be use.

1. the z-axis is based on pressure level
2. we designed and tested for surface to 20km around. so it should cover stratosphere.
3. make sure for model you have more than 1 input file, this code will do time interpolation from closest 2 time step.

hope that helps.

and by the way, what observation data and model data you used?

[DWesl]

I'm trying to use this branch for stratospheric ozone and keep running into problems resulting from the model data being on pressure levels with a 1-D coordinate variable for model pressure. Is this an expected use-case, and am I missing something simple?

I am NOT exactly can solve your question. but we can look into some details of my code to make it more clear to be use.

1. the z-axis is based on pressure level
2. we designed and tested for surface to 20km around. so it should cover stratosphere.
3. make sure for model you have more than 1 input file, this code will do time interpolation from closest 2 time step.

hope that helps.

and by the way, what observation data and model data you used?

Post-processed GFS O3MR 1-degree grib2 files for model data: my specific example is from my own run, but the operational output and the archive data here with inventory here should be close. I ran it through noaa-oar-arl/fv3grib2nc4, tried to follow examples for the configuration file, and ran into errors relating to pressure_model being a 1-D coordinate variable instead of a 3-D data variable.

Obs data is 100m SPO ozonesonde files from the GML archive

[btang1]

I'm trying to use this branch for stratospheric ozone and keep running into problems resulting from the model data being on pressure levels with a 1-D coordinate variable for model pressure. Is this an expected use-case, and am I missing something simple?

I am NOT exactly can solve your question. but we can look into some details of my code to make it more clear to be use.

1. the z-axis is based on pressure level
2. we designed and tested for surface to 20km around. so it should cover stratosphere.
3. make sure for model you have more than 1 input file, this code will do time interpolation from closest 2 time step.

hope that helps.
and by the way, what observation data and model data you used?

Post-processed GFS O3MR 1-degree grib2 files. My specific example is from my own run, but the operational output should be close. I ran it through noaa-oar-arl/fv3grib2nc4, tried to follow examples for the configuration file, and ran into errors relating to pressure_model being a 1-D coordinate variable instead of a 3-D data variable.

em. I think this could be the problem.

I the ozone sonder code we use 3D pressure variable for vertical interpolation. so the pressure variable we define in yaml file is 3D (lat, lon, vertical). our pressure is directly come from UFS-AQM from dyn.f* files. and if in your output file's pressure is 1-D I guess that could be the reason that yaml did not recognize it.

you mention using noaa-oar-arl/fv3grib2nc4 I wonder if it extract only 1D pressure or 3D pressure.

[DWesl]

That script seems to be a wrapper around wgrib2, which produces only a 1-D pressure variable when the data is on pressure levels.

From your description, the solution would seem to be adding a variable pressure_model to the output file that is horizontally constant?

[btang1]

That script seems to be a wrapper around wgrib2, which produces only a 1-D pressure variable when the data is on pressure levels.

From your description, the solution would seem to be adding a variable pressure_model to the output file that is horizontally constant?

yes. let is try that. and the pressure variable horizontal dimension need to match model horizontal dimension, try to see if that works.

[DWesl]

Okay, will try adding something like this to fv3grib2nc4 and run again

data = xr.open_dataset(out_file)
pressure_model, o3mr = xr.broadcast(data.coords["plevel"], data["O3MR"].isel(time=0))
pressure_model.name = "pressure_model"
pressure_model.to_netcdf(out_file, mode="a")

Thanks for the ideas

[btang1]

Okay, will try adding something like this to fv3grib2nc4 and run again

data = xr.open_dataset(out_file)
pressure_model, o3mr = xr.broadcast(data.coords["plevel"], data["O3MR"].isel(time=0))
pressure_model.name = "pressure_model"
pressure_model.to_netcdf(out_file, mode="a")

Thanks for the ideas

Any updates? is adding 3d pressure working for your case?

[DWesl]

3D pressure got the axes in the wrong order, so the vertical interpolation failed.
4D pressure got through vertical interpolation, and is now complaining that it can't add the old pressure coordinate to the new dataset.

I'm going to check whether renaming the current pressure coordinate plevel to z helps any. I'll need to change the loop where the renaming happens, but that's straightforward.

EDIT: Renaming plevel to z leads to it trying to interpolate z to the new pressure levels, which fails (it has a use for those values, but that's what pressure_model)
Renaming plevel to y leads to a different error about lat/lon mismatch.
Line 288 of melodies_monet/util/tools.py implies z is the right choice, so I'll work on figuring that out.

EDIT2: Told it to stop trying to interpolate plevel/z onto the pressure levels. It's now complaining about the observation index not being sorted (I think it's sorted by altitude, for whatever reason). After sorting the observation and model dataframes by time, I get an "After pairing" message, so I think that worked.

[btang1]

@zmoon looks like you finished the requirement @rschwant asked.
could we merge this now? @rschwant @zmoon

[rschwant]

Yes, I can approve this. @zmoon do you know why the docs are failing now? I think this is independent of the changes here.

[zmoon]

@rschwant there was a bug, which the CI caught for us, I think I've fixed it.

[rschwant]

Great, glad the CI found the bug. Already very useful! This is approved! Glad this is getting into the tool.

[btang1]

Great, glad the CI found the bug. Already very useful! This is approved! Glad this is getting into the tool.

thanks! @rschwant @zmoon it has been a long and wonderful journey!

I saw it been approved. has it been merged to develop yet?