From 78b9e8895faaaebedae20db41fd9e066ca54be1a Mon Sep 17 00:00:00 2001 From: emilyhcliu <36091766+emilyhcliu@users.noreply.github.com> Date: Thu, 21 Mar 2024 14:41:31 -0400 Subject: [PATCH 1/3] Add updates for ozone and end-to-end results (#827) GDAS assimilates two types of ozone: - Total Column Ozone - Profile Ozone For the focus cycle: 2021080100Z, ozone data includes: - ompsnp npp (profile ozone) - ompstc npp (total ozone) - omi aura (total ozone) This PR includes the following updates for the three ozone data above: - BUFR converters - ObsSpace YAML files - Validation result Here is the updated QC flowchart for ozone: ![GSI-ozone-flowchart](https://github.com/NOAA-EMC/GDASApp/assets/36091766/102ee04a-30f8-4198-b211-023549f2cac8) --------- Co-authored-by: RussTreadon-NOAA Co-authored-by: Cory Martin Co-authored-by: BrettHoover-NOAA <98188219+BrettHoover-NOAA@users.noreply.github.com> Co-authored-by: Brett Hoover Co-authored-by: Xuanli Li <101414760+XuanliLi-NOAA@users.noreply.github.com> Co-authored-by: RussTreadon-NOAA <26926959+RussTreadon-NOAA@users.noreply.github.com> --- parm/atm/obs/config/omi_aura.yaml.j2 | 104 +++- parm/atm/obs/config/ompsnp_npp.yaml.j2 | 77 +-- ...ompstc8_npp.yaml.j2 => ompstc_npp.yaml.j2} | 58 ++- parm/atm/obs/lists/gdas_prototype_3d.yaml.j2 | 3 + parm/atm/obs/testing/omi_aura.yaml | 4 +- parm/atm/obs/testing/omi_aura_noqc.yaml | 5 +- parm/atm/obs/testing/ompsnp_npp.yaml | 4 +- parm/atm/obs/testing/ompsnp_npp_noqc.yaml | 5 +- .../{ompstc8_npp.yaml => ompstc_npp.yaml} | 16 +- ...tc8_npp_noqc.yaml => ompstc_npp_noqc.yaml} | 15 +- parm/ioda/bufr2ioda/bufr2ioda_ozone_omi.json | 20 + .../bufr2ioda/bufr2ioda_ozone_ompsnp.json | 21 + .../bufr2ioda/bufr2ioda_ozone_ompstc.json | 21 + ush/ioda/bufr2ioda/bufr2ioda_ozone_omi.py | 337 +++++++++++++ ush/ioda/bufr2ioda/bufr2ioda_ozone_ompsnp.py | 459 ++++++++++++++++++ ush/ioda/bufr2ioda/bufr2ioda_ozone_ompstc.py | 328 +++++++++++++ 16 files changed, 1388 insertions(+), 89 deletions(-) rename parm/atm/obs/config/{ompstc8_npp.yaml.j2 => ompstc_npp.yaml.j2} (61%) rename parm/atm/obs/testing/{ompstc8_npp.yaml => ompstc_npp.yaml} (79%) rename parm/atm/obs/testing/{ompstc8_npp_noqc.yaml => ompstc_npp_noqc.yaml} (66%) create mode 100644 parm/ioda/bufr2ioda/bufr2ioda_ozone_omi.json create mode 100644 parm/ioda/bufr2ioda/bufr2ioda_ozone_ompsnp.json create mode 100644 parm/ioda/bufr2ioda/bufr2ioda_ozone_ompstc.json create mode 100755 ush/ioda/bufr2ioda/bufr2ioda_ozone_omi.py create mode 100755 ush/ioda/bufr2ioda/bufr2ioda_ozone_ompsnp.py create mode 100755 ush/ioda/bufr2ioda/bufr2ioda_ozone_ompstc.py diff --git a/parm/atm/obs/config/omi_aura.yaml.j2 b/parm/atm/obs/config/omi_aura.yaml.j2 index 9cfd90dc9..427b34f05 100644 --- a/parm/atm/obs/config/omi_aura.yaml.j2 +++ b/parm/atm/obs/config/omi_aura.yaml.j2 @@ -3,7 +3,7 @@ obsdatain: engine: type: H5File - obsfile: '{{ DATA }}/obs/{{ OPREFIX }}omi_aura.{{ current_cycle | to_YMDH }}.nc' + obsfile: '{{ DATA }}/obs/{{ OPREFIX }}omi_aura.tm00.nc' obsdataout: engine: type: H5File @@ -11,13 +11,23 @@ io pool: max pool size: 1 simulated variables: [ozoneTotal] - + + #obs operator: + # name: AtmVertInterpLay + # geovals: [mole_fraction_of_ozone_in_air] + # coefficients: [0.007886131] # convert from ppmv to DU + # nlevels: [1] + obs operator: - name: AtmVertInterpLay - geovals: [mole_fraction_of_ozone_in_air] - coefficients: [0.007886131] # convert from ppmv to DU - nlevels: [1] - + name: ColumnRetrieval + nlayers_retrieval: 1 + tracer variables: [mole_fraction_of_ozone_in_air] + isApriori: false + isAveragingKernel: false + totalNoVertice: true + stretchVertices: topbottom #options: top, bottom, topbottom, none + model units coeff: 2.241398632746E-3 + obs pre filters: - filter: Perform Action filter variables: @@ -25,7 +35,33 @@ action: name: assign error error parameter: 6.0 - + + - filter: Create Diagnostic Flags + filter variables: + - name: ozoneTotal + flags: + - name: ObsValueSanityCheck + initial value: false + force reinitialization: false + - name: RowAnomaly + initial value: false + force reinitialization: false + - name: BadScan + initial value: false + force reinitialization: false + - name: Thinning + initial value: false + force reinitialization: false + - name: RetrievalQualityCodeFlag + initial value: false + force reinitialization: false + - name: RetrievalQualityAlgorithmFlag + initial value: false + force reinitialization: false + - name: GrossCheck + initial value: false + force reinitialization: false + obs prior filters: # GSI read routine QC # range sanity check @@ -34,9 +70,11 @@ - name: ozoneTotal minvalue: 0 maxvalue: 10000 - action: - name: reject - + actions: + - name: set + flag: ObsValueSanityCheck + - name: reject + # Do not use the data if row anomaly (bit 10)is 1 - filter: RejectList filter variables: @@ -45,7 +83,11 @@ - variable: name: MetaData/totalOzoneQualityFlag any_bit_set_of: 9 - + actions: + - name: set + flag: RowAnomaly + - name: reject + # Scan position check: reject scan position >= 25 - filter: RejectList filter variables: @@ -54,7 +96,11 @@ - variable: name: MetaData/sensorScanPosition minvalue: 25 - + actions: + - name: set + flag: BadScan + - name: reject + # Accept total_ozone_error_flag values of 0 and 1, but not any others. - filter: RejectList filter variables: @@ -63,7 +109,11 @@ - variable: name: MetaData/totalOzoneQualityCode is_not_in: 0, 1 - + actions: + - name: set + flag: RetrievalQualityCodeFlag + - name: reject + # Use data with best ozone algorighm - filter: RejectList filter variables: @@ -72,7 +122,21 @@ - variable: name: MetaData/bestOzoneAlgorithmFlag is_in: 3, 13 - + actions: + - name: set + flag: RetrievalQualityAlgorithmFlag + - name: reject + + # Data Thinning + - filter: Gaussian Thinning + horizontal_mesh: 150 + use_reduced_horizontal_grid: true + distance_norm: geodesic + actions: + - name: set + flag: Thinning + - name: reject + obs post filters: # GSI setup routine QC # Gross check @@ -81,9 +145,9 @@ - name: ozoneTotal threshold: 10.0 absolute threshold: 300.0 - action: - name: reject - + actions: + - name: set + flag: GrossCheck + - name: reject + # End of Filters - - diff --git a/parm/atm/obs/config/ompsnp_npp.yaml.j2 b/parm/atm/obs/config/ompsnp_npp.yaml.j2 index 82108dcd4..c95cc0274 100644 --- a/parm/atm/obs/config/ompsnp_npp.yaml.j2 +++ b/parm/atm/obs/config/ompsnp_npp.yaml.j2 @@ -3,7 +3,7 @@ obsdatain: engine: type: H5File - obsfile: '{{ DATA }}/obs/{{ OPREFIX }}ompsnp_npp.{{ current_cycle | to_YMDH }}.nc' + obsfile: '{{ DATA }}/obs/{{ OPREFIX }}ompsnp_npp.tm00.nc' obsgrouping: group variables: ["latitude"] sort variable: "pressure" @@ -15,13 +15,29 @@ io pool: max pool size: 1 simulated variables: [ozoneLayer] - + + #obs operator: + # name: AtmVertInterpLay + # geovals: [mole_fraction_of_ozone_in_air] + # coefficients: [0.007886131] # convert from ppmv to DU + # nlevels: [22] + obs operator: - name: AtmVertInterpLay - geovals: [mole_fraction_of_ozone_in_air] - coefficients: [0.007886131] # convert from ppmv to DU - nlevels: [22] - + name: ColumnRetrieval + nlayers_retrieval: 1 + tracer variables: [mole_fraction_of_ozone_in_air] + isApriori: false + isAveragingKernel: false + totalNoVertice: false + stretchVertices: none #options: top, bottom, topbottom, none + # model units coeff: 2.240013904035E-3 # this number to match the gsihofx values + model units coeff: 2.241398632746E-3 # this number to match the gsihofx values (use this) + # the actual scientific conversion factor is + # 2.1415E-3 kg[O3]/m-2 to DU + # so the name of the geovals + # is also likely wrong, as it apprears to be a mass + # fraction given the conversion factor needed + obs pre filters: # Observation error assignment - filter: Perform Action @@ -35,8 +51,9 @@ xvar: name: MetaData/pressure xvals: [0.001, 10.1325, 16.00935, 25.43258, 40.32735, 63.93607, 101.325, 160.0935, 254.3257, 403.2735, 639.3608, 1013.25, 1600.935, 2543.258, 4032.735, 6393.607, 10132.5, 16009.35, 25432.57, 40327.35, 63936.07, 101325] - errors: [7.7236, 0.02, 0.02, 0.025, 0.08, 0.15, 0.056, 0.125, 0.2, 0.299, 0.587, 0.864, 1.547, 2.718, 3.893, 4.353, 3.971, 4.407, 4.428, 3.312, 2.198, 2.285] - + errors: [7.7236, 0.020, 0.020, 0.025, 0.080, 0.150, 0.056, 0.125, 0.200, 0.299, 0.587, 0.864, 1.547, 2.718, 3.893, 4.353, 3.971, 4.407, 4.428, 3.312, 2.198, 2.285] + # errors: [7.7236, 0.020, 0.020, 0.025, 0.040, 0.080, 0.156, 0.245, 0.510, 1.098, 3.917, 6.124, 6.347, 5.798, 6.843, 9.253,10.091,10.967, 8.478, 5.572, 2.638, 3.525] # operational from gfs.v16.3.9 (late 2023) + obs prior filters: # Do not assimilation where pressure is zero # Zero pressure indicates the data is total column ozone @@ -47,7 +64,7 @@ - variable: name: MetaData/pressure maxvalue: 0.0001 - + # Sanity check on observaton values - filter: Bounds Check filter variables: @@ -56,7 +73,7 @@ maxvalue: 1000 action: name: reject - + # Total Ozone Quality Check (keeps 0, 2) # 0 indentifies good data # 2 identifies good data with a solar zenith angle > 84 degrees @@ -67,7 +84,7 @@ - variable: name: MetaData/totalOzoneQuality is_not_in: 0, 2 - + # Profile Ozone Quality Check (keeps 0, 1, 7) # 0 : good data # 1 : good data with a solar zenith angle > 84 degrees @@ -79,7 +96,7 @@ - variable: name: MetaData/profileOzoneQuality is_not_in: 0, 1, 7 - + obs post filters: # Gross error check - filter: Background Check @@ -92,7 +109,7 @@ - variable: name: MetaData/pressure maxvalue: 0.001 - + - filter: Background Check filter variables: - name: ozoneLayer @@ -104,7 +121,7 @@ name: MetaData/pressure minvalue: 30000.0 maxvalue: 110000.0 - + - filter: Background Check filter variables: - name: ozoneLayer @@ -116,7 +133,7 @@ name: MetaData/pressure minvalue: 20000.0 maxvalue: 30000.0 - + - filter: Background Check filter variables: - name: ozoneLayer @@ -128,7 +145,7 @@ name: MetaData/pressure minvalue: 10100.0 maxvalue: 20000.0 - + - filter: Background Check filter variables: - name: ozoneLayer @@ -140,7 +157,7 @@ name: MetaData/pressure minvalue: 6400.0 maxvalue: 10100.0 - + - filter: Background Check filter variables: - name: ozoneLayer @@ -152,7 +169,7 @@ name: MetaData/pressure minvalue: 4000.0 maxvalue: 6400.0 - + - filter: Background Check filter variables: - name: ozoneLayer @@ -164,7 +181,7 @@ name: MetaData/pressure minvalue: 2600.0 maxvalue: 4000.0 - + - filter: Background Check filter variables: - name: ozoneLayer @@ -176,7 +193,7 @@ name: MetaData/pressure minvalue: 1600.0 maxvalue: 2600.0 - + - filter: Background Check filter variables: - name: ozoneLayer @@ -188,7 +205,7 @@ name: MetaData/pressure minvalue: 1100.0 maxvalue: 1600.0 - + - filter: Background Check filter variables: - name: ozoneLayer @@ -200,7 +217,7 @@ name: MetaData/pressure minvalue: 700.0 maxvalue: 1100.0 - + - filter: Background Check filter variables: - name: ozoneLayer @@ -212,7 +229,7 @@ name: MetaData/pressure minvalue: 400.0 maxvalue: 700.0 - + - filter: Background Check filter variables: - name: ozoneLayer @@ -224,7 +241,7 @@ name: MetaData/pressure minvalue: 300.0 maxvalue: 400.0 - + - filter: Background Check filter variables: - name: ozoneLayer @@ -236,7 +253,7 @@ name: MetaData/pressure minvalue: 200.0 maxvalue: 300.0 - + - filter: Background Check filter variables: - name: ozoneLayer @@ -248,7 +265,7 @@ name: MetaData/pressure maxvalue: 70.0 maxvalue: 200.0 - + - filter: Background Check filter variables: - name: ozoneLayer @@ -260,7 +277,7 @@ name: MetaData/pressure minvalue: 40.0 maxvalue: 70.0 - + - filter: Background Check filter variables: - name: ozoneLayer @@ -272,7 +289,7 @@ name: MetaData/pressure minvalue: 30.0 maxvalue: 40.0 - + - filter: Background Check filter variables: - name: ozoneLayer @@ -283,5 +300,5 @@ - variable: name: MetaData/pressure maxvalue: 30.0 - + # End of Filters diff --git a/parm/atm/obs/config/ompstc8_npp.yaml.j2 b/parm/atm/obs/config/ompstc_npp.yaml.j2 similarity index 61% rename from parm/atm/obs/config/ompstc8_npp.yaml.j2 rename to parm/atm/obs/config/ompstc_npp.yaml.j2 index e22eee7c8..8fd6b69cd 100644 --- a/parm/atm/obs/config/ompstc8_npp.yaml.j2 +++ b/parm/atm/obs/config/ompstc_npp.yaml.j2 @@ -1,23 +1,33 @@ - obs space: - name: ompstc8_npp + name: ompstc_npp obsdatain: engine: type: H5File - obsfile: '{{ DATA }}/obs/{{ OPREFIX }}ompstc8_npp.{{ current_cycle | to_YMDH }}.nc' + obsfile: '{{ DATA }}/obs/{{ OPREFIX }}ompstc_npp.tm00.nc' obsdataout: engine: type: H5File - obsfile: '{{ DATA }}/diags/diag_ompstc8_npp_{{ current_cycle | to_YMDH }}.nc' + obsfile: '{{ DATA }}/diags/diag_ompstc_npp_{{ current_cycle | to_YMDH }}.nc' io pool: max pool size: 1 simulated variables: [ozoneTotal] - + + #obs operator: + # name: AtmVertInterpLay + # geovals: [mole_fraction_of_ozone_in_air] + # coefficients: [0.007886131] # convert from ppmv to DU + # nlevels: [1] + obs operator: - name: AtmVertInterpLay - geovals: [mole_fraction_of_ozone_in_air] - coefficients: [0.007886131] # convert from ppmv to DU - nlevels: [1] - + name: ColumnRetrieval + nlayers_retrieval: 1 + tracer variables: [mole_fraction_of_ozone_in_air] + isApriori: false + isAveragingKernel: false + totalNoVertice: true + stretchVertices: topbottom # options: top, bottom, topbottom, none + model units coeff: 2.241398632746E-3 + obs pre filters: - filter: Perform Action filter variables: @@ -25,7 +35,7 @@ action: name: assign error error parameter: 6.0 - + obs prior filters: # GSI read routine QC # range sanity check @@ -36,16 +46,23 @@ maxvalue: 1000 action: name: reject - + + #- filter: Gaussian Thinning + # horizontal_mesh: 150 + # use_reduced_horizontal_grid: true + # distance_norm: geodesic + # action: + # name: reject + # Accept total_ozone_error_flag values of 0 and 1, but not any others. - filter: RejectList filter variables: - name: ozoneTotal where: - variable: - name: MetaData/totalOzoneQualityCode + name: MetaData/totalOzoneQualityCode is_not_in: 0, 1 - + - filter: RejectList filter variables: - name: ozoneTotal @@ -53,7 +70,7 @@ - variable: name: MetaData/bestOzoneAlgorithmFlag is_in: 3, 13 - + # GSI setup routine QC - filter: RejectList filter variables: @@ -65,7 +82,7 @@ - variable: name: MetaData/latitude minvalue: 50.0 - + - filter: RejectList filter variables: - name: ozoneTotal @@ -76,7 +93,14 @@ - variable: name: MetaData/latitude maxvalue: -50.0 - + + - filter: Gaussian Thinning + horizontal_mesh: 150 + use_reduced_horizontal_grid: true + distance_norm: geodesic + action: + name: reject + obs post filters: - filter: Background Check filter variables: @@ -85,5 +109,5 @@ absolute threshold: 300.0 action: name: reject - + # End of Filters diff --git a/parm/atm/obs/lists/gdas_prototype_3d.yaml.j2 b/parm/atm/obs/lists/gdas_prototype_3d.yaml.j2 index a8bc25019..f576fc25e 100644 --- a/parm/atm/obs/lists/gdas_prototype_3d.yaml.j2 +++ b/parm/atm/obs/lists/gdas_prototype_3d.yaml.j2 @@ -4,4 +4,7 @@ observers: {% include 'atm/obs/config/scatwind_ascat_metop-b.yaml.j2' %} {% include 'atm/obs/config/conv_ps.yaml.j2' %} {% include 'atm/obs/config/gnssro.yaml.j2' %} +{% include 'atm/obs/config/ompsnp_npp.yaml.j2' %} +{% include 'atm/obs/config/ompstc_npp.yaml.j2' %} +{% include 'atm/obs/config/omi_aura.yaml.j2' %} {% endfilter %} diff --git a/parm/atm/obs/testing/omi_aura.yaml b/parm/atm/obs/testing/omi_aura.yaml index 8c5bb52da..b409700ad 100644 --- a/parm/atm/obs/testing/omi_aura.yaml +++ b/parm/atm/obs/testing/omi_aura.yaml @@ -88,4 +88,6 @@ obs post filters: name: reject # End of Filters -passedBenchmark: 1170 +passedBenchmark: 1182 # GSI modified : keep bad data and then thinned ----- this is to check QC in the read routine + # GSI original: remove bad data and then thinned + # original: thinned: 2402; unthinned: 53768 diff --git a/parm/atm/obs/testing/omi_aura_noqc.yaml b/parm/atm/obs/testing/omi_aura_noqc.yaml index d32954a8f..0cabf0c9a 100644 --- a/parm/atm/obs/testing/omi_aura_noqc.yaml +++ b/parm/atm/obs/testing/omi_aura_noqc.yaml @@ -11,12 +11,15 @@ obs space: io pool: max pool size: 1 simulated variables: [ozoneTotal] + geovals: filename: !ENV omi_aura_geoval_${CDATE}.nc4 + obs operator: name: AtmVertInterpLay geovals: [mole_fraction_of_ozone_in_air] coefficients: [0.007886131] # convert from ppmv to DU +# coefficients: [0.0078976797] # convert from ppmv to DU nlevels: [1] obs pre filters: @@ -27,6 +30,6 @@ obs pre filters: name: assign error error parameter: 6.0 -passedBenchmark: 4927 # total:6082; missing:1155 +passedBenchmark: 4980 # total:6140; missing:1180 #vector ref: GsiHofXBc #tolerance: 1.e-5 diff --git a/parm/atm/obs/testing/ompsnp_npp.yaml b/parm/atm/obs/testing/ompsnp_npp.yaml index b7974f5cd..0071d3284 100644 --- a/parm/atm/obs/testing/ompsnp_npp.yaml +++ b/parm/atm/obs/testing/ompsnp_npp.yaml @@ -38,8 +38,8 @@ obs pre filters: xvar: name: MetaData/pressure xvals: [0.001, 10.1325, 16.00935, 25.43258, 40.32735, 63.93607, 101.325, 160.0935, 254.3257, 403.2735, 639.3608, 1013.25, 1600.935, 2543.258, 4032.735, 6393.607, 10132.5, 16009.35, 25432.57, 40327.35, 63936.07, 101325] - errors: [7.7236, 0.02, 0.02, 0.025, 0.08, 0.15, 0.056, 0.125, 0.2, 0.299, 0.587, 0.864, 1.547, 2.718, 3.893, 4.353, 3.971, 4.407, 4.428, 3.312, 2.198, 2.285] - + errors: [7.7236, 0.020, 0.020, 0.025, 0.080, 0.150, 0.056, 0.125, 0.200, 0.299, 0.587, 0.864, 1.547, 2.718, 3.893, 4.353, 3.971, 4.407, 4.428, 3.312, 2.198, 2.285] +# errors: [7.7236, 0.020, 0.020, 0.025, 0.040, 0.080, 0.156, 0.245, 0.510, 1.098, 3.917, 6.124, 6.347, 5.798, 6.843, 9.253, 10.091, 10.967, 8.478, 5.572, 2.638, 3.525] # operational from gfs.v16.3.9 (late 2023) obs prior filters: # Do not assimilation where pressure is zero # Zero pressure indicates the data is total column ozone diff --git a/parm/atm/obs/testing/ompsnp_npp_noqc.yaml b/parm/atm/obs/testing/ompsnp_npp_noqc.yaml index 639a8dda5..9a2d6c6b4 100644 --- a/parm/atm/obs/testing/ompsnp_npp_noqc.yaml +++ b/parm/atm/obs/testing/ompsnp_npp_noqc.yaml @@ -15,8 +15,10 @@ obs space: io pool: max pool size: 1 simulated variables: [ozoneLayer] + geovals: filename: !ENV ompsnp_npp_geoval_${CDATE}.nc4 + obs operator: name: AtmVertInterpLay geovals: [mole_fraction_of_ozone_in_air] @@ -36,7 +38,8 @@ obs pre filters: xvar: name: MetaData/pressure xvals: [0.001, 10.1325, 16.00935, 25.43258, 40.32735, 63.93607, 101.325, 160.0935, 254.3257, 403.2735, 639.3608, 1013.25, 1600.935, 2543.258, 4032.735, 6393.607, 10132.5, 16009.35, 25432.57, 40327.35, 63936.07, 101325] - errors: [7.7236, 0.02, 0.02, 0.025, 0.08, 0.15, 0.056, 0.125, 0.2, 0.299, 0.587, 0.864, 1.547, 2.718, 3.893, 4.353, 3.971, 4.407, 4.428, 3.312, 2.198, 2.285] + errors: [7.7236, 0.020, 0.020, 0.025, 0.080, 0.150, 0.056, 0.125, 0.200, 0.299, 0.587, 0.864, 1.547, 2.718, 3.893, 4.353, 3.971, 4.407, 4.428, 3.312, 2.198, 2.285] +# errors: [7.7236, 0.020, 0.020, 0.025, 0.040, 0.080, 0.156, 0.245, 0.510, 1.098, 3.917, 6.124, 6.347, 5.798, 6.843, 9.253, 10.091, 10.967, 8.478, 5.572, 2.638, 3.525] # operational from gfs.v16.3.9 (late 2023) passedBenchmark: 6314 # total:6314; missing:0 #vector ref: GsiHofXBc diff --git a/parm/atm/obs/testing/ompstc8_npp.yaml b/parm/atm/obs/testing/ompstc_npp.yaml similarity index 79% rename from parm/atm/obs/testing/ompstc8_npp.yaml rename to parm/atm/obs/testing/ompstc_npp.yaml index 8fb3d6165..8f10bb5d2 100644 --- a/parm/atm/obs/testing/ompstc8_npp.yaml +++ b/parm/atm/obs/testing/ompstc_npp.yaml @@ -1,23 +1,19 @@ obs space: - name: ompstc8_npp + name: ompstc_npp obsdatain: engine: type: H5File - obsfile: !ENV ompstc8_npp_obs_${CDATE}.nc4 - obsgrouping: - group variables: ["latitude"] - sort variable: "pressure" - sort order: "ascending" + obsfile: !ENV ompstc_npp_obs_${CDATE}.nc4 obsdataout: engine: type: H5File - obsfile: !ENV ompstc8_npp_diag_${CDATE}.nc4 + obsfile: !ENV ompstc_npp_diag_${CDATE}.nc4 io pool: max pool size: 1 simulated variables: [ozoneTotal] geovals: - filename: !ENV ompstc8_npp_geoval_${CDATE}.nc4 + filename: !ENV ompstc_npp_geoval_${CDATE}.nc4 obs operator: name: AtmVertInterpLay @@ -94,4 +90,6 @@ obs post filters: name: reject # End of Filters -passedBenchmark: 6130 +passedBenchmark: 6130 # GSI modified : keep bad data and then thinned ----- this is to check QC in the read routine + # GSI original: remove bad data and then thinned + # original: thinned: 6214; unthinned: 44899 diff --git a/parm/atm/obs/testing/ompstc8_npp_noqc.yaml b/parm/atm/obs/testing/ompstc_npp_noqc.yaml similarity index 66% rename from parm/atm/obs/testing/ompstc8_npp_noqc.yaml rename to parm/atm/obs/testing/ompstc_npp_noqc.yaml index de9ef7839..98fa0062f 100644 --- a/parm/atm/obs/testing/ompstc8_npp_noqc.yaml +++ b/parm/atm/obs/testing/ompstc_npp_noqc.yaml @@ -1,27 +1,26 @@ obs space: - name: ompstc8_npp + name: ompstc_npp obsdatain: engine: type: H5File - obsfile: !ENV ompstc8_npp_obs_${CDATE}.nc4 - obsgrouping: - group variables: ["latitude"] - sort variable: "pressure" - sort order: "ascending" + obsfile: !ENV ompstc_npp_obs_${CDATE}.nc4 obsdataout: engine: type: H5File - obsfile: !ENV ompstc8_npp_diag_${CDATE}.nc4 + obsfile: !ENV ompstc_npp_diag_${CDATE}.nc4 io pool: max pool size: 1 simulated variables: [ozoneTotal] + geovals: - filename: !ENV ompstc8_npp_geoval_${CDATE}.nc4 + filename: !ENV ompstc_npp_geoval_${CDATE}.nc4 + obs operator: name: AtmVertInterpLay geovals: [mole_fraction_of_ozone_in_air] coefficients: [0.007886131] # convert from ppmv to DU nlevels: [1] + obs pre filters: - filter: Perform Action filter variables: diff --git a/parm/ioda/bufr2ioda/bufr2ioda_ozone_omi.json b/parm/ioda/bufr2ioda/bufr2ioda_ozone_omi.json new file mode 100644 index 000000000..ff851fba1 --- /dev/null +++ b/parm/ioda/bufr2ioda/bufr2ioda_ozone_omi.json @@ -0,0 +1,20 @@ +{ + "data_format" : "bufr_d", + "data_type" : "omi", + "cycle_type" : "{{ RUN }}", + "cycle_datetime" : "{{ current_cycle | to_YMDH }}", + "dump_directory" : "{{ DMPDIR }}", + "ioda_directory" : "{{ COM_OBS }}", + "ioda_type" : "omi", + "subsets" : "NC008013", + "source" : "MTYP 008-013 OZONE MONITORING INSTRUMENT (OMI) Total Column Ozone", + "process_level" : "Level-2", + "platform_description" : "NASA EOS program for atmospheric chemistry monitoriny - 3rd flight unit since 2004", + "sensor_description" : "UV/VIS grating imaging spectrometer, three bands, 1560 channels total", + "data_description" : "Level-2 retrieved total zone from Ozone Monitoring Instrument (OMI)", + "data_provider" : "U.S. NOAA/NESDIS", + "sensor_info" : { "sensor_name": "OMI", "sensor_full_name": "Ozone Monitoring Instrument", "sensor_id": 394 }, + "satellite_info" : [ + { "satellite_name": "AURA", "satellite_full_name": "Earth Observation System - Aura", "satellite_id": 785, "launch time": "20040715" } + ] +} diff --git a/parm/ioda/bufr2ioda/bufr2ioda_ozone_ompsnp.json b/parm/ioda/bufr2ioda/bufr2ioda_ozone_ompsnp.json new file mode 100644 index 000000000..e177548f8 --- /dev/null +++ b/parm/ioda/bufr2ioda/bufr2ioda_ozone_ompsnp.json @@ -0,0 +1,21 @@ +{ + "data_format" : "bufr_d", + "data_type" : "ompsn8", + "cycle_type" : "{{ RUN }}", + "cycle_datetime" : "{{ current_cycle | to_YMDH }}", + "dump_directory" : "{{ DMPDIR }}", + "ioda_directory" : "{{ COM_OBS }}", + "ioda_type" : "ompsnp", + "subsets" : "NC008017", + "source" : "MTYP 008-017 OMPS NADIR PROFILE OZONE (VERSION 8 BUFR)", + "process_level" : "Level-2", + "platform_description" : "Joint Polar Satellite System operaring since 2011 from 833 km altitude in sunsychronous orbits", + "sensor_description" : "UV grating spectrometer of spectral resolution 1 nm (sampling at 0.42 nm); Nadir-viewing ozone profiler: 250-310 nm; one along-track sounding every 38 s at 250 km resolution", + "data_description" : "OMPS NP is the ozone profile derived from Nadir Profile Spectrometeor in a single ground pixel of 250x250 km at nadir from SNPP; 5 ground pixel of 50x50 km at nadir from JPSS", + "data_provider" : "U.S. NOAA/NESDIS", + "sensor_info" : { "sensor_name": "ompsnp", "sensor_full_name": "Ozone Mapping and Profiler Suite", "sensor_id": 947 }, + "satellite_info" : [ + { "satellite_name": "NPP", "satellite_full_name": "Suomi National Polar-orbiting Partnership", "satellite_id": 224, "launch time": "20111028" }, + { "satellite_name": "N20", "satellite_full_name": "National Oceanic and Atmospheric Administration - 20", "satellite_id": 225, "launch time": "20171118" } + ] +} diff --git a/parm/ioda/bufr2ioda/bufr2ioda_ozone_ompstc.json b/parm/ioda/bufr2ioda/bufr2ioda_ozone_ompstc.json new file mode 100644 index 000000000..a86404a27 --- /dev/null +++ b/parm/ioda/bufr2ioda/bufr2ioda_ozone_ompstc.json @@ -0,0 +1,21 @@ +{ + "data_format" : "bufr_d", + "data_type" : "ompst8", + "cycle_type" : "{{ RUN }}", + "cycle_datetime" : "{{ current_cycle | to_YMDH }}", + "dump_directory" : "{{ DMPDIR }}", + "ioda_directory" : "{{ COM_OBS }}", + "ioda_type" : "ompstc", + "subsets" : "NC008018", + "source" : "MTYP 008-018 OMPS TOTAL COLUMN OZONE (VERSION 8 BUFR) ", + "process_level" : "Level-2", + "platform_description" : "Joint Polar Satellite System operaring since 2011 from 833 km altitude in sunsychronous orbits", + "sensor_description" : "UV grating spectrometer of spectral resolution 1 nm (sampling at 0.42 nm); corss-track mapper for total ozone: 300-380 nm; corss-track swath 2800 km; along-track one 50-km strip in 7.6 s", + "data_description" : "OMPS TC is the Total Column Ozone derived from UV and near-infrared observations (between 300-380 nm) made by the OMPS Nadir Total Column Mapper S-NPP (50 km) and N-21 (10 km)", + "data_provider" : "U.S. NOAA/NESDIS", + "sensor_info" : { "sensor_name": "OMPS-nadir Mapper", "sensor_full_name": "Ozone Mapping and Profiler Suite", "sensor_id": 947 }, + "satellite_info" : [ + { "satellite_name": "NPP", "satellite_full_name": "Suomi National Polar-orbiting Partnership", "satellite_id": 224, "launch time": "20111028" }, + { "satellite_name": "N20", "satellite_full_name": "National Oceanic and Atmospheric Administration - 20", "satellite_id": 225, "launch time": "20171118" } + ] +} diff --git a/ush/ioda/bufr2ioda/bufr2ioda_ozone_omi.py b/ush/ioda/bufr2ioda/bufr2ioda_ozone_omi.py new file mode 100755 index 000000000..fba51a9ce --- /dev/null +++ b/ush/ioda/bufr2ioda/bufr2ioda_ozone_omi.py @@ -0,0 +1,337 @@ +#!/usr/bin/env python3 +import argparse +import numpy as np +import numpy.ma as ma +from pyiodaconv import bufr +import calendar +import json +import time +import math +import datetime +import os +from datetime import datetime +from pyioda import ioda_obs_space as ioda_ospace +from wxflow import Logger + +# ======================================================================= +# Subset | Description | +# ----------------------------------------------------------------------- +# NC008013 | Level-2 total column ozone retrievals from OMI on Earth | +# | Observation System - Aura | +# ======================================================================= + + +def bufr_to_ioda(config, logger): + + # ================================== + # Get parameters from configuration + # ================================== + subsets = config["subsets"] + source = config["source"] + data_format = config["data_format"] + data_type = config["data_type"] + data_description = config["data_description"] + data_provider = config["data_provider"] + cycle_type = config["cycle_type"] + dump_dir = config["dump_directory"] + ioda_dir = config["ioda_directory"] + ioda_type = config["ioda_type"] + cycle = config["cycle_datetime"] + yyyymmdd = cycle[0:8] + hh = cycle[8:10] + + satellite_info_array = config["satellite_info"] + sensor_name = config["sensor_info"]["sensor_name"] + sensor_full_name = config["sensor_info"]["sensor_full_name"] + sensor_id = config["sensor_info"]["sensor_id"] + + # Get derived parameters + yyyymmdd = cycle[0:8] + hh = cycle[8:10] + reference_time = datetime.strptime(cycle, "%Y%m%d%H") + reference_time = reference_time.strftime("%Y-%m-%dT%H:%M:%SZ") + + # General informaton for global attribute + converter = 'BUFR Converter' + process_level = config["process_level"] + platform_description = config["platform_description"] + sensor_description = config["sensor_description"] + + logger.info(f'Processing {data_description}') + logger.info(f'reference_time = {reference_time}') + + bufrfile = f"{cycle_type}.t{hh}z.{data_type}.tm00.{data_format}" + DATA_PATH = os.path.join(dump_dir, f"{cycle_type}.{yyyymmdd}", str(hh), 'atmos', bufrfile) + if not os.path.isfile(DATA_PATH): + logger.info(f"DATA_PATH {DATA_PATH} does not exist") + return + logger.debug(f"The DATA_PATH is: {DATA_PATH}") + + # ============================================ + # Make the QuerySet for all the data we want + # ============================================ + start_time = time.time() + + logger.info('Making QuerySet') + q = bufr.QuerySet() + + # MetaData + q.add('latitude', '*/CLAT') + q.add('longitude', '*/CLON') + q.add('satelliteId', '*/SAID') + q.add('year', '*/YEAR') + q.add('month', '*/MNTH') + q.add('day', '*/DAYS') + q.add('hour', '*/HOUR') + q.add('minute', '*/MINU') + q.add('second', '*/SECO') + q.add('solarZenithAngle', '*/SOZA') + q.add('sensorScanPosition', '*/FOVN') + + # Quality + q.add('bestOzoneAlgorithmFlag', '*/AFBO') + q.add('totalOzoneQualityFlag', '*/TOQF') + q.add('totalOzoneQualityCode', '*/TOQC') + + # ObsValue + q.add('ozoneTotal', '*/OZON') + + end_time = time.time() + running_time = end_time - start_time + logger.debug(f'Processing time for making QuerySet : {running_time} seconds') + + # ============================================================== + # Open the BUFR file and execute the QuerySet to get ResultSet + # Use the ResultSet returned to get numpy arrays of the data + # ============================================================== + start_time = time.time() + + logger.info('Executing QuerySet to get ResultSet') + with bufr.File(DATA_PATH) as f: + try: + r = f.execute(q) + except Exception as err: + logger.info(f'Return with {err}') + return + + # MetaData + satid = r.get('satelliteId') + year = r.get('year') + month = r.get('month') + day = r.get('day') + hour = r.get('hour') + minute = r.get('minute') + second = r.get('second') + lat = r.get('latitude') + lon = r.get('longitude') + solzenang = r.get('solarZenithAngle') + scanpos = r.get('sensorScanPosition') + + # Quality Information + toqc = r.get('totalOzoneQualityCode') + toqf = r.get('totalOzoneQualityFlag') + afbo = r.get('bestOzoneAlgorithmFlag') + + # ObsValue + # Total Ozone + o3val = r.get('ozoneTotal', type='float') + + # DateTime: seconds since Epoch time + # IODA has no support for numpy datetime arrays dtype=datetime64[s] + timestamp = r.get_datetime('year', 'month', 'day', 'hour', 'minute', 'second').astype(np.int64) + + end_time = time.time() + running_time = end_time - start_time + logger.debug(f'Processing time for executing QuerySet to get ResultSet : {running_time} seconds') + + # ========================= + # Create derived variables + # ========================= + start_time = time.time() + + # Create pressure variable and fill in with zeros + pressure = np.ma.array(np.zeros(lat.shape[0], dtype=np.float32)) + pressure.mask = lat.mask + pressure.fill_value = lat.fill_value + + end_time = time.time() + running_time = end_time - start_time + print('Running time for creating derived variables : ', running_time, 'seconds') + + # ===================================== + # Split output based on satellite id + # Create IODA ObsSpace + # Write IODA output + # ===================================== + logger.info('Split data based on satellite id, Create IODA ObsSpace and Write IODA output') + + # Find unique satellite identifiers in data to process + unique_satids = np.unique(satid) + logger.info(f'Number of Unique satellite identifiers : {len(unique_satids)}') + logger.info(f'Unique satellite identifiers: {unique_satids}') + + logger.debug(f'Loop through unique satellite identifier : {unique_satids}') + total_ob_processed = 0 + for sat in unique_satids.tolist(): + start_time = time.time() + + matched = False + for satellite_info in satellite_info_array: + if (satellite_info["satellite_id"] == sat): + matched = True + satellite_id = satellite_info["satellite_id"] + satellite_name = satellite_info["satellite_name"] + satinst = sensor_name.lower()+'_'+satellite_name.lower() + logger.debug(f'Split data for {satinst} satid = {sat}') + + if matched: + + # Define a boolean mask to subset data from the original data object + mask = satid == sat + # MetaData + lon2 = lon[mask] + lat2 = lat[mask] + satid2 = satid[mask] + solzenang2 = solzenang[mask] + scanpos2 = scanpos[mask] + timestamp2 = timestamp[mask] + pressure2 = pressure[mask] + + # QC Info + toqc2 = toqc[mask] + toqf2 = toqf[mask] + afbo2 = afbo[mask] + + # ObsValue + o3val2 = o3val[mask] + + timestamp2_min = datetime.fromtimestamp(timestamp2.min()) + timestamp2_max = datetime.fromtimestamp(timestamp2.max()) + + # Create the dimensions + dims = { + 'Location': np.arange(0, lat2.shape[0]), + } + + # Create IODA ObsSpace + sat = satellite_name.lower() + iodafile = f"{cycle_type}.t{hh}z.{ioda_type}_{sat}.tm00.nc" + OUTPUT_PATH = os.path.join(ioda_dir, iodafile) + logger.info(f'Create output file : {OUTPUT_PATH}') + obsspace = ioda_ospace.ObsSpace(OUTPUT_PATH, mode='w', dim_dict=dims) + + # Create Global attributes + logger.debug('Create global attributes') + obsspace.write_attr('sourceFiles', bufrfile) + obsspace.write_attr('source', source) + obsspace.write_attr('description', data_description) + obsspace.write_attr('datetimeReference', reference_time) + obsspace.write_attr('Converter', converter) + obsspace.write_attr('platformLongDescription', platform_description) + obsspace.write_attr('platformCommonName', satellite_name) + obsspace.write_attr('platform', satellite_id) + obsspace.write_attr('sensorLongDescription', sensor_description) + obsspace.write_attr('sensorCommonName', sensor_name) + obsspace.write_attr('sensor', sensor_id) + obsspace.write_attr('dataProviderOrigin', data_provider) + obsspace.write_attr('processingLevel', process_level) + obsspace.write_attr('datetimeRange', [str(timestamp2_min), str(timestamp2_max)]) + + # Create IODA variables + logger.debug('Create variables: name, type, units, and attributes') + # Longitude + obsspace.create_var('MetaData/longitude', dtype=lon2.dtype, fillval=lon2.fill_value) \ + .write_attr('units', 'degrees_east') \ + .write_attr('valid_range', np.array([-180, 180], dtype=np.float32)) \ + .write_attr('long_name', 'Longitude') \ + .write_data(lon2) + + # Latitude + obsspace.create_var('MetaData/latitude', dtype=lat2.dtype, fillval=lat2.fill_value) \ + .write_attr('units', 'degrees_north') \ + .write_attr('valid_range', np.array([-90, 90], dtype=np.float32)) \ + .write_attr('long_name', 'Latitude') \ + .write_data(lat2) + + # Datetime + obsspace.create_var('MetaData/dateTime', dtype=timestamp2.dtype, fillval=timestamp2.fill_value) \ + .write_attr('units', 'seconds since 1970-01-01T00:00:00Z') \ + .write_attr('long_name', 'Unix Epoch') \ + .write_data(timestamp2) + + # Satellite Identifier + obsspace.create_var('MetaData/satelliteIdentifier', dtype=satid2.dtype, fillval=satid2.fill_value) \ + .write_attr('long_name', 'Satellite Identifier') \ + .write_data(satid2) + + # Total Ozone Quality Flag + obsspace.create_var('MetaData/totalOzoneQualityFlag', dtype=toqf2.dtype, fillval=toqf2.fill_value) \ + .write_attr('long_name', 'Total Ozone Quality Flag ') \ + .write_data(toqf2) + + # Total Ozone Quality Code + obsspace.create_var('MetaData/totalOzoneQualityCode', dtype=toqc2.dtype, fillval=toqc2.fill_value) \ + .write_attr('long_name', 'OMI Total Ozone Quality Code') \ + .write_data(toqc2) + + # Pressure + obsspace.create_var('MetaData/pressure', dtype=pressure2.dtype, fillval=pressure2.fill_value) \ + .write_attr('units', 'pa') \ + .write_attr('long_name', 'Pressure') \ + .write_data(pressure2) + + # Algorithm Flag for Best Ozone + obsspace.create_var('MetaData/bestOzoneAlgorithmFlag', dtype=afbo2.dtype, fillval=afbo2.fill_value) \ + .write_attr('long_name', 'Algorithm Flag for Best Ozone') \ + .write_data(afbo2) + + # Solar Zenith Angle + obsspace.create_var('MetaData/solarZenithAngle', dtype=solzenang2.dtype, fillval=solzenang2.fill_value) \ + .write_attr('units', 'm') \ + .write_attr('long_name', 'Solar Zenith Angle') \ + .write_data(solzenang2) + + # Sensor Scan Position + obsspace.create_var('MetaData/sensorScanPosition', dtype=scanpos2.dtype, fillval=scanpos2.fill_value) \ + .write_attr('long_name', 'Sensor Scan Position') \ + .write_data(scanpos2) + + # Total Ozone + obsspace.create_var('ObsValue/ozoneTotal', dtype=o3val2.dtype, fillval=o3val2.fill_value) \ + .write_attr('units', 'DU') \ + .write_attr('long_name', 'Total Column Ozone') \ + .write_data(o3val2) + + end_time = time.time() + running_time = end_time - start_time + total_ob_processed += len(satid2) + logger.debug(f'Processing time for splitting and output IODA for {satinst} : {running_time}, seconds') + + else: + logger.info(f'Do not find this satellite id in the configuration: satid = {sat}') + + logger.info('All Done!') + logger.info(f'Total number of observation processed : {total_ob_processed}') + + +if __name__ == '__main__': + + start_time = time.time() + + parser = argparse.ArgumentParser() + parser.add_argument('-c', '--config', type=str, help='Input JSON configuration', required=True) + parser.add_argument('-v', '--verbose', help='print debug logging information', + action='store_true') + args = parser.parse_args() + + log_level = 'DEBUG' if args.verbose else 'INFO' + logger = Logger('BUFR2IODA_ozone_omi.py', level=log_level, colored_log=True) + + with open(args.config, "r") as json_file: + config = json.load(json_file) + + bufr_to_ioda(config, logger) + + end_time = time.time() + running_time = end_time - start_time + logger.info(f"Total running time : {running_time} seconds") diff --git a/ush/ioda/bufr2ioda/bufr2ioda_ozone_ompsnp.py b/ush/ioda/bufr2ioda/bufr2ioda_ozone_ompsnp.py new file mode 100755 index 000000000..025d81216 --- /dev/null +++ b/ush/ioda/bufr2ioda/bufr2ioda_ozone_ompsnp.py @@ -0,0 +1,459 @@ +#!/usr/bin/env python3 +import argparse +import numpy as np +import numpy.ma as ma +from pyiodaconv import bufr +import calendar +import json +import time +import math +import datetime +import os +from datetime import datetime +from pyioda import ioda_obs_space as ioda_ospace +from wxflow import Logger + +# ================================================================================================== +# Subset | Description (OMPS NP) | +# -------------------------------------------------------------------------------------------------- +# NC008017 | OMPS NP is the Ozone Nadir Profile derived from the measurements made by the OMPS | +# | Nadir Profiler and the OMPS Nadir Mapper sensors on Suomi-NPP and JPSS. The data | +# | contains layer ozone amounts derived from the ratio of the observed backscattered | +# | spectral radiance to the incoming solar spectral irradiance for all solar zenith | +# | angle viewing conditions less than or equal to 80 degrees. | +# ================================================================================================== + + +def format_element(x): + + field_width = 15 + decimal_places = 5 + scientific_notation_threshold = 1000 + + if abs(x) < scientific_notation_threshold: + return f"{x:>{field_width}.{decimal_places}f}" + else: + return f"{x:>{field_width}.5e}" + + +def bufr_to_ioda(config, logger): + + # Get parameters from configuration + subsets = config["subsets"] + source = config["source"] + data_format = config["data_format"] + data_type = config["data_type"] + data_description = config["data_description"] + data_provider = config["data_provider"] + cycle_type = config["cycle_type"] + dump_dir = config["dump_directory"] + ioda_dir = config["ioda_directory"] + ioda_type = config["ioda_type"] + cycle = config["cycle_datetime"] + yyyymmdd = cycle[0:8] + hh = cycle[8:10] + + satellite_info_array = config["satellite_info"] + sensor_name = config["sensor_info"]["sensor_name"] + sensor_full_name = config["sensor_info"]["sensor_full_name"] + sensor_id = config["sensor_info"]["sensor_id"] + + # Get derived parameters + yyyymmdd = cycle[0:8] + hh = cycle[8:10] + reference_time = datetime.strptime(cycle, "%Y%m%d%H") + reference_time = reference_time.strftime("%Y-%m-%dT%H:%M:%SZ") + + # General informaton for global attribute + converter = 'BUFR to IODA Converter' + process_level = config["process_level"] + platform_description = config["platform_description"] + sensor_description = config["sensor_description"] + + logger.info(f'Processing {data_description}') + logger.info(f'reference_time = {reference_time}') + + bufrfile = f"{cycle_type}.t{hh}z.{data_type}.tm00.{data_format}" + DATA_PATH = os.path.join(dump_dir, f"{cycle_type}.{yyyymmdd}", str(hh), 'atmos', bufrfile) + if not os.path.isfile(DATA_PATH): + logger.info(f"DATA_PATH {DATA_PATH} does not exist") + return + logger.debug(f"The DATA_PATH is: {DATA_PATH}") + + # ============================================ + # Make the QuerySet for all the data we want + # ============================================ + start_time = time.time() + + logger.info('Making QuerySet') + q = bufr.QuerySet() + + # MetaData + q.add('latitude', '*/CLAT') + q.add('longitude', '*/CLON') + q.add('satelliteId', '*/SAID') + q.add('year', '*/YEAR') + q.add('month', '*/MNTH') + q.add('day', '*/DAYS') + q.add('hour', '*/HOUR') + q.add('minute', '*/MINU') + q.add('second', '*/SECO') + q.add('solarZenithAngle', '*/SOZA') + q.add('topLevelPressure', '*/OZOPQLSQ/PRLC[2]') + q.add('bottomLevelPressure', '*/OZOPQLSQ/PRLC[1]') + + # Quality + q.add('totalOzoneQuality', '*/SBUVTOQ') + q.add('profileOzoneQuality', '*/SBUVPOQ') + + # ObsValue + q.add('ozoneLayer', '*/OZOPQLSQ/OZOP[2]') + + end_time = time.time() + running_time = end_time - start_time + logger.debug(f'Processing time for making QuerySet : {running_time} seconds') + + # ============================================================== + # Open the BUFR file and execute the QuerySet to get ResultSet + # Use the ResultSet returned to get numpy arrays of the data + # ============================================================== + start_time = time.time() + + logger.info('Executing QuerySet to get ResultSet') + with bufr.File(DATA_PATH) as f: + try: + r = f.execute(q) + except Exception as err: + logger.info(f'Return with {err}') + return + + # MetaData + satid = r.get('satelliteId', 'ozoneLayer') + year = r.get('year', 'ozoneLayer') + month = r.get('month', 'ozoneLayer') + day = r.get('day', 'ozoneLayer') + hour = r.get('hour', 'ozoneLayer') + minute = r.get('minute', 'ozoneLayer') + second = r.get('second', 'ozoneLayer') + lat = r.get('latitude', 'ozoneLayer') + lon = r.get('longitude', 'ozoneLayer') + solzenang = r.get('solarZenithAngle', 'ozoneLayer') + ptop = r.get('topLevelPressure', 'ozoneLayer', type='float') + pbot = r.get('bottomLevelPressure', 'ozoneLayer', type='float') + + # Quality Information + poqc = r.get('profileOzoneQuality', 'ozoneLayer') + toqc = r.get('totalOzoneQuality', 'ozoneLayer') + + # ObsValue + # Total Ozone + o3val = r.get('ozoneLayer', 'ozoneLayer', type='float') + + # DateTime: seconds since Epoch time + # IODA has no support for numpy datetime arrays dtype=datetime64[s] + timestamp = r.get_datetime('year', 'month', 'day', 'hour', 'minute', 'second', 'ozoneLayer').astype(np.int64) + + end_time = time.time() + running_time = end_time - start_time + logger.debug(f'Processing time for executing QuerySet to get ResultSet : {running_time} seconds') + + # ========================= + # Create derived variables + # ========================= + start_time = time.time() + + # Set reference pressure + nlevs = 21 + nprofs = int(len(o3val)/nlevs) + pref0 = np.array([1000.000, 631.000, 398.000, 251.000, 158.000, 100.000, 63.100, + 39.800, 25.100, 15.800, 10.000, 6.310, 3.980, 2.510, + 1.580, 1.000, 0.631, 0.398, 0.251, 0.158, 0.100]) + pref0 = pref0 * 100. * 1.01325 + + nprofs = int(len(o3val)/nlevs) + pressure = np.tile(pref0, nprofs).astype(np.float32) + + # Set reference pressure vertices + prestop = np.array([0, 631.000, 398.000, 251.000, 158.000, 100.000, 63.100, + 39.800, 25.100, 15.800, 10.000, 6.310, 3.980, 2.510, + 1.580, 1.000, 0.631, 0.398, 0.251, 0.158, 0.100, 0]) + presbot = np.array([1000.000, 1000.000, 631.000, 398.000, 251.000, 158.000, 100.000, + 63.100, 39.800, 25.100, 15.800, 10.000, 6.310, 3.980, + 2.510, 1.580, 1.000, 0.631, 0.398, 0.251, 0.158, 0.100]) + prestop = prestop * 100. * 1.01325 + presbot = presbot * 100. * 1.01325 + + pres = np.zeros((nlevs+1, 2), dtype=np.float32) + pres[:, 0] = prestop + pres[:, 1] = presbot + presv = np.tile(pres, (nprofs, 1)).astype(np.float32) + + # Calculate total column ozone + o3tot = np.sum(o3val.reshape(-1, nlevs), axis=1) + o3tot[o3tot < 1.0e-20] = o3val.fill_value + + # Append total column zoone at the end of each profile + # Create a dictionary to store the arrays + arrays_dict = { + 'o3val': o3val, + 'lat': lat, + 'lon': lon, + 'timestamp': timestamp, + 'solzenang': solzenang, + 'satid': satid, + 'poqc': poqc, + 'toqc': toqc, + 'pressure': pressure, + 'ptop': ptop, + 'pbot': pbot, + } + + for key in arrays_dict.keys(): + array = arrays_dict[key] + + # Reshape the array (1D to 2D) + reshaped_array = array.reshape((nprofs, nlevs)) + + # Create a temporary array with zeros +# tmp = np.zeros((nprofs, 22)) + tmp = np.empty((nprofs, 22), dtype=array.dtype) + + # Assign the reshaped array values to the temporary array + tmp[:, 1:] = reshaped_array + + if key == 'o3val': + tmp[:, 0] = o3tot + elif key == 'pressure': + tmp[:, 0] = 0.0 + else: + tmp[:, 0] = tmp[:, 1] + + # Update the array in the dictionary + arrays_dict[key] = tmp + + # Update the variables with the modified arrays + lat1 = arrays_dict['lat'] + lon1 = arrays_dict['lon'] + o3val1 = arrays_dict['o3val'] + timestamp1 = arrays_dict['timestamp'] + solzenang1 = arrays_dict['solzenang'] + satid1 = arrays_dict['satid'] + poqc1 = arrays_dict['poqc'] + toqc1 = arrays_dict['toqc'] + pressure1 = arrays_dict['pressure'] + ptop1 = arrays_dict['ptop'] + pbot1 = arrays_dict['pbot'] + presv1 = presv.reshape(nprofs, 22, 2) + + # Print the profiles for debugging + for i in range(o3val1.shape[0]): + for k in reversed(range(o3val1.shape[1])): + row = np.array([i+1, k+1, lat1[i, k], lon1[i, k], ptop1[i, k], pbot1[i, k], pressure1[i, k], presv1[i, k, 0], presv1[i, k, 1], + o3val1[i, k], toqc1[i, k], poqc1[i, k], solzenang1[i, k]]) + frow = [format_element(x) for x in row.flatten()] + row = ' '.join(frow) + logger.debug(row) + + # Update the variables with the modified arrays + lat1 = arrays_dict['lat'].flatten() + lon1 = arrays_dict['lon'].flatten() + o3val1 = arrays_dict['o3val'].flatten() + timestamp1 = arrays_dict['timestamp'].flatten() + solzenang1 = arrays_dict['solzenang'].flatten() + satid1 = arrays_dict['satid'].flatten() + poqc1 = arrays_dict['poqc'].flatten() + toqc1 = arrays_dict['toqc'].flatten() + pressure1 = arrays_dict['pressure'].flatten() + ptop1 = arrays_dict['ptop'].flatten() + pbot1 = arrays_dict['pbot'].flatten() + presv1 = presv1.reshape(-1, presv1.shape[-1]) + + end_time = time.time() + running_time = end_time - start_time + logger.debug(f'Running time for creating derived variables : {running_time} seconds') + + # ===================================== + # Split output based on satellite id + # Create IODA ObsSpace + # Write IODA output + # ===================================== + logger.info('Split data based on satellite id, Create IODA ObsSpace and Write IODA output') + + # Find unique satellite identifiers in data to process + unique_satids = np.unique(satid) + logger.info(f'Number of Unique satellite identifiers : {len(unique_satids)}') + logger.info(f'Unique satellite identifiers: {unique_satids}') + + logger.debug(f'Loop through unique satellite identifier : {unique_satids}') + total_ob_processed = 0 + for sat in unique_satids.tolist(): + start_time = time.time() + + matched = False + for satellite_info in satellite_info_array: + if (satellite_info["satellite_id"] == sat): + matched = True + satellite_id = satellite_info["satellite_id"] + satellite_name = satellite_info["satellite_name"] + satinst = sensor_name.lower()+'_'+satellite_name.lower() + logger.debug(f'Split data for {satinst} satid = {sat}') + + if matched: + + # Define a boolean mask to subset data from the original data object + mask = satid1 == sat + # MetaData + lon2 = np.flip(lon1[mask], axis=0) + lat2 = np.flip(lat1[mask], axis=0) + satid2 = np.flip(satid1[mask], axis=0) + solzenang2 = np.flip(solzenang1[mask], axis=0) + timestamp2 = np.flip(timestamp1[mask], axis=0) + pressure2 = np.flip(pressure1[mask], axis=0) + presv2 = np.flip(presv1[mask], axis=0) + ptop2 = ptop1[mask] + pbot2 = pbot1[mask] + + # QC Info + toqc2 = np.flip(toqc1[mask], axis=0) + poqc2 = np.flip(poqc1[mask], axis=0) + + # ObsValue + o3val2 = np.flip(o3val1[mask], axis=0) + + timestamp2_min = datetime.fromtimestamp(timestamp1.min()) + timestamp2_max = datetime.fromtimestamp(timestamp1.max()) + + # Check unique observation time + unique_timestamp = np.unique(timestamp1) + + # Create the dimensions + dims = { + 'Location': np.arange(0, lat2.shape[0]), + 'Vertice': np.array([2, 1]) + } + + # Create IODA ObsSpace + sat = satellite_name.lower() + iodafile = f"{cycle_type}.t{hh}z.{ioda_type}_{sat}.tm00.nc" + OUTPUT_PATH = os.path.join(ioda_dir, iodafile) + obsspace = ioda_ospace.ObsSpace(OUTPUT_PATH, mode='w', dim_dict=dims) + logger.info(f'Create output file : {OUTPUT_PATH}') + + # Create Global attributes + logger.debug('Create global attributes') + obsspace.write_attr('sourceFiles', bufrfile) + obsspace.write_attr('source', source) + obsspace.write_attr('description', data_description) + obsspace.write_attr('datetimeReference', reference_time) + obsspace.write_attr('Converter', converter) + obsspace.write_attr('platformLongDescription', platform_description) + obsspace.write_attr('platformCommonName', satellite_name) + obsspace.write_attr('platform', satellite_id) + obsspace.write_attr('sensorLongDescription', sensor_description) + obsspace.write_attr('sensorCommonName', sensor_name) + obsspace.write_attr('sensor', sensor_id) + obsspace.write_attr('dataProviderOrigin', data_provider) + obsspace.write_attr('processingLevel', process_level) + obsspace.write_attr('datetimeRange', [str(timestamp2_min), str(timestamp2_max)]) + + # Create IODA variables + # Longitude + obsspace.create_var('MetaData/longitude', dtype=lon2.dtype, fillval=lon.fill_value) \ + .write_attr('units', 'degrees_east') \ + .write_attr('valid_range', np.array([-180, 180], dtype=np.float32)) \ + .write_attr('long_name', 'Longitude') \ + .write_data(lon2) + + # Latitude + obsspace.create_var('MetaData/latitude', dtype=lat2.dtype, fillval=lat.fill_value) \ + .write_attr('units', 'degrees_north') \ + .write_attr('valid_range', np.array([-90, 90], dtype=np.float32)) \ + .write_attr('long_name', 'Latitude') \ + .write_data(lat2) + + # Datetime + obsspace.create_var('MetaData/dateTime', dtype=timestamp2.dtype, fillval=timestamp.fill_value) \ + .write_attr('units', 'seconds since 1970-01-01T00:00:00Z') \ + .write_attr('long_name', 'Unix Epoch') \ + .write_data(timestamp2) + + # Satellite Identifier + obsspace.create_var('MetaData/satelliteIdentifier', dtype=satid2.dtype, fillval=satid.fill_value) \ + .write_attr('long_name', 'Satellite Identifier') \ + .write_data(satid2) + + # Total Ozone Quality + obsspace.create_var('MetaData/totalOzoneQuality', dtype=toqc2.dtype, fillval=toqc.fill_value) \ + .write_attr('long_name', 'Total Ozone Quality') \ + .write_data(toqc2) + + # Profile Ozone Quality + obsspace.create_var('MetaData/profileOzoneQuality', dtype=poqc2.dtype, fillval=poqc.fill_value) \ + .write_attr('long_name', 'Layer Ozone Quality') \ + .write_data(poqc2) + + # Pressure + obsspace.create_var('MetaData/pressure', dtype=pressure2.dtype, fillval=lat.fill_value) \ + .write_attr('units', 'pa') \ + .write_attr('long_name', 'Pressure') \ + .write_data(pressure2) + + # Pressure vertices + obsspace.create_var('RetrievalAncillaryData/pressureVertice', dim_list=['Location', 'Vertice'], dtype=pbot2.dtype, fillval=lat.fill_value) \ + .write_attr('units', 'pa') \ + .write_attr('long_name', 'Retrieval Pressure Vertices') \ + .write_data(presv2) + + # Top-level Pressure + obsspace.create_var('MetaData/topLevelPressure', dtype=ptop2.dtype, fillval=lat.fill_value) \ + .write_attr('units', 'pa') \ + .write_attr('long_name', 'Top Level Pressure') \ + .write_data(ptop2) + + # Bottom-level Pressure + obsspace.create_var('MetaData/bottomLevelPressure', dtype=pbot2.dtype, fillval=lat.fill_value) \ + .write_attr('units', 'pa') \ + .write_attr('long_name', 'Bottom Level Pressure') \ + .write_data(pbot2) + + # Solar Zenith Angle + obsspace.create_var('MetaData/solarZenithAngle', dtype=solzenang2.dtype, fillval=solzenang.fill_value) \ + .write_attr('units', 'degree') \ + .write_attr('valid_range', np.array([0, 180], dtype=np.float32)) \ + .write_attr('long_name', 'Solar Zenith Angle') \ + .write_data(solzenang2) + + # Layer Ozone and Total Ozone + obsspace.create_var('ObsValue/ozoneLayer', dtype=o3val2.dtype, fillval=o3val.fill_value) \ + .write_attr('units', 'DU') \ + .write_attr('long_name', 'Layer Ozone') \ + .write_data(o3val2) + + else: + logger.info(f'Do not find this satellite id in the configuration: satid = {sat}') + + logger.info('All Done!') + logger.info(f'Total number of observation processed : {total_ob_processed}') + + +if __name__ == '__main__': + + start_time = time.time() + + parser = argparse.ArgumentParser() + parser.add_argument('-c', '--config', type=str, help='Input JSON configuration', required=True) + parser.add_argument('-v', '--verbose', help='print debug logging information', + action='store_true') + args = parser.parse_args() + + log_level = 'DEBUG' if args.verbose else 'INFO' + logger = Logger('BUFR2IODA_ozone_ompsnp.py', level=log_level, colored_log=True) + + with open(args.config, "r") as json_file: + config = json.load(json_file) + + bufr_to_ioda(config, logger) + + end_time = time.time() + running_time = end_time - start_time + logger.info(f"Total running time : {running_time} seconds") diff --git a/ush/ioda/bufr2ioda/bufr2ioda_ozone_ompstc.py b/ush/ioda/bufr2ioda/bufr2ioda_ozone_ompstc.py new file mode 100755 index 000000000..12359e9a2 --- /dev/null +++ b/ush/ioda/bufr2ioda/bufr2ioda_ozone_ompstc.py @@ -0,0 +1,328 @@ +#!/usr/bin/env python3 +import argparse +import numpy as np +import numpy.ma as ma +from pyiodaconv import bufr +import calendar +import json +import time +import math +import datetime +import os +from datetime import datetime +from pyioda import ioda_obs_space as ioda_ospace +from wxflow import Logger + +# ================================================================================================== +# Subset | Description (OMPS TC) | +# -------------------------------------------------------------------------------------------------- +# NC008018 | OMPS TC is the Total Column Ozone derived from the ultraviolet, visible and | +# | near-infrared observations (between 300-380 nm) made by the OMPS Nadir Total Column | +# | Mapper on Suomi-NPP and JPSS at a spatial resolution of 50 km at nadir. | +# ================================================================================================== + + +def bufr_to_ioda(config, logger): + + # Get parameters from configuration + subsets = config["subsets"] + source = config["source"] + data_format = config["data_format"] + data_type = config["data_type"] + data_description = config["data_description"] + data_provider = config["data_provider"] + cycle_type = config["cycle_type"] + dump_dir = config["dump_directory"] + ioda_dir = config["ioda_directory"] + ioda_type = config["ioda_type"] + cycle = config["cycle_datetime"] + yyyymmdd = cycle[0:8] + hh = cycle[8:10] + + satellite_info_array = config["satellite_info"] + sensor_name = config["sensor_info"]["sensor_name"] + sensor_full_name = config["sensor_info"]["sensor_full_name"] + sensor_id = config["sensor_info"]["sensor_id"] + + # Get derived parameters + yyyymmdd = cycle[0:8] + hh = cycle[8:10] + reference_time = datetime.strptime(cycle, "%Y%m%d%H") + reference_time = reference_time.strftime("%Y-%m-%dT%H:%M:%SZ") + + # General informaton for global attribute + converter = 'BUFR to IODA Converter' + process_level = config["process_level"] + platform_description = config["platform_description"] + sensor_description = config["sensor_description"] + + logger.info(f'Processing {data_description}') + logger.info(f'reference_time = {reference_time}') + + bufrfile = f"{cycle_type}.t{hh}z.{data_type}.tm00.{data_format}" + DATA_PATH = os.path.join(dump_dir, f"{cycle_type}.{yyyymmdd}", str(hh), 'atmos', bufrfile) + if not os.path.isfile(DATA_PATH): + logger.info(f"DATA_PATH {DATA_PATH} does not exist") + return + logger.debug(f"The DATA_PATH is: {DATA_PATH}") + + # ============================================ + # Make the QuerySet for all the data we want + # ============================================ + start_time = time.time() + + logger.info('Making QuerySet') + q = bufr.QuerySet() + + # MetaData + q.add('latitude', '*/CLAT') + q.add('longitude', '*/CLON') + q.add('satelliteId', '*/SAID') + q.add('year', '*/YEAR') + q.add('month', '*/MNTH') + q.add('day', '*/DAYS') + q.add('hour', '*/HOUR') + q.add('minute', '*/MINU') + q.add('second', '*/SECO') + q.add('solarZenithAngle', '*/SOZA') + q.add('sensorScanPosition', '*/FOVN') + + # Quality + q.add('bestOzoneAlgorithmFlag', '*/AFBO') + q.add('qualityFlags', '*/TOQC') + + # ObsValue + q.add('ozoneTotal', '*/OZON') + + end_time = time.time() + running_time = end_time - start_time + logger.debug(f'Processing time for making QuerySet : {running_time} seconds') + + # ============================================================== + # Open the BUFR file and execute the QuerySet to get ResultSet + # Use the ResultSet returned to get numpy arrays of the data + # ============================================================== + start_time = time.time() + + logger.info('Executing QuerySet to get ResultSet') + with bufr.File(DATA_PATH) as f: + try: + r = f.execute(q) + except Exception as err: + logger.info(f'Return with {err}') + return + + # MetaData + satid = r.get('satelliteId') + year = r.get('year') + month = r.get('month') + day = r.get('day') + hour = r.get('hour') + minute = r.get('minute') + second = r.get('second') + lat = r.get('latitude') + lon = r.get('longitude') + solzenang = r.get('solarZenithAngle') + scanpos = r.get('sensorScanPosition') + + # Quality Information + toqc = r.get('qualityFlags') + afbo = r.get('bestOzoneAlgorithmFlag') + + # ObsValue + # Total Ozone + o3val = r.get('ozoneTotal', type='float') + + # DateTime: seconds since Epoch time + # IODA has no support for numpy datetime arrays dtype=datetime64[s] + timestamp = r.get_datetime('year', 'month', 'day', 'hour', 'minute', 'second').astype(np.int64) + + end_time = time.time() + running_time = end_time - start_time + logger.debug(f'Processing time for executing QuerySet to get ResultSet : {running_time} seconds') + + # ========================= + # Create derived variables + # ========================= + start_time = time.time() + + # Create pressure variable and fill in with zeros + pressure = np.ma.array(np.zeros(lat.shape[0], dtype=np.float32)) + pressure.mask = lat.mask + pressure.fill_value = lat.fill_value + + end_time = time.time() + running_time = end_time - start_time + logger.debug(f'Running time for creating derived variables : {running_time} seconds') + + # ===================================== + # Split output based on satellite id + # Create IODA ObsSpace + # Write IODA output + # ===================================== + logger.info('Split data based on satellite id, Create IODA ObsSpace and Write IODA output') + + # Find unique satellite identifiers in data to process + unique_satids = np.unique(satid) + logger.info(f'Number of Unique satellite identifiers : {len(unique_satids)}') + logger.info(f'Unique satellite identifiers: {unique_satids}') + + logger.debug(f'Loop through unique satellite identifier : {unique_satids}') + total_ob_processed = 0 + for sat in unique_satids.tolist(): + start_time = time.time() + + matched = False + for satellite_info in satellite_info_array: + if (satellite_info["satellite_id"] == sat): + matched = True + satellite_id = satellite_info["satellite_id"] + satellite_name = satellite_info["satellite_name"] + satinst = sensor_name.lower()+'_'+satellite_name.lower() + logger.debug(f'Split data for {satinst} satid = {sat}') + + if matched: + + # Define a boolean mask to subset data from the original data object + mask = satid == sat + # MetaData + lon2 = lon[mask] + lat2 = lat[mask] + satid2 = satid[mask] + solzenang2 = solzenang[mask] + scanpos2 = scanpos[mask] + timestamp2 = timestamp[mask] + pressure2 = pressure[mask] + + # QC Info + toqc2 = toqc[mask] + afbo2 = afbo[mask] + + # ObsValue + o3val2 = o3val[mask] + + timestamp2_min = datetime.fromtimestamp(timestamp2.min()) + timestamp2_max = datetime.fromtimestamp(timestamp2.max()) + + # Create the dimensions + dims = { + 'Location': np.arange(0, lat2.shape[0]), + } + + # Create IODA ObsSpace + sat = satellite_name.lower() + iodafile = f"{cycle_type}.t{hh}z.{ioda_type}_{sat}.tm00.nc" + OUTPUT_PATH = os.path.join(ioda_dir, iodafile) + obsspace = ioda_ospace.ObsSpace(OUTPUT_PATH, mode='w', dim_dict=dims) + logger.info(f'Create output file : {OUTPUT_PATH}') + + # Create Global attributes + logger.debug('Create global attributes') + obsspace.write_attr('sourceFiles', bufrfile) + obsspace.write_attr('source', source) + obsspace.write_attr('description', data_description) + obsspace.write_attr('datetimeReference', reference_time) + obsspace.write_attr('Converter', converter) + obsspace.write_attr('platformLongDescription', platform_description) + obsspace.write_attr('platformCommonName', satellite_name) + obsspace.write_attr('platform', satellite_id) + obsspace.write_attr('sensorLongDescription', sensor_description) + obsspace.write_attr('sensorCommonName', sensor_name) + obsspace.write_attr('sensor', sensor_id) + obsspace.write_attr('dataProviderOrigin', data_provider) + obsspace.write_attr('processingLevel', process_level) + obsspace.write_attr('datetimeRange', [str(timestamp2_min), str(timestamp2_max)]) + + # Create IODA variables + logger.debug('Create variables: name, type, units, and attributes') + # Longitude + obsspace.create_var('MetaData/longitude', dtype=lon2.dtype, fillval=lon2.fill_value) \ + .write_attr('units', 'degrees_east') \ + .write_attr('valid_range', np.array([-180, 180], dtype=np.float32)) \ + .write_attr('long_name', 'Longitude') \ + .write_data(lon2) + + # Latitude + obsspace.create_var('MetaData/latitude', dtype=lat2.dtype, fillval=lat2.fill_value) \ + .write_attr('units', 'degrees_north') \ + .write_attr('valid_range', np.array([-90, 90], dtype=np.float32)) \ + .write_attr('long_name', 'Latitude') \ + .write_data(lat2) + + # Datetime + obsspace.create_var('MetaData/dateTime', dtype=timestamp2.dtype, fillval=timestamp2.fill_value) \ + .write_attr('units', 'seconds since 1970-01-01T00:00:00Z') \ + .write_attr('long_name', 'Unix Epoch') \ + .write_data(timestamp2) + + # Satellite Identifier + obsspace.create_var('MetaData/satelliteIdentifier', dtype=satid2.dtype, fillval=satid2.fill_value) \ + .write_attr('long_name', 'Satellite Identifier') \ + .write_data(satid2) + + # Total Ozone Quality Code + obsspace.create_var('MetaData/totalOzoneQualityCode', dtype=toqc2.dtype, fillval=toqc2.fill_value) \ + .write_attr('long_name', 'OMI Total Ozone Quality Code') \ + .write_data(toqc2) + + # Pressure + obsspace.create_var('MetaData/pressure', dtype=pressure2.dtype, fillval=pressure2.fill_value) \ + .write_attr('units', 'pa') \ + .write_attr('long_name', 'Pressure') \ + .write_data(pressure2) + + # Algorithm Flag for Best Ozone + obsspace.create_var('MetaData/bestOzoneAlgorithmFlag', dtype=afbo2.dtype, fillval=afbo2.fill_value) \ + .write_attr('long_name', 'Algorithm Flag for Best Ozone') \ + .write_data(afbo2) + + # Solar Zenith Angle + obsspace.create_var('MetaData/solarZenithAngle', dtype=solzenang2.dtype, fillval=solzenang2.fill_value) \ + .write_attr('units', 'm') \ + .write_attr('long_name', 'Solar Zenith Angle') \ + .write_data(solzenang2) + + # Sensor Scan Position + obsspace.create_var('MetaData/sensorScanPosition', dtype=scanpos2.dtype, fillval=scanpos2.fill_value) \ + .write_attr('long_name', 'Sensor Scan Position') \ + .write_data(scanpos2) + + # Total Ozone + obsspace.create_var('ObsValue/ozoneTotal', dtype=o3val2.dtype, fillval=o3val2.fill_value) \ + .write_attr('units', 'DU') \ + .write_attr('long_name', 'Total Column Ozone') \ + .write_data(o3val2) + + end_time = time.time() + running_time = end_time - start_time + total_ob_processed += len(satid2) + logger.debug(f'Processing time for splitting and output IODA for {satinst} : {running_time}, seconds') + + else: + logger.info(f'Do not find this satellite id in the configuration: satid = {sat}') + + logger.info('All Done!') + logger.info(f'Total number of observation processed : {total_ob_processed}') + + +if __name__ == '__main__': + + start_time = time.time() + + parser = argparse.ArgumentParser() + parser.add_argument('-c', '--config', type=str, help='Input JSON configuration', required=True) + parser.add_argument('-v', '--verbose', help='print debug logging information', + action='store_true') + args = parser.parse_args() + + log_level = 'DEBUG' if args.verbose else 'INFO' + logger = Logger('BUFR2IODA_ozone_ompstc8.py', level=log_level, colored_log=True) + + with open(args.config, "r") as json_file: + config = json.load(json_file) + + bufr_to_ioda(config, logger) + + end_time = time.time() + running_time = end_time - start_time + logger.info(f"Total running time : {running_time} seconds") From a575b46b204eb111fbd012f9a138f87169f12f59 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Fri, 22 Mar 2024 11:52:35 -0400 Subject: [PATCH 2/3] Fix to marine DA prep (#989) --- scripts/exgdas_global_marine_analysis_prep.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/exgdas_global_marine_analysis_prep.py b/scripts/exgdas_global_marine_analysis_prep.py index 5e37d3114..171ea8154 100755 --- a/scripts/exgdas_global_marine_analysis_prep.py +++ b/scripts/exgdas_global_marine_analysis_prep.py @@ -271,7 +271,7 @@ def find_clim_ens(input_date): varconfig = Template.substitute_structure(varconfig, TemplateConstants.DOLLAR_PARENTHESES, envconfig.get) # Remove empty obs spaces in var_yaml -ufsda.yamltools.save_check(varconfig.as_dict(), target=var_yaml, app='var') +ufsda.yamltools.save_check(varconfig, target=var_yaml, app='var') ################################################################################ # Prepare the yamls for the "checkpoint" jjob From 861ae2b86d3cfbaeec9866f1a41d94f83170757c Mon Sep 17 00:00:00 2001 From: RussTreadon-NOAA <26926959+RussTreadon-NOAA@users.noreply.github.com> Date: Fri, 22 Mar 2024 14:28:57 -0400 Subject: [PATCH 3/3] restore test_gdasapp_atm_jjob functionality (#981) `test_gdasapp_atm_jjob` tests fail due to inconsistencies between the GDASApp tests and recent g-w updates. This PR restores _Passed_ result for `test_gdasapp_atm_jjob`. Fixes #980 --- test/atm/global-workflow/config.atmanl | 18 +++++++++++------- test/atm/global-workflow/config.atmensanl | 7 +++++-- test/atm/global-workflow/config.yaml | 7 +++---- 3 files changed, 19 insertions(+), 13 deletions(-) diff --git a/test/atm/global-workflow/config.atmanl b/test/atm/global-workflow/config.atmanl index 4aa60d76e..5649f70bb 100755 --- a/test/atm/global-workflow/config.atmanl +++ b/test/atm/global-workflow/config.atmanl @@ -5,17 +5,22 @@ echo "BEGIN: config.atmanl" +export OBS_LIST=@OBS_LIST@ +export JEDIYAML="${PARMgfs}/gdas/atm/variational/3dvar_drpcg.yaml.j2" +export STATICB_TYPE=@STATICB_TYPE@ +export INTERP_METHOD='barycentric' + if [[ ${DOHYBVAR} = "YES" ]]; then # shellcheck disable=SC2153 export CASE_ANL=${CASE_ENS} + export BERROR_YAML="${PARMgfs}/gdas/atm/berror/hybvar_${STATICB_TYPE}.yaml.j2" else export CASE_ANL=${CASE} + export BERROR_YAML="${PARMgfs}/gdas/atm/berror/staticb_${STATICB_TYPE}.yaml.j2" fi -export OBS_LIST=@OBS_LIST@ -export JEDIYAML="${HOMEgfs}/parm/gdas/atm/variational/3dvar_drpcg.yaml.j2" -export STATICB_TYPE="identity" -export BERROR_YAML="${HOMEgfs}/parm/gdas/atm/berror/staticb_${STATICB_TYPE}.yaml.j2" -export INTERP_METHOD='barycentric' + +export CRTM_FIX_YAML="${PARMgfs}/gdas/atm_crtm_coeff.yaml.j2" +export JEDI_FIX_YAML="${PARMgfs}/gdas/atm_jedi_fix.yaml.j2" export layout_x_atmanl=@LAYOUT_X_ATMANL@ export layout_y_atmanl=@LAYOUT_Y_ATMANL@ @@ -23,7 +28,6 @@ export layout_y_atmanl=@LAYOUT_Y_ATMANL@ export io_layout_x=@IO_LAYOUT_X@ export io_layout_y=@IO_LAYOUT_Y@ -export JEDIEXE=${HOMEgfs}/exec/fv3jedi_var.x -export crtm_VERSION="2.4.0" +export JEDIEXE=${EXECgfs}/fv3jedi_var.x echo "END: config.atmanl" diff --git a/test/atm/global-workflow/config.atmensanl b/test/atm/global-workflow/config.atmensanl index 86640782c..c7bddbdf6 100644 --- a/test/atm/global-workflow/config.atmensanl +++ b/test/atm/global-workflow/config.atmensanl @@ -6,15 +6,18 @@ echo "BEGIN: config.atmensanl" export OBS_LIST=@OBS_LIST@ -export JEDIYAML="${HOMEgfs}/parm/gdas/atm/lgetkf/lgetkf.yaml.j2" +export JEDIYAML="${PARMgfs}/gdas/atm/lgetkf/lgetkf.yaml.j2" export INTERP_METHOD='barycentric' +export CRTM_FIX_YAML="${PARMgfs}/gdas/atm_crtm_coeff.yaml.j2" +export JEDI_FIX_YAML="${PARMgfs}/gdas/atm_jedi_fix.yaml.j2" + export layout_x_atmensanl=@LAYOUT_X_ATMENSANL@ export layout_y_atmensanl=@LAYOUT_Y_ATMENSANL@ export io_layout_x=@IO_LAYOUT_X@ export io_layout_y=@IO_LAYOUT_Y@ -export JEDIEXE=${HOMEgfs}/exec/fv3jedi_letkf.x +export JEDIEXE=${EXECgfs}/fv3jedi_letkf.x echo "END: config.atmensanl" diff --git a/test/atm/global-workflow/config.yaml b/test/atm/global-workflow/config.yaml index b7594f4bf..73f2852c9 100644 --- a/test/atm/global-workflow/config.yaml +++ b/test/atm/global-workflow/config.yaml @@ -1,22 +1,21 @@ base: HOMEgfs: "@topdir@" + DOHYBVAR: "NO" DO_JEDIATMVAR: "YES" DO_JEDIATMENS: "YES" DATAPATH: "@bindir@/test/atm/global-workflow/testrun" DUMPDIR: "@dumpdir@" STMP: "@bindir@/test/atm/global-workflow/testrun" PTMP: "@bindir@/test/atm/global-workflow/testrun" + atmanl: OBS_LIST: "@srcdir@/test/atm/global-workflow/gdas_prototype.yaml.j2" + STATICB_TYPE: "identity" ATMRES_ANL: "C48" LAYOUT_X_ATMANL: 1 LAYOUT_Y_ATMANL: 1 - IO_LAYOUT_X: 1 - IO_LAYOUT_Y: 1 atmensanl: OBS_LIST: "@srcdir@/test/atm/global-workflow/lgetkf_prototype.yaml.j2" LAYOUT_X_ATMENSANL: 1 LAYOUT_Y_ATMENSANL: 1 - IO_LAYOUT_X: 1 - IO_LAYOUT_Y: 1