Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

YAML Configuration | Radiosondes #78

Closed
4 tasks
danholdaway opened this issue May 12, 2022 · 48 comments
Closed
4 tasks

YAML Configuration | Radiosondes #78

danholdaway opened this issue May 12, 2022 · 48 comments
Assignees
Labels
ufo done Observation operator evaluation in UFO completer waiting on jedi merge

Comments

@danholdaway
Copy link
Contributor

danholdaway commented May 12, 2022

Agreement with GSI using GeoVaLs in UFO (Hofx)

  • Done

Agreement with GSI using GeoVaLs in UFO (Quality control)

  • Done

Agreement with GSI using Swell

  • Done

Current issues:

@gmao-msienkie gmao-msienkie self-assigned this May 18, 2022
@gmao-msienkie
Copy link
Contributor

Original yaml file fails to be processed - the input data is missing 'LaunchTime' values. If a station has more than one launch in a synoptic time, the launch time would be used to group the obs from each launch separately. In the GSI the obs error inflation when multiple obs exist inside vertical model layers is done report-by-report, hence, sounding by sounding. Without launch time the multiple launches from a station would be grouped together for that filter.

@danholdaway
Copy link
Contributor Author

sondes_gsihofxbc_vs_hofx_air_temperature
sondes_gsihofxbc_vs_hofx_eastward_wind
sondes_gsihofxbc_vs_hofx_northward_wind
sondes_gsihofxbc_vs_hofx_specific_humidity
sondes_gsihofxbc_vs_hofx_surface_pressure
sondes_gsihofxbc_vs_hofx_virtual_temperature

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Nov 1, 2022

corscat_eferr_vs_finaler
Observational errors in JEDI and GSI. There is no duplication inflation in JEDI yet. A question, how to re-set observational error range [100-300 Pa] after the ObsErrorFactorSfcPressure@ObsFunction function? The following bounds check filters after it have no impact. cory.r.martin@noaa.gov

  # error inflation based on Ps correction
  - filter: Perform Action
    filter variables:
    - name: surface_pressure
    action:
      name: inflate error
      inflation variable:
        name: ObsErrorFactorSfcPressure@ObsFunction
        options:
          geovar_geomz: geopotential_height
          geovar_sfc_geomz: surface_geopotential_height
  # ObsError >= 100
  - filter: Bounds Check
    filter variables:
    - name: surface_pressure@ObsErrorData
    where:
    - variable: surface_pressure@ObsErrorData
      maxvalue: 100
    action:
      name: assign error
      error parameter: 100
  # ObsError <= 300
  - filter: Bounds Check
    filter variables:
    - name: surface_pressure@ObsErrorData
    where:
    - variable: surface_pressure@ObsErrorData
      minvalue: 300
    action:
      name: assign error
      error parameter: 300

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Nov 1, 2022

I tried another way to re-set the obs error range:

  # ObsError >= 100
  - filter: Variable Assignment
    where:
    - variable:
        #name: surface_pressure@ObsError
        name: ObsErrorData/surface_pressure
      maxvalue: 100
    assignments:
    #- name: surface_pressure@ObsError
    - name: ObsErrorData/surface_pressure
      value: 100
  # ObsError <= 300
  - filter: Variable Assignment
    where:
    - variable:
        #name: surface_pressure@ObsError
        name: ObsErrorData/surface_pressure
      minvalue: 300
    assignments:
    #- name: surface_pressure@ObsError
    - name: ObsErrorData/surface_pressure
      value: 300

However, JEDI this time thinks "ObsErrorData" is a new variable and asks an input of data 'type'. However, the UFO document says it is the updated obs error:
https://jointcenterforsatellitedataassimilation-jedi-docs.readthedocs-hosted.com/en/latest/inside/jedi-components/ufo/qcfilters/introduction.html#observation-errors

Anyway, Here are steps to handle a surface pressure error in GSI:
1, Assign an initial error (from the error table)
2, ObsErrorFactorSfcPressure@ObsFunction
3, Reset ObsError between Min ~ Max
4, Gross error check: (O-F)/ObsError
5, Inflate the error, which is the value before Step 3, because of duplication, and make it the final obs error.

@gmao-jjin3
Copy link
Contributor

corscat_omb_vs_gsiombbc
These O-B in UFO are close to those in GSI except for a few spots, likely because of some missing QC.

@gmao-jjin3
Copy link
Contributor

Here is an issue to re-set ObsError within a range after initial inflation.
https://github.com/JCSDA-internal/ufo/issues/2363

@gmao-jjin3
Copy link
Contributor

corscat_eferr_vs_finaler
Now the observational errors can be re-bounded correctly after initial inflation. See Yaml settings in https://github.com/JCSDA-internal/ufo/issues/2363.

@gmao-jjin3
Copy link
Contributor

corscat_qc_omb_vs_gsiombbc
Now UFO can reproduce GSI surface pressure Hofx after a bug fix in surface pressure correction. Sondes' "height", not "station_elevation" should be used in the correction. See the UFO pull-request: https://github.com/JCSDA-internal/ufo/pull/2493

@danholdaway
Copy link
Contributor Author

Great stuff @gmao-jjin3! Thanks! When making these scatter plots can you make them for all observations rather than what passed QC? There should be a match for all obs regardless of whether each system used the observation or not.

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Dec 1, 2022

@danholdaway Now I can see that the combining programs within the original IODA-converter package messed up observations. Here are surface_pressure and virtual_temperature in those separated and combined files.

For_others/conv_obs_check% python ck_ps.py
/discover/nobackup/drholdaw/JediData/GMAOReferenceRun/geos_to_r2d2/observations/conventional/obs_sep/sondes_ps_obs_2020121500.nc4
--Obs with z=3650m:
i,type,subtype,sta_id,lat,lon,z(3650m),ps: 203 120 0 b'55591   ' 29.67 91.13 3650.0 65110.0
For_others/conv_obs_check% python ck_tv.py
/discover/nobackup/drholdaw/JediData/GMAOReferenceRun/geos_to_r2d2/observations/conventional/obs_sep/sondes_tv_obs_2020121500.nc4
--Obs with z=3650m:
i,type,subtype,sta_id,lat,lon,z(3650m),tv: 6570 120 0 b'55591   ' 29.67 91.13 3650.0 269.85
For_others/conv_obs_check% python ck_tv_combined.py
/discover/nobackup/drholdaw/JediData/GMAOReferenceRun/geos_to_r2d2/observations/conventional/obs/sondes_obs_.nc4
--Obs with z=3650m:
i,type,subtype,sta_id,lat,lon,z(3650m),tv,ps: 28708 -2147483647 -2147483647 b'55591   ' 29.67 91.13 3650.0 **9.96921e+36 65110.0**
i,type,subtype,sta_id,lat,lon,z(3650m),tv,ps: 28710 120 0 b'55591   ' 29.67 91.13 3650.0 **269.85 9.96921e+36**

@danholdaway
Copy link
Contributor Author

Is this what is leading to the difference between GSI and UFO? Good to know if you've found the cause. Do you have a sense for what is going wrong in the combine tool? This is coming from ioda converters.

@gmao-jjin3
Copy link
Contributor

It is a reason. Not sure what is wrong in the combine Python programs at the moment.

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Dec 2, 2022

Okay. There is a very tiny difference in "air_pressure" in separated conventional IODA files. This values is saved in hPa in GSI' surface_pressure output and is multiplied by 100 in the IODA converter, while it is saved in Pa in outputs for surface temperature. They are supposed to be same in converted IODA files. However, the combine Python program doesn't think they are same because of the tiny difference! As a result, the combine program messed up temperature and pressure at this location.

For_others/conv_obs_check% python  ck_tv.py 
/discover/nobackup/drholdaw/JediData/GMAOReferenceRun/geos_to_r2d2/observations/conventional/obs_sep/sondes_tv_obs_2020121500.nc4
--Obs with z=3650m:
air_pressure 65110.0
i,type,subtype,sta_id,lat,lon,z(3650m),tv: 6570 120 0 b'55591   ' 29.67 91.13 3650.0 269.85
For_others/conv_obs_check% python ck_ps.py
/discover/nobackup/drholdaw/JediData/GMAOReferenceRun/geos_to_r2d2/observations/conventional/obs_sep/sondes_ps_obs_2020121500.nc4
--Obs with z=3650m:
air_pressure 65109.996
i,type,subtype,sta_id,lat,lon,z(3650m),ps: 203 120 0 b'55591   ' 29.67 91.13 3650.0 65110.0
For_others/conv_obs_check% python ck_tv_combined.py 
/discover/nobackup/drholdaw/JediData/GMAOReferenceRun/geos_to_r2d2/observations/conventional/obs/sondes_obs_.nc4
--Obs with z=3650m:
air_pressure 65109.996
i,type,subtype,sta_id,lat,lon,z(3650m),tv,ps: 28708 -2147483647 -2147483647 b'55591   ' 29.67 91.13 3650.0 9.96921e+36 65110.0
air_pressure 65110.0
i,type,subtype,sta_id,lat,lon,z(3650m),tv,ps: 28710 120 0 b'55591   ' 29.67 91.13 3650.0 269.85 9.96921e+36

@danholdaway
Copy link
Contributor Author

OK. So the combine program needs to be fitted with some tolerance for deciding if the obs are the same then. Or perhaps using some smarter logic to decide if the obs are the same. Tagging @CoryMartin-NOAA. I believe EMC worked on this combine program so may have a quicker idea of how to modify it.

@CoryMartin-NOAA
Copy link

@danholdaway I think this will be difficult to do because it is using (if I recall correctly) a python method that uses exact values to determine if it is the same location. I'm not sure how easy this will be to implement.

@danholdaway
Copy link
Contributor Author

OK, thanks for letting us know. I don't think we'll be able to replicate the GSI without the combine step since there are filters requiring multiple variables so we'll have to think on it. For example I think temperature is used on one of the surface pressure filters for example. Thanks for the work to establish what's happening Jianjun. We'll have to think a bit on the correct solution.

@CoryMartin-NOAA
Copy link

We can probably only check lat/lon/time and ignore pressure?

@gmao-jjin3
Copy link
Contributor

Noop, we cannot ignore pressure because sonde data are profiles. It may work for surface observations but we very likely get into troubles for other observations. I am re-making those IODA data after changing GSI. However, it is hard to guarantee there is no machine errors like these and it still is a risk not to have any tolerance in co-locating data even my new data work.

@gmao-jjin3
Copy link
Contributor

corscat_gsihofxbc_vs_hofx
This figure shows JEDI HofX are identical to GSI Hofx before QC. It looks like that there is no issue with my new data. However, we don't really solve the combine issue.

@gmao-jjin3
Copy link
Contributor

corscat_eferr_x_GsiDup_vs_finaler

This figure shows UFO surface pressure errors timed by those GSI duplicate factors are equal to GSI final observation errors. I think this figure confirms that surface pressure errors are set up correctly before we finalize the duplicate factor function.

@gmao-jjin3
Copy link
Contributor

corscat_EffectiveQC_vs_GsiEffectiveQC

@danholdaway This figure shows GSI effective QC flags vs UFO effective QC flags for all of the data. It puzzles me because that I find those GSI and UFO flags have the same number of zero values in the output file /discover/nobackup/jjin3/jedi/JediUfoTests/Output/sondes_qc_2020121500_0000.nc4, which is not what this figure shows. This figure is made with configuration between lines 341-365 in /discover/nobackup/jjin3/jedi/JediUfoTests/Eva/Config/sondes.yaml. My short Python program to check those GSI and UFO QC flags is /discover/nobackup/jjin3/jedi/For_others/conv_obs_check/ioda_sondes_ps/ck_ps_qc.py . Can you do me a favor to take a look at the yaml file and the python program and see why they show different results? Thanks.

@gmao-jjin3
Copy link
Contributor

@danholdaway Never mind. I think I need remove undefined values during plotting.

@gmao-jjin3
Copy link
Contributor

corscat_EffectiveQC_vs_GsiEffectiveQC
This figure suggests that there are same numbers zero values in GSI effective QC flags and UFO effective QC flags.

@gmao-jjin3
Copy link
Contributor

density_qc_omb_vs_gsiombbc

Histograms of surface pressure O-F after QC. They are identical in UFO and GSI.

@gmao-jjin3
Copy link
Contributor

@CoryMartin-NOAA I agree. I don't see a solution for this.

@gmao-msienkie
Copy link
Contributor

I was looking at some output that I had produced for the sondes running the ObsErrorFactorConventional filter and came across a sounding from station 50953 that looks like it was duplicated. The obs are slightly different on some levels though. In any case the two soundings were processed separately through errormod, but were treated as all one sounding in UFO by ObsErrorFactorConventional because they had the same launch time (time of first sounding level was the same for both, though they had different DHR sounding times). In any case, that's another reason for the output from UFO not matching GSI.
SID=50953_duplicate_sounding.txt

@gmao-msienkie
Copy link
Contributor

I made a list of the ADPUPA data in the gdas.prepbufr file from 20201215/00z and found there were quite a few (87!) duplicated soundings on that date. (From the station IDs they appear to be Chinese sondes.) Possibly some of the soundings were transmitted both in BUFR and text format. In any case, the GSI 'errormod' filter would treat the soundings separately while the UFO filter would filter the soundings together if the reported release time were the same. While the DHR values given for the sounding pairs were different, the drift time for the first level of the soundings was the same. Actually the drift lat,lon,time was the same for quite a few levels (if not most) and so would trigger down-weighting by the duplicate check.
dup_sonde_adpupa.txt

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Feb 8, 2023

density_qc_omb_vs_gsiombbc b
density_qc_omb_vs_gsiombbc

corscat_eferr_x_GsiDup_vs_finaler

corscat_eferr_vs_finaler

UFO O-F and observational errors for airTemperature match GSI values when observational errors made by UFO's function "ObsErrorFactorConventional" are replaced by GSI adjusted errors are used.

@gmao-jjin3
Copy link
Contributor

density_qc_omb_vs_gsiombbc b
density_qc_omb_vs_gsiombbc
corscat_eferr_x_GsiDup_vs_finaler

UFO O-F and observational errors for virtualTemperature do not match those values in GSI.

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Feb 8, 2023

Jianjun,

Myself and Hui Shao reviewed the GSI code quite a lot, yes.
Using error inflation as a means to deal with too many data points in a single grid volume is very questionable as a technique instead of alternate means, but, yes, we're aware of the GSI code for doing error inflation. We are trying our best to match GSI, but ultimately, that is not our goal. Rather, we desire to make the best DA use possible. We strive to make a first attempt to match GSI. However, this is not always guaranteed as sufficient time to produce an identical match is difficult.

Please feel free to open a GitHub/Zenhub issue if things are not optimal.

Greg

On Tue, Feb 7, 2023 at 2:03 PM Jin, Jianjun (GSFC-610.1)[SCIENCE SYSTEMS AND APPLICATIONS INC] <jianjun.jin@nasa.gov> wrote:
Hi Cory and Greg

Are you aware of any UFO codes to deal with temperature obs error inflation in GSI’s setupt.f90:
https://github.com/CoryMartin-NOAA/GSI/blob/ea04c69c4a3c29fc2109b6d7cdd21dfd049b1cf0/src/gsi/setupt.f90#L819

Thank you,

Jianjun

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Feb 13, 2023

density_qc_omb_vs_gsiombbc b
density_qc_omb_vs_gsiombbc
corscat_eferr_x_GsiDup_vs_finaler

corscat_eferr_vs_finaler

Now UFO O-F and observational errors match GSI values using (1) Wei's pressure-check inflator, and (2) GSI adjust error in Yaml configuration.

@gmao-jjin3
Copy link
Contributor

Many thanks to @gmao-wgu for writing the pressure check function. Here are a few updates in my directory: /discover/nobackup/jjin3/jedi/jedi-work2/Wei_Gu/ioda_v3_temperature

  1. Variable names were updated as in ioda_v3.
  2. Added the function for temperature observations as in setupt.f90 in GSI.
  3. Removed "GSI adjust error" from the function. We use it as a temporary solution and we can configure it in YAML files, so that we can remove it easily in future. BTW, I also copied my configuration files into that directory.

@danholdaway
Copy link
Contributor Author

Thanks for uploading these plots @gmao-jjin3. In each case can you also provide the plots without the multiplication by GsiDup? For posterity it would be good to just get a sense of how many observations this is impacting.

@gmao-jjin3
Copy link
Contributor

@danholdaway Added a figure showing comparisons between UFO final error without GSI duplicator versus GSI final observation error (see above). There are not much difference between them in sonde observations, but there are a lot differences in ship observations.

@asewnath
Copy link
Collaborator

Screen Shot 2023-02-15 at 3 34 04 PM

With duplicate check filter added, some duplicates not accounted for
Screen Shot 2023-02-15 at 3 34 30 PM

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Feb 15, 2023

Duplicate factor is applied.

corscat_eferr_vs_finaler
corscat_eferr_vs_finaler

@danholdaway danholdaway added waiting on jedi merge ufo done Observation operator evaluation in UFO completer labels Feb 15, 2023
@gmao-wgu
Copy link
Contributor

gmao-wgu commented Feb 22, 2023

The consistency check for Radiosonde wind (kx=220, 221) is completed, which can be demonstrated from the following plots:

The number of observations passing QC in the GSI:

sondes_density_qced_omb_vs_gsiomb_gsi_stats

The number of observations passing QC in Jedi:

sondes_density_qced_omb_vs_gsiomb_jedi_stats

The final observation error comparison:

sondes_corscat_qced_eferr_vs_gsifinal

@gmao-msienkie
Copy link
Contributor

I made a list of the ADPUPA data in the gdas.prepbufr file from 20201215/00z and found there were quite a few (87!) duplicated soundings on that date. (From the station IDs they appear to be Chinese sondes.) Possibly some of the soundings were transmitted both in BUFR and text format. In any case, the GSI 'errormod' filter would treat the soundings separately while the UFO filter would filter the soundings together if the reported release time were the same. While the DHR values given for the sounding pairs were different, the drift time for the first level of the soundings was the same. Actually the drift lat,lon,time was the same for quite a few levels (if not most) and so would trigger down-weighting by the duplicate check.

I was looking at this issue again because it came up when I was checking how the wind obs were handled by ObsErrorFactorConventional versus the GSIAdjustObsError. It seems to me that the agreement would be improved if we changed the definition of 'time_launch' in the GSI for raobs so that it uses the DHR (like the aircraft) rather than the first time level of the drfdat. It looks like the pairs of soundings have different DHR values (for whatever reason) while they often have the same values for the initial drfdat time. That puts the two reports in the same sounding for ObsErrorFactorConventional while the GSI processes the sounding reports in 'errormod' separately. It would be better to process those together (or actually better to screen for duplicates better in the preprocessing) but for the purpose of a match with GSI I think the change in the definition would help.

@gmao-jjin3
Copy link
Contributor

Thank you @gmao-msienkie . I made the change in read_prepbufr.f90. Let's see if we can get better results tomorrow.

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Apr 12, 2023

corscat_eferr_vs_finaler

The final observation errors do match after duplication is considered only for the same observation type such as sondes, surface, etc. in GSI. See the previous comparison when duplicate factor is considered for surface pressure in all conventional data: #78 (comment)

@gmao-jjin3
Copy link
Contributor

corscat_eferr_vs_finaler

@gmao-msienkie 's suggestion about launch_time does make a difference, removed one of the out of line two spots in the last figure posted on Feb 15. #78 (comment)

@gmao-jjin3
Copy link
Contributor

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Aug 16, 2023

density_qc_omb_vs_gsiomb
corscat_qc_omb_vs_hofx-GsiHofXBc
corscat_eferr_x_GsiDup_vs_finaler
corscat_eferr_vs_finalerr
Sonde surface pressures from x0048. Final observational error cannot be reproduced by UFO.

@gmao-jjin3
Copy link
Contributor

Just add a note that surface pressure hofx values in UFO match those in GSI after a bug fix which is merged in UFO developed branch: https://github.com/JCSDA-internal/ufo/pull/2987
corscat_qc_omb_vs_hofx-GsiHofXBc

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Aug 30, 2023

Validation against GEOS x0048 results (/discover/nobackup/jjin3/jedi/For_others/x0048.jj_20230818/conv/combine_conv/obs/sondes_obs_2021121200.nc4)

density_qc_omb_vs_gsiomb
corscat_qc_omb_vs_hofx-GsiHofXBc
corscat_EffectiveQC_vs_GsiEffectiveQC
corscat_eferr_x_GsiDup_vs_finaler
corscat_eferr_vs_finalerr

@gmao-jjin3
Copy link
Contributor

Validation against GEOS x0048 results (/discover/nobackup/jjin3/jedi/For_others/x0048.jj_20230818/conv/combine_conv/obs/sondes_obs_2021121200.nc4)
density_qc_omb_vs_gsiomb
corscat_qc_omb_vs_hofx-GsiHofXBc
corscat_EffectiveQC_vs_GsiEffectiveQC
corscat_eferr_x_GsiDup_vs_finaler
corscat_eferr_vs_finalerr

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
ufo done Observation operator evaluation in UFO completer waiting on jedi merge
Projects
None yet
Development

No branches or pull requests

6 participants