Evaluate Wind QC Procedures in Aircraft YAML File

[CoryMartin-NOAA]

An aircraft YAML file exists in GDASApp copied from UFO/ewok/jedi-gdas on June 28, 2022.

Before diving into replicating the PrepBUFR preprocessing for aircraft observations, we should first check that the QC procedures implemented in UFO match those being performed in GSI on PrepBUFR aircraft observations.

• [ ] Run GDASApp for first week of August 2021 for aircraft observations
• [ ] Compare EffectiveQC against GsiEffectiveQC for all variables (T, q, winds, etc.)
• [ ] Make any modifications to YAML as necessary

This issue will be used to track progress, including figures, statistics, and paths to files.

[nicholasesposito]

I've generated two images of EffectiveQC and GsiEffectiveQC for air_temperature. Next is to

1. Look at the log to see if anything is rejected
2. Plots QC and see if there is a reason for that red area in the Atlantic (and other red areas). Perhaps they are on a reject list?
3. Uncomment Pre-QC section of the yaml. Re-run and see what the results are.

[nicholasesposito]

Per the most recent post, the results of each of the next 3 steps are as follows:

1. Look at the log to see if anything is rejected QC Aircraft air_temperature: 8 missing values. QC Aircraft air_temperature: 16 out of bounds. QC Aircraft air_temperature: 37 rejected by first-guess check. QC Aircraft air_temperature: 120326 passed out of 120387 observations. QC Aircraft eastward_wind: 690 missing values. QC Aircraft eastward_wind: 4685 out of bounds. QC Aircraft eastward_wind: 149 rejected by first-guess check. QC Aircraft eastward_wind: 114863 passed out of 120387 observations. QC Aircraft northward_wind: 690 missing values. QC Aircraft northward_wind: 4685 out of bounds. QC Aircraft northward_wind: 65 rejected by first-guess check. QC Aircraft northward_wind: 114947 passed out of 120387 observations. QC Aircraft specific_humidity: 108054 missing values. QC Aircraft specific_humidity: 529 rejected by first-guess check. QC Aircraft specific_humidity: 11804 passed out of 120387 observations.

2. Plots QC and see if there is a reason for that red area in the Atlantic (and other red areas). Perhaps they are on a reject list? I checked all of the 14 values that preQC could be, and none of them were 37, 16, 8 (all taken from above for air_temperature) or any combination of the sums.

3. Uncomment Pre-QC section of the yaml. Re-run and see what the results are. Currently building. WIll update here.

[nicholasesposito]

I uncommented the PreQC sections. I had to rebuild GDASApp because in one of the Obsfunctions, it was looking for ObsValue/pressure (which is the new way) where it needed to look for MetaData/air_pressure (because we're using the old way). I updated that, re-built, and re-ran and got the following plots:

This passes the "eyeball" test. Next thing I'm going to do is a scatterplot to see if the data truly matches.

[nicholasesposito]

I plotted GsiEffectiveQC and EffectiveQC for air_temperature. I couldn't at first because len(GsiEQC) was 8 fewer. I noticed that 8 was also the number of missing values from earlier. I decided to plot entire arrays. The lengths were now the same so a plot was doable. The plot wasn't perfect. I did a point-to-point comparison of the arrays and 5863 of 120387 (4.8% of points) data points did not match up.

NOTE: there are 27 data points where GSIEffectiveQC = 0 and EffectiveQC = 19.

[PraveenKumar-NOAA]

Created a document listing all the Aircraft QC procedures implemented in UFO and GSI/prepBUFR: https://docs.google.com/document/d/1NsapgN9IT_j7OZLWh9PaOnJ6VwVMOaByT_qv_s-RmQI/edit?pli=1

[nicholasesposito]

Here are the jedi hofx and GsiHofX for air_temperature plots:

[nicholasesposito]

And here is the air_Pressure vs (GsiHofX - jediHofX) NOTE: updated to invert y-axis

[CoryMartin-NOAA]

Here is a histogram of @nicholasesposito 's test, while the above scatter looks worrisome, the histogram shows much closer agreement

Mean of JEDI - GSI (non bias corrected) is -0.001K and the std. dev. is 0.123K

[PraveenKumar-NOAA]

I have created plots of EffectiveQC and GsiEffectiveQC for both the eastward_wind and northward_wind, with commented PreQC section in yaml. These plots look like the air-temperature one that Nick created. Next step is to create and check these for uncommented PreQC section.

[nicholasesposito]

Here are the plots for Specific Humidity. The first four are how the code current is.

The following 4 plots are for what happens when the PreQC check is added into the yaml.

[CoryMartin-NOAA]

@nicholasesposito those air_pressure Y axis plots look incorrect. It looks like it is plotting GsiEffectiveQC-EffectiveQC and not GsiHofX-hofx...

[nicholasesposito]

@CoryMartin-NOAA Updated!

[CoryMartin-NOAA]

@nicholasesposito I'm a little curious as to why the scatter plots are different for the QC runs. Does the scatter plot code filter out data somehow? I would think the H(x) differences should be the same regardless of QC choices

[nicholasesposito]

@CoryMartin-NOAA The scatterplots do no filtering. If you'd like to take a look, they're here. They're quite simple: /scratch1/NCEPDEV/da/Nicholas.Esposito/gdas_single_test_hofx3d_orig/diags/scatterplot.py /scratch1/NCEPDEV/da/Nicholas.Esposito/gdas_single_test_hofx3d_orig/diags/scatterplot_P_vs_G-J.py

[CoryMartin-NOAA]

@nicholasesposito thanks, I see that. Very interesting, So I wonder why the scatter plots are so much different... The most minimum values differ by an order of magnitude.

[nicholasesposito]

I re-built the GDASApp, ran it, and plotted everything again. Results don't look too different.

air_temperature WITHOUT preqc

air_temperature WITH preqc

specific_humidity WITHOUT preqc

specific_humidity WITH preqc

[CoryMartin-NOAA]

Are the specific humidity plots correct? The # of points are equal but the GSI plots have them just mainly over CONUS whereas the JEDI plot has points globally...

[nicholasesposito]

I made a change in my plotting script. The data is right, but the number of datapoints was being calculated incorrectly. I updated the images above.

[PraveenKumar-NOAA]

I built and ran GDASApp for a yaml with uncommented PreQCs section and plotted EffectiveQC vs GsiEffectiveQC for the eastWind and northWind. Plots look visually similar, but below is full details of the actual number of missing values, out of bounds, rejected by first-guess check, and rejected by pre-QC reports w/o pre QC filter:

eastward_wind and northward_wind without pre-QC: _QC Aircraft eastward_wind: 690 missing values. QC Aircraft eastward_wind: 4685 out of bounds. QC Aircraft eastward_wind: 149 rejected by first-guess check. QC Aircraft eastward_wind: 114863 passed out of 120387 observations. QC Aircraft northward_wind: 690 missing values. QC Aircraft northward_wind: 4685 out of bounds. QC Aircraft northward_wind: 65 rejected by first-guess check. QC Aircraft northwardwind: 114947 passed out of 120387 observations.

eastward_wind and northward_wind with pre-QC: _QC Aircraft eastward_wind: 690 missing values. QC Aircraft eastward_wind: 5985 rejected by pre QC. QC Aircraft eastward_wind: 4296 out of bounds. QC Aircraft eastward_wind: 109416 passed out of 120387 observations. QC Aircraft northward_wind: 690 missing values. QC Aircraft northward_wind: 5985 rejected by pre QC. QC Aircraft northward_wind: 4296 out of bounds. QC Aircraft northwardwind: 109416 passed out of 120387 observations.

[nicholasesposito]

Using the sfc.yaml that Cory graciously wrote as a basis, I wrote an aircraft yaml that does similar things and ran it with ush/ufoeval/run_ufo_hofx_test.yaml. Most plots looked pretty good, and some that Cory and I discussed are below:

For this plot, ideally we would like a vertical line at 0, but the near-surface is quite scattered. It's likely a 'near surface' correction. It should be noted that numerous mountainous areas can have pressures into the 800s of hPas, and so that is likely why there are differences in those pressure ranges.

For this plot, the differences are VERY small (1e-8 kg/kg) so those differences are likely associated with machine-precision.

These plots are quite interesting. Northward wind has an offshoot to the toward negative Obs Values with increasing H(x) and appears to be mirroring the main data which shoots to the top-right (positive ObsValues) once ObsValues get to around 0. It's likely this is QC-related. Eastward wind instead has some data points that increase H(x) once ObsValues hit 0, and it's likely this is QC-related as well. Presumably, we will use QC to throw out these obs. For both plots, the Jedi H(x) vs Obs plots look the same.

Lastly, this is the HofX difference for air_temperature. We can't really say anything about it until the bias correction is finished, which is what Ben is working on.

[nicholasesposito]

Today, I ran it with Greg's QC that is already in there.

The histogram for all 4 variables are similar to this and look great.

Winds Eastward and Northward were similar in that some of the differences were very large, probably drowning out the smaller differences.

The air temperature hofxdiff histogram doesn't really mean anything until the bias correction is finished.

Specific humidity hofxdiff values are around 10e-8, which is fine. There were no measurable differences for the wind.

The GSI O-B vs JEDI O-B were interesting, as they showed the QC was working and many of the outer data points were discarded.

The comparisons of GSI HofX also allowed us to see that the QC was working, and many of the data points that were thought to be problematic are no longer there.

[nicholasesposito]

A few days ago, an option to apply near surface wind scaling for the VertInterp ObsOperator was added. The newest wind plots look much better. Starting off the errordiff vs Pressure did not change, per below:

The most striking difference was with the hofxdiff vs pressure. There is no longer the "tree" of data, with near-surface data points having large hofxdiff values. They are now all within about 1e-5 of each other.

The hotfx histogram reflects this.

The results are the same for windNorthward as well

[CoryMartin-NOAA]

@nicholasesposito fantastic! 1e-5 m/s is good enough for me!

[nicholasesposito]

The most recent task was to investigate the GSI setupw.f90 script and see what the QC performed for winds were. I found a small section, which was already in the aircraft.yaml, in addition to numerous other winds QC. I took out the additional QC from the yaml to determine what effect that would have on the results. What I found is that we had much better results when all of the QC remained in.

For example, you can see below that when comparing the Observations to GSI H(x) a number of data points to the top-left (for windNorthward) or vertically upward (for windEastward) of the plot have now returned, in comparison to keeping them all in (see the plots in the comment above).

The errordiff vs Pressure also got a little worse (R^2 increased from .0003 to .0009)

JEDI Omb vs GSI Omb plots don't look great either.

On the bright side, the hofxdiff are very small (10e-8) while errordiff looks fine.

[nicholasesposito]

The next thing to do was to see if the number of passing JEDI points (where EffectiveQC=0) and passing GSI points (where GsiEffectiveQC=0) are the same. They were not, but of all the tests performed so far, the yaml with only GSI QC performed the best. It was later discovered that of all the points for the GSIEffectiveQC=1 and EffectiveQC=0, the preQC = 1 at all 95 points where the difference was. To determine why the 95 was occurring, we also looked at GsiFinalObsError and found that all points were 1e8. We then looked at GsiUseFlag, and all values were -1, meaning they were not used. Looked at ObsType, and that was also pretty random.

At this point, we're not sure where the 95 differences are coming. There is a section of setupw.f90 that is a little unclear, so we will likely look into that in the future.

[nicholasesposito]

Recompiled GDASApp. Ran with old yamls and found results were different. Determined new results were better. Decided we want to determine why the 95 points are different. Now adding the gross check information into yaml. Went through prepobserrtable to add error information for pressure levels to ObsErrorModelStepwiseLinear@ObsFunction. Then added the satwindsspdbcheck using cermin/carmax from the convinfo table. This decreased the number of points where effqc=0 and gsieff=1 down to 52.