From ee2b266e9bfc7c30216eb88284093ce68ad7b71f Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Wed, 24 Apr 2024 10:56:31 -0400 Subject: [PATCH 01/10] No mpi for the ascii stats (#1070) Maybe it's an easy fix ... But it takes 2 seconds to sequentially compute all the stats for all the obs spaces that we currently run with. Good enough! --- utils/soca/gdassoca_obsstats.h | 28 ++++++---------------------- 1 file changed, 6 insertions(+), 22 deletions(-) diff --git a/utils/soca/gdassoca_obsstats.h b/utils/soca/gdassoca_obsstats.h index 98e70bf3b..5f95574a1 100644 --- a/utils/soca/gdassoca_obsstats.h +++ b/utils/soca/gdassoca_obsstats.h @@ -52,32 +52,18 @@ namespace gdasapp { std::vector obsSpaces; fullConfig.get("obs spaces", obsSpaces); - // check if the number of workers is correct. Only the serial and 1 diag file - // per pe works for now - ASSERT(getComm().size() <= 1 || getComm().size() >= obsSpaces.size()); - - // initialize pe's color and communicator's name - int color(0); - std::string commNameStr = "comm_idle"; - if (this->getComm().rank() < obsSpaces.size()) { - color = 1; - std::string commNameStr = "comm_work"; - } - - // Create the new communicator () - char const *commName = commNameStr.c_str(); - eckit::mpi::Comm & commObsSpace = this->getComm().split(color, commName); + // only the serial case works for now. + ASSERT(getComm().size() == 1); - // spread out the work over pes - if (color > 0) { + for (int i = 0; i < obsSpaces.size(); i++) { // get the obs space configuration - auto obsSpace = obsSpaces[commObsSpace.rank()]; + auto obsSpace = obsSpaces[i]; eckit::LocalConfiguration obsConfig(obsSpace, "obs space"); // get the obs diag file std::string obsFile; obsConfig.get("obsdatain.engine.obsfile", obsFile); - oops::Log::info() << "========= Processing" << obsFile + oops::Log::info() << "========= Processing " << obsFile << " date: " << extractDateFromFilename(obsFile) << std::endl; @@ -86,8 +72,7 @@ namespace gdasapp { obsSpace.get("variable", variable); // read the obs space - ioda::ObsSpace ospace(obsConfig, commObsSpace, timeWindow, - oops::mpi::myself()); + ioda::ObsSpace ospace(obsConfig, getComm(), timeWindow, getComm()); const size_t nlocs = ospace.nlocs(); oops::Log::info() << "nlocs =" << nlocs << std::endl; std::vector var(nlocs); @@ -125,7 +110,6 @@ namespace gdasapp { // Close the file outputFile.close(); } - getComm().barrier(); return EXIT_SUCCESS; } From c9029839500368f1afdd5bd02de14f00ecc5a8cf Mon Sep 17 00:00:00 2001 From: BrettHoover-NOAA <98188219+BrettHoover-NOAA@users.noreply.github.com> Date: Thu, 25 Apr 2024 10:10:47 -0500 Subject: [PATCH 02/10] Added YAML, JSON, python files for assimilating VIIRS satwinds (#1055) Adding satwinds from the Visible Infrared Imaging Radiometer Suite (VIIRS) from SNPP/NOAA-20 to GDASApp end-to-end testing new files include: parm/atm/obs/config/satwind_viirs_npp.yaml.j2: QC filter YAML for VIIRS SNPP satwinds (jinja2 standard) parm/atm/obs/config/satwind_viirs_n20.yaml.j2: QC filter YAML for VIIRS NOAA-20 satwinds (jinja2 standard) parm/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.json: JSON containing data format, sensor, and satellite information for VIIRS SNPP/NOAA-20 satwinds ush/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.py: bufr2ioda code for extracting VIIRS SNPP/NOAA-20 satwinds from BUFR See #1054 for end-to-end testing results. Acceptance and ob-errors agree well with some expected deviation. No thinning is applied to these winds by regular convention. Note: We are still using `qualityInformationWithoutForecast` as the variable-name for QI in the IODA converter. This variable-name is currently not registered in the IODA ObsSpace.yaml, but there is an ongoing discussion with JCSDA to have it added, issue is here: https://github.com/JCSDA-internal/ioda/issues/1233 Co-authored-by: Brett Hoover --- parm/atm/obs/config/satwind_viirs_n20.yaml.j2 | 289 +++++++++++ parm/atm/obs/config/satwind_viirs_npp.yaml.j2 | 289 +++++++++++ .../bufr2ioda_satwind_amv_viirs.json | 16 + .../bufr2ioda/bufr2ioda_satwind_amv_viirs.py | 470 ++++++++++++++++++ 4 files changed, 1064 insertions(+) create mode 100644 parm/atm/obs/config/satwind_viirs_n20.yaml.j2 create mode 100644 parm/atm/obs/config/satwind_viirs_npp.yaml.j2 create mode 100644 parm/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.json create mode 100755 ush/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.py diff --git a/parm/atm/obs/config/satwind_viirs_n20.yaml.j2 b/parm/atm/obs/config/satwind_viirs_n20.yaml.j2 new file mode 100644 index 000000000..2a0c2258a --- /dev/null +++ b/parm/atm/obs/config/satwind_viirs_n20.yaml.j2 @@ -0,0 +1,289 @@ +- obs space: + name: satwind_viirs_n20 + obsdatain: + engine: + type: H5File + obsfile: '{{ DATA }}/obs/{{ OPREFIX }}satwnd.viirs_n20.tm00.nc' + obsdataout: + engine: + type: H5File + obsfile: '{{ DATA }}/diags/diag_satwind_viirs_n20_{{ current_cycle | to_YMDH }}.nc' + io pool: + max pool size: 1 + simulated variables: [windEastward, windNorthward] + + obs operator: + name: VertInterp + hofx scaling field: SurfaceWindScalingPressure + hofx scaling field group: DerivedVariables + + linear obs operator: + name: VertInterp + + # NOTE: Tests using the Gaussian Thinning filter (below) to duplicate GSI's thinning of AHI/Himawari-8 satwinds + # results in more JEDI satwinds in the diag file than in GSI, but far fewer JEDI satwinds assimilated than + # GSI. JEDI under-counts assimilated winds by roughly 25-40%, relative to GSI, and this under-count is not + # even including the temporal thinning which is applied in GSI but not JEDI (by this filter below). See + # GDASApp Issue #741 for details: https://github.com/NOAA-EMC/GDASApp/issues/741 + #obs pre filters: + #- filter: Gaussian Thinning + # horizontal_mesh: 200 + # vertical_mesh: 10000 + # use_reduced_horizontal_grid: true + # round_horizontal_bin_count_to_nearest: true + # partition_longitude_bins_using_mesh: true + + obs prior filters: + # Apply variable changes needed for wind scaling + # For wind observations with pressure provided + - filter: Variable Transforms + Transform: SurfaceWindScalingPressure + SkipWhenNoObs: False + + # Calculate error inflation factor for duplicate observations + #- filter: Variable Assignment + # assignments: + # - name: ObsErrorFactorDuplicateCheck/windEastward + # type: float + # function: + # name: ObsFunction/ObsErrorFactorDuplicateCheck + # options: + # use_air_pressure: true + # variable: windEastward + + #- filter: Variable Assignment + # assignments: + # - name: ObsErrorFactorDuplicateCheck/windNorthward + # type: float + # function: + # name: ObsFunction/ObsErrorFactorDuplicateCheck + # options: + # use_air_pressure: true + # variable: windNorthward + + obs post filters: + # Assign the initial observation error, based on height/pressure + # Hard-wiring to prepobs_errtable.global by Type + # ObsError is currently not updating in diag file, but passes directly to EffectiveError when no inflation is specified in YAML + + # Type 260 (VIIRS LWIR) + - filter: Perform Action + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 260 + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + xvar: + name: MetaData/pressure + xvals: [110000.,105000.,100000.,95000.,90000.,85000.,80000.,75000.,70000.,65000.,60000.,55000.,50000.,45000.,40000.,35000.,30000.,25000.,20000.,15000.,10000.,7500.,5000.,4000.,3000.,2000.,1000.,500.,400.,300.,200.,100.,0.] + errors: [3.8,3.8,3.8,3.8,3.8,3.8,3.8,3.8,3.9,3.9,4.,4.,4.1,5.,6.,6.3,6.6,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.] + + # sanity-check criteria + # Observation Range Sanity Check + # NOT EXPLICITLY CLEARED: No obs in this range in file, so 0 Bounds Check rejects (which is correct) but essentially untested + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + minvalue: -130. + maxvalue: 130. + action: + name: reject + + # Velocity Sanity Check + # NOT EXPLICITLY CLEARED: No obs in this range in file, so 0 Bounds Check rejects (which is correct) but essentially untested + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/Velocity + maxvalue: 130. + action: + name: reject + + # Reject any observation with a /=0 surface type (non-water surface) within + # 200 hPa of the surface pressure (as part of the LNVD check). + # CLEARED: maxvalue is rejecting >, not >= as per a Perform Action, so threshold is unchanged + - filter: Difference Check + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: + name: GeoVaLs/water_area_fraction + maxvalue: 0.99 + reference: GeoVaLs/surface_pressure + value: MetaData/pressure + maxvalue: -20000. # within 200 hPa above surface pressure, negative p-diff + action: + name: reject + + # LNVD check + # CLEARED: maxvalue is rejecting >, not >= as per a Perform Action, so threshold is unchanged + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/SatWindsLNVDCheck + maxvalue: 3. + action: + name: reject + + # GSI setupw routine QC + # Reject any ob Type [240–260] when pressure greater than 950 mb. + # CLEARED: minvalue/maxvalue are >=/<=, not >/<, so editing range by 1 Pa + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 240-260 + test variables: + - name: MetaData/pressure + maxvalue: 95001. + action: + name: reject + + # Multiple satellite platforms, reject when pressure is more than 50 mb above tropopause. + # CLEARED: minvalue is rejecting <, not <= as per a Perform Action, so threshold is unchanged + # Notes (eliu): This tropopause check reject too many obs; probably due to tropopause pressure estimation + # Turn this check off for now. + # Need to check if troposphere pressure was implemented correctly in fv3-jed + - filter: Difference Check + filter variables: + - name: windEastward + - name: windNorthward + reference: GeoVaLs/tropopause_pressure + value: MetaData/pressure + minvalue: -5000. # 50 hPa above tropopause level, negative p-diff + action: + name: reject + + # All satwinds must adjust errors based on ObsErrorFactorPressureCheck + # prior to the SPDB check (i.e. the gross-error check). The gross-error + # check uses the adjusted errors for error-bound tightening and rejection, + # so this check has to come first. This check will inflate errors for obs + # that are too close to either the model top or bottom. + # Notes (eliu): GMAO added a required parameter: adjusted_error_name. + - filter: Perform Action + filter variables: + - name: windEastward + where: + - variable: + name: ObsType/windEastward + is_in: 240-260 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck + options: + surface_obs: false + variable: windEastward + inflation factor: 4.0 + + - filter: Perform Action + filter variables: + - name: windNorthward + where: + - variable: + name: ObsType/windNorthward + is_in: 240-260 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck + options: + variable: windNorthward + inflation factor: 4.0 + + # All satwinds subject to a gross-error check that contains significant + # modifiers for satwinds with a negative speed-bias. ALL wind gross-error + # checks are currently being done by the SatWindsSPDBCheck. + # CLEARED + - filter: Background Check + filter variables: + - name: windEastward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260] + cgross: [ 2.5, 2.5, 2.5, 1.5, 2.5, 1.3, 1.3, 2.5, 2.5, 2.5, 2.5, 1.3, 2.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5] + error_min: [1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4] + error_max: [6.1, 6.1, 15.0, 15.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.1, 20.1, 20.1, 20.1, 20.1, 20.1] + variable: windEastward + action: + name: reject + + - filter: Background Check + filter variables: + - name: windNorthward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260] + cgross: [ 2.5, 2.5, 2.5, 1.5, 2.5, 1.3, 1.3, 2.5, 2.5, 2.5, 2.5, 1.3, 2.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5] + error_min: [1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4] + error_max: [6.1, 6.1, 15.0, 15.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.1, 20.1, 20.1, 20.1, 20.1, 20.1] + variable: windNorthward + action: + name: reject + + # The last error inflation check is for duplicate observations. This one needs + # to come last, because we don't want to inflate errors for duplication if one + # of the duplicates should be rejected. + # Notes (eliu): ObsErrorFactorDuplicateCheck obsfunction requires PreUseFlag (usage parameter from read_satwnd.f90). + # : Turn off duplicate check for now. + #- filter: Perform Action + # filter variables: + # - name: windEastward + # action: + # name: inflate error + # inflation variable: + # name: ObsErrorFactorDuplicateCheck/windEastward + + #- filter: Perform Action + # filter variables: + # - name: windNorthward + # action: + # name: inflate error + # inflation variable: + # name: ObsErrorFactorDuplicateCheck/windNorthward + + # We are extending this to an additional filter that inflates final ob-errors across-the-board by + # 1/0.8 = 1.25. This is caused by the GSI value of nvqc being set to .true. in the global operational + # configuration, see: https://github.com/NOAA-EMC/global-workflow/blob/d5ae3328fa4041b177357b1133f6b92e81c859d7/scripts/exglobal_atmos_analysis.sh#L750 + # This setting activates Line 1229 of setupw.f90 to scale ratio_errors by 0.8, which is applied in + # the denominator of the final ob-error, so 1/0.8 = 1.25 factor of ob-error inflation. + # + # If this nvqc functionality were to be switched off (i.e. if variational qc were to be turned off), + # you would want to remove this last inflation filter. + #- filter: Perform Action + # filter variables: + # - name: windEastward + # where: + # - variable: ObsType/windEastward + # is_in: 240-260 + # action: + # name: inflate error + # inflation factor: 1.25 + + #- filter: Perform Action + # filter variables: + # - name: windNorthward + # where: + # - variable: ObsType/windNorthward + # is_in: 240-260 + # action: + # name: inflate error + # inflation factor: 1.25 + + # End of Filters diff --git a/parm/atm/obs/config/satwind_viirs_npp.yaml.j2 b/parm/atm/obs/config/satwind_viirs_npp.yaml.j2 new file mode 100644 index 000000000..ce8a377a2 --- /dev/null +++ b/parm/atm/obs/config/satwind_viirs_npp.yaml.j2 @@ -0,0 +1,289 @@ +- obs space: + name: satwind_viirs_npp + obsdatain: + engine: + type: H5File + obsfile: '{{ DATA }}/obs/{{ OPREFIX }}satwnd.viirs_npp.tm00.nc' + obsdataout: + engine: + type: H5File + obsfile: '{{ DATA }}/diags/diag_satwind_viirs_npp_{{ current_cycle | to_YMDH }}.nc' + io pool: + max pool size: 1 + simulated variables: [windEastward, windNorthward] + + obs operator: + name: VertInterp + hofx scaling field: SurfaceWindScalingPressure + hofx scaling field group: DerivedVariables + + linear obs operator: + name: VertInterp + + # NOTE: Tests using the Gaussian Thinning filter (below) to duplicate GSI's thinning of AHI/Himawari-8 satwinds + # results in more JEDI satwinds in the diag file than in GSI, but far fewer JEDI satwinds assimilated than + # GSI. JEDI under-counts assimilated winds by roughly 25-40%, relative to GSI, and this under-count is not + # even including the temporal thinning which is applied in GSI but not JEDI (by this filter below). See + # GDASApp Issue #741 for details: https://github.com/NOAA-EMC/GDASApp/issues/741 + #obs pre filters: + #- filter: Gaussian Thinning + # horizontal_mesh: 200 + # vertical_mesh: 10000 + # use_reduced_horizontal_grid: true + # round_horizontal_bin_count_to_nearest: true + # partition_longitude_bins_using_mesh: true + + obs prior filters: + # Apply variable changes needed for wind scaling + # For wind observations with pressure provided + - filter: Variable Transforms + Transform: SurfaceWindScalingPressure + SkipWhenNoObs: False + + # Calculate error inflation factor for duplicate observations + #- filter: Variable Assignment + # assignments: + # - name: ObsErrorFactorDuplicateCheck/windEastward + # type: float + # function: + # name: ObsFunction/ObsErrorFactorDuplicateCheck + # options: + # use_air_pressure: true + # variable: windEastward + + #- filter: Variable Assignment + # assignments: + # - name: ObsErrorFactorDuplicateCheck/windNorthward + # type: float + # function: + # name: ObsFunction/ObsErrorFactorDuplicateCheck + # options: + # use_air_pressure: true + # variable: windNorthward + + obs post filters: + # Assign the initial observation error, based on height/pressure + # Hard-wiring to prepobs_errtable.global by Type + # ObsError is currently not updating in diag file, but passes directly to EffectiveError when no inflation is specified in YAML + + # Type 260 (VIIRS LWIR) + - filter: Perform Action + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 260 + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + xvar: + name: MetaData/pressure + xvals: [110000.,105000.,100000.,95000.,90000.,85000.,80000.,75000.,70000.,65000.,60000.,55000.,50000.,45000.,40000.,35000.,30000.,25000.,20000.,15000.,10000.,7500.,5000.,4000.,3000.,2000.,1000.,500.,400.,300.,200.,100.,0.] + errors: [3.8,3.8,3.8,3.8,3.8,3.8,3.8,3.8,3.9,3.9,4.,4.,4.1,5.,6.,6.3,6.6,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.] + + # sanity-check criteria + # Observation Range Sanity Check + # NOT EXPLICITLY CLEARED: No obs in this range in file, so 0 Bounds Check rejects (which is correct) but essentially untested + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + minvalue: -130. + maxvalue: 130. + action: + name: reject + + # Velocity Sanity Check + # NOT EXPLICITLY CLEARED: No obs in this range in file, so 0 Bounds Check rejects (which is correct) but essentially untested + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/Velocity + maxvalue: 130. + action: + name: reject + + # Reject any observation with a /=0 surface type (non-water surface) within + # 200 hPa of the surface pressure (as part of the LNVD check). + # CLEARED: maxvalue is rejecting >, not >= as per a Perform Action, so threshold is unchanged + - filter: Difference Check + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: + name: GeoVaLs/water_area_fraction + maxvalue: 0.99 + reference: GeoVaLs/surface_pressure + value: MetaData/pressure + maxvalue: -20000. # within 200 hPa above surface pressure, negative p-diff + action: + name: reject + + # LNVD check + # CLEARED: maxvalue is rejecting >, not >= as per a Perform Action, so threshold is unchanged + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/SatWindsLNVDCheck + maxvalue: 3. + action: + name: reject + + # GSI setupw routine QC + # Reject any ob Type [240–260] when pressure greater than 950 mb. + # CLEARED: minvalue/maxvalue are >=/<=, not >/<, so editing range by 1 Pa + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 240-260 + test variables: + - name: MetaData/pressure + maxvalue: 95001. + action: + name: reject + + # Multiple satellite platforms, reject when pressure is more than 50 mb above tropopause. + # CLEARED: minvalue is rejecting <, not <= as per a Perform Action, so threshold is unchanged + # Notes (eliu): This tropopause check reject too many obs; probably due to tropopause pressure estimation + # Turn this check off for now. + # Need to check if troposphere pressure was implemented correctly in fv3-jed + - filter: Difference Check + filter variables: + - name: windEastward + - name: windNorthward + reference: GeoVaLs/tropopause_pressure + value: MetaData/pressure + minvalue: -5000. # 50 hPa above tropopause level, negative p-diff + action: + name: reject + + # All satwinds must adjust errors based on ObsErrorFactorPressureCheck + # prior to the SPDB check (i.e. the gross-error check). The gross-error + # check uses the adjusted errors for error-bound tightening and rejection, + # so this check has to come first. This check will inflate errors for obs + # that are too close to either the model top or bottom. + # Notes (eliu): GMAO added a required parameter: adjusted_error_name. + - filter: Perform Action + filter variables: + - name: windEastward + where: + - variable: + name: ObsType/windEastward + is_in: 240-260 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck + options: + surface_obs: false + variable: windEastward + inflation factor: 4.0 + + - filter: Perform Action + filter variables: + - name: windNorthward + where: + - variable: + name: ObsType/windNorthward + is_in: 240-260 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck + options: + variable: windNorthward + inflation factor: 4.0 + + # All satwinds subject to a gross-error check that contains significant + # modifiers for satwinds with a negative speed-bias. ALL wind gross-error + # checks are currently being done by the SatWindsSPDBCheck. + # CLEARED + - filter: Background Check + filter variables: + - name: windEastward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260] + cgross: [ 2.5, 2.5, 2.5, 1.5, 2.5, 1.3, 1.3, 2.5, 2.5, 2.5, 2.5, 1.3, 2.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5] + error_min: [1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4] + error_max: [6.1, 6.1, 15.0, 15.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.1, 20.1, 20.1, 20.1, 20.1, 20.1] + variable: windEastward + action: + name: reject + + - filter: Background Check + filter variables: + - name: windNorthward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260] + cgross: [ 2.5, 2.5, 2.5, 1.5, 2.5, 1.3, 1.3, 2.5, 2.5, 2.5, 2.5, 1.3, 2.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5] + error_min: [1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4] + error_max: [6.1, 6.1, 15.0, 15.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.1, 20.1, 20.1, 20.1, 20.1, 20.1] + variable: windNorthward + action: + name: reject + + # The last error inflation check is for duplicate observations. This one needs + # to come last, because we don't want to inflate errors for duplication if one + # of the duplicates should be rejected. + # Notes (eliu): ObsErrorFactorDuplicateCheck obsfunction requires PreUseFlag (usage parameter from read_satwnd.f90). + # : Turn off duplicate check for now. + #- filter: Perform Action + # filter variables: + # - name: windEastward + # action: + # name: inflate error + # inflation variable: + # name: ObsErrorFactorDuplicateCheck/windEastward + + #- filter: Perform Action + # filter variables: + # - name: windNorthward + # action: + # name: inflate error + # inflation variable: + # name: ObsErrorFactorDuplicateCheck/windNorthward + + # We are extending this to an additional filter that inflates final ob-errors across-the-board by + # 1/0.8 = 1.25. This is caused by the GSI value of nvqc being set to .true. in the global operational + # configuration, see: https://github.com/NOAA-EMC/global-workflow/blob/d5ae3328fa4041b177357b1133f6b92e81c859d7/scripts/exglobal_atmos_analysis.sh#L750 + # This setting activates Line 1229 of setupw.f90 to scale ratio_errors by 0.8, which is applied in + # the denominator of the final ob-error, so 1/0.8 = 1.25 factor of ob-error inflation. + # + # If this nvqc functionality were to be switched off (i.e. if variational qc were to be turned off), + # you would want to remove this last inflation filter. + #- filter: Perform Action + # filter variables: + # - name: windEastward + # where: + # - variable: ObsType/windEastward + # is_in: 240-260 + # action: + # name: inflate error + # inflation factor: 1.25 + + #- filter: Perform Action + # filter variables: + # - name: windNorthward + # where: + # - variable: ObsType/windNorthward + # is_in: 240-260 + # action: + # name: inflate error + # inflation factor: 1.25 + + # End of Filters diff --git a/parm/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.json b/parm/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.json new file mode 100644 index 000000000..eab19c6b2 --- /dev/null +++ b/parm/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.json @@ -0,0 +1,16 @@ +{ + "data_format" : "bufr_d", + "data_type" : "satwnd", + "cycle_type" : "{{ RUN }}", + "cycle_datetime" : "{{ current_cycle | to_YMDH }}", + "dump_directory" : "{{ DMPDIR }}", + "ioda_directory" : "{{ COM_OBS }}", + "subsets" : [ "NC005091" ], + "data_description" : "NC005091 S-NPP/NOAA-20 SATWIND VIIRS IR(LW)", + "data_provider" : "U.S. NOAA/NESDIS", + "sensor_info" : { "sensor_name": "VIIRS", "sensor_full_name": "Visible Infrared Imaging Radiometer Suite", "sensor_id": 999 }, + "satellite_info" : [ + { "satellite_name": "npp", "satellite_full_name": "SUOMI-NPP", "satellite_id": 224, "launch time": "YYYYMMDD" }, + { "satellite_name": "n20", "satellite_full_name": "NOAA-20", "satellite_id": 225, "launch time": "YYYYMMDD" }, + ] +} diff --git a/ush/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.py b/ush/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.py new file mode 100755 index 000000000..18f87db04 --- /dev/null +++ b/ush/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.py @@ -0,0 +1,470 @@ +#!/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 + +# ==================================================================== +# Satellite Winds (AMV) BUFR dump file for VIIRS/S-NPP,NOAA-20 +# ==================================================================== +# Subset | Spectral Band | Code (002023) | ObsType +# -------------------------------------------------------------------- +# NC005091 | IRLW (Freq < 5E+13) | Method 1 | 260 +# ==================================================================== + +# Define and initialize global variables +global float32_fill_value +global int32_fill_value +global int64_fill_value + +float32_fill_value = np.float32(0) +int32_fill_value = np.int32(0) +int64_fill_value = np.int64(0) + + +def Compute_WindComponents_from_WindDirection_and_WindSpeed(wdir, wspd): + + uob = (-wspd * np.sin(np.radians(wdir))).astype(np.float32) + vob = (-wspd * np.cos(np.radians(wdir))).astype(np.float32) + + return uob, vob + + +def Get_ObsType(swcm, chanfreq): + + obstype = swcm.copy() + + # Use numpy vectorized operations + obstype = np.where(swcm == 1, 260, obstype) # IRLW + + if not np.any(np.isin(obstype, [260])): + raise ValueError("Error: Unassigned ObsType found ... ") + + return obstype + + +def bufr_to_ioda(config, logger): + + subsets = config["subsets"] + logger.debug(f"Checking subsets = {subsets}") + + # Get parameters from configuration + subsets = config["subsets"] + 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"] + 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 + converter = 'BUFR to IODA Converter' + process_level = 'Level-2' + platform_description = 'TERRA/AQUA' + sensor_description = 'Moderate Resolution Imaging Spectroradiometer' + + logger.info(f'sensor_name = {sensor_name}') + logger.info(f'sensor_full_name = {sensor_full_name}') + logger.info(f'sensor_id = {sensor_id}') + 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(subsets) + + # MetaData + q.add('latitude', '*/CLATH') + q.add('longitude', '*/CLONH') + 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('satelliteZenithAngle', '*/SAZA') + q.add('sensorCentralFrequency', '*/SCCF') + q.add('pressure', '*/PRLC[1]') + + # Processing Center + q.add('dataProviderOrigin', '*/OGCE[1]') + q.add('windGeneratingApplication', '*/AMVQIC/GNAPS') + +# # Quality Infomation (Quality Indicator w/o forecast) + q.add('qualityInformationWithoutForecast', '*/AMVQIC/PCCF') + + # Wind Retrieval Method Information + q.add('windComputationMethod', '*/SWCM') + + # ObsValue + q.add('windDirection', '*/WDIR') + q.add('windSpeed', '*/WSPD') + + 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') + satzenang = r.get('satelliteZenithAngle') + pressure = r.get('pressure', type='float') + chanfreq = r.get('sensorCentralFrequency', type='float') + + # Processing Center + ogce = r.get('dataProviderOrigin') + ga = r.get('windGeneratingApplication') + + # Quality Information + qi = r.get('qualityInformationWithoutForecast', type='float') + # For TERRA/AQUA MODIS data, qi w/o forecast (qifn) is packaged in same + # vector of qi with ga = 5 (QI without forecast), and EE is + # packaged in same vector of qi with ga=7 (Estimated Error (EE) in m/s + # converted to a percent confidence)shape (4,nobs). Must conduct a + # search and extract the correct vector for gnap and qi + # 1. Find dimension-sizes of ga and qi (should be the same!) + gDim1, gDim2 = np.shape(ga) + qDim1, qDim2 = np.shape(qi) + logger.info(f'Generating Application and Quality Information SEARCH:') + logger.info(f'Dimension size of GNAP ({gDim1},{gDim2})') + logger.info(f'Dimension size of PCCF ({qDim1},{qDim2})') + # 2. Initialize gnap and qifn as None, and search for dimension of + # ga with values of 5. If the same column exists for qi, assign + # gnap to ga[:,i] and qifn to qi[:,i], else raise warning that no + # appropriate GNAP/PCCF combination was found + gnap = None + qifn = None + for i in range(gDim2): + if np.unique(ga[:, i].squeeze()) == 5: + if i <= qDim2: + logger.info(f'GNAP/PCCF found for column {i}') + gnap = ga[:, i].squeeze() + qifn = qi[:, i].squeeze() + else: + logger.info(f'ERROR: GNAP column {i} outside of PCCF dimension {qDim2}') + if (gnap is None) & (qifn is None): + logger.info(f'ERROR: GNAP == 5 NOT FOUND OR OUT OF PCCF DIMENSION-RANGE, WILL FAIL!') + # If EE is needed, key search on np.unique(ga[:,i].squeeze()) == 7 instead + + # Wind Retrieval Method Information + swcm = r.get('windComputationMethod') + # ObsValue + # Wind direction and Speed + wdir = r.get('windDirection', type='float') + wspd = r.get('windSpeed') + + # 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) + + # Check BUFR variable generic dimension and type + + # Global variables declaration + # Set global fill values + float32_fill_value = satzenang.fill_value + int32_fill_value = satid.fill_value + int64_fill_value = timestamp.fill_value.astype(np.int64) + + end_time = time.time() + running_time = end_time - start_time + logger.info(f'Processing time for executing QuerySet to get ResultSet : {running_time} seconds') + + # ========================= + # Create derived variables + # ========================= + start_time = time.time() + + logger.info('Creating derived variables') + logger.debug('Creating derived variables - wind components (uob and vob)') + + uob, vob = Compute_WindComponents_from_WindDirection_and_WindSpeed(wdir, wspd) + + logger.debug(f' uob min/max = {uob.min()} {uob.max()}') + logger.debug(f' vob min/max = {vob.min()} {vob.max()}') + + obstype = Get_ObsType(swcm, chanfreq) + + height = np.full_like(pressure, fill_value=pressure.fill_value, dtype=np.float32) + stnelev = np.full_like(pressure, fill_value=pressure.fill_value, dtype=np.float32) + + end_time = time.time() + running_time = end_time - start_time + logger.info(f'Processing time for creating derived variables : {running_time} seconds') + + # ===================================== + # Split output based on satellite id + # Create IODA ObsSpace + # Write IODA output + # ===================================== + logger.info('Create IODA ObsSpace and Write IODA output based on satellite ID') + + # 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] + timestamp2 = timestamp[mask] + satid2 = satid[mask] + satzenang2 = satzenang[mask] + chanfreq2 = chanfreq[mask] + obstype2 = obstype[mask] + pressure2 = pressure[mask] + height2 = height[mask] + stnelev2 = stnelev[mask] + + # Processing Center + ogce2 = ogce[mask] + + # QC Info + qifn2 = qifn[mask] + + # Method + swcm2 = swcm[mask] + + # ObsValue + wdir2 = wdir[mask] + wspd2 = wspd[mask] + uob2 = uob[mask] + vob2 = vob[mask] + + # Timestamp Range + timestamp2_min = datetime.fromtimestamp(timestamp2.min()) + timestamp2_max = datetime.fromtimestamp(timestamp2.max()) + + # Check unique observation time + unique_timestamp2 = np.unique(timestamp2) + logger.debug(f'Processing output for satid {sat}') + + # Create the dimensions + dims = { + 'Location': np.arange(0, wdir2.shape[0]) + } + + # Create IODA ObsSpace + iodafile = f"{cycle_type}.t{hh}z.{data_type}.{satinst}.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('Write global attributes') + obsspace.write_attr('Converter', converter) + obsspace.write_attr('sourceFiles', bufrfile) + obsspace.write_attr('dataProviderOrigin', data_provider) + obsspace.write_attr('description', data_description) + obsspace.write_attr('datetimeReference', reference_time) + obsspace.write_attr('datetimeRange', [str(timestamp2_min), str(timestamp2_max)]) + obsspace.write_attr('sensor', sensor_id) + obsspace.write_attr('platform', satellite_id) + obsspace.write_attr('platformCommonName', satellite_name) + obsspace.write_attr('sensorCommonName', sensor_name) + obsspace.write_attr('processingLevel', process_level) + obsspace.write_attr('platformLongDescription', platform_description) + obsspace.write_attr('sensorLongDescription', sensor_description) + + # Create IODA variables + logger.debug('Write 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=lat.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=np.int64, fillval=int64_fill_value) \ + .write_attr('units', 'seconds since 1970-01-01T00:00:00Z') \ + .write_attr('long_name', 'Datetime') \ + .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) + + # Sensor Zenith Angle + obsspace.create_var('MetaData/satelliteZenithAngle', dtype=satzenang2.dtype, fillval=satzenang2.fill_value) \ + .write_attr('units', 'degree') \ + .write_attr('valid_range', np.array([0, 90], dtype=np.float32)) \ + .write_attr('long_name', 'Satellite Zenith Angle') \ + .write_data(satzenang2) + + # Sensor Centrall Frequency + obsspace.create_var('MetaData/sensorCentralFrequency', dtype=chanfreq2.dtype, fillval=chanfreq2.fill_value) \ + .write_attr('units', 'Hz') \ + .write_attr('long_name', 'Satellite Channel Center Frequency') \ + .write_data(chanfreq2) + + # Data Provider + obsspace.create_var('MetaData/dataProviderOrigin', dtype=ogce2.dtype, fillval=ogce2.fill_value) \ + .write_attr('long_name', 'Identification of Originating/Generating Center') \ + .write_data(ogce2) + + # Quality: Percent Confidence - Quality Information Without Forecast + obsspace.create_var('MetaData/qualityInformationWithoutForecast', dtype=qifn2.dtype, fillval=qifn2.fill_value) \ + .write_attr('long_name', 'Quality Information Without Forecast') \ + .write_data(qifn2) + + # Wind Computation Method + obsspace.create_var('MetaData/windComputationMethod', dtype=swcm2.dtype, fillval=swcm2.fill_value) \ + .write_attr('long_name', 'Satellite-derived Wind Computation Method') \ + .write_data(swcm2) + + # 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) + + # Height (mimic prepbufr) + obsspace.create_var('MetaData/height', dtype=height2.dtype, fillval=height2.fill_value) \ + .write_attr('units', 'm') \ + .write_attr('long_name', 'Height of Observation') \ + .write_data(height2) + + # Station Elevation (mimic prepbufr) + obsspace.create_var('MetaData/stationElevation', dtype=stnelev2.dtype, fillval=stnelev2.fill_value) \ + .write_attr('units', 'm') \ + .write_attr('long_name', 'Station Elevation') \ + .write_data(stnelev2) + + # ObsType based on computation method/spectral band + obsspace.create_var('ObsType/windEastward', dtype=obstype2.dtype, fillval=swcm2.fill_value) \ + .write_attr('long_name', 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band') \ + .write_data(obstype2) + + # ObsType based on computation method/spectral band + obsspace.create_var('ObsType/windNorthward', dtype=obstype2.dtype, fillval=swcm2.fill_value) \ + .write_attr('long_name', 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band') \ + .write_data(obstype2) + + # U-Wind Component + obsspace.create_var('ObsValue/windEastward', dtype=uob2.dtype, fillval=uob2.fill_value) \ + .write_attr('units', 'm s-1') \ + .write_attr('long_name', 'Eastward Wind Component') \ + .write_data(uob2) + + # V-Wind Component + obsspace.create_var('ObsValue/windNorthward', dtype=vob2.dtype, fillval=vob2.fill_value) \ + .write_attr('units', 'm s-1') \ + .write_attr('long_name', 'Northward Wind Component') \ + .write_data(vob2) + + end_time = time.time() + running_time = end_time - start_time + total_ob_processed += len(satid2) + logger.debug(f'Number of observation 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_satwind_amv_ahi.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 1b07517a22cd569d35ee24d341c15a97fc6ad932 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Thu, 25 Apr 2024 11:43:39 -0400 Subject: [PATCH 03/10] Addition of a switch for the cycling type (#1072) Draft because I still need to test this while cycling. - fixes #1071 --- scripts/exgdas_global_marine_analysis_post.py | 9 +++++++-- scripts/exgdas_global_marine_analysis_prep.py | 19 ++++++++++++++----- 2 files changed, 21 insertions(+), 7 deletions(-) diff --git a/scripts/exgdas_global_marine_analysis_post.py b/scripts/exgdas_global_marine_analysis_post.py index 54cce2363..1376bdd52 100755 --- a/scripts/exgdas_global_marine_analysis_post.py +++ b/scripts/exgdas_global_marine_analysis_post.py @@ -50,6 +50,7 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]): cyc = str(os.getenv('cyc')).zfill(2) bcyc = str((int(cyc) - 3) % 24).zfill(2) gcyc = str((int(cyc) - 6) % 24).zfill(2) # previous cycle +cdatedt = datetime.strptime(cdate, '%Y%m%d%H') bdatedt = datetime.strptime(cdate, '%Y%m%d%H') - timedelta(hours=3) bdate = datetime.strftime(bdatedt, '%Y-%m-%dT%H:00:00Z') mdate = datetime.strftime(datetime.strptime(cdate, '%Y%m%d%H'), '%Y-%m-%dT%H:00:00Z') @@ -93,8 +94,12 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]): os.path.join(com_ocean_analysis, f'{RUN}.t{bcyc}z.ocngrid.nc')]) # Copy the CICE analysis restart -cdateice = pdy + '.' + cyc + '0000' -post_file_list.append([os.path.join(anl_dir, 'Data', f'{cdateice}.cice_model.res.nc'), +if os.getenv('DOIAU') == "YES": + cice_rst_date = bdatedt.strftime('%Y%m%d.%H%M%S') +else: + cice_rst_date = cdatedt.strftime('%Y%m%d.%H%M%S') + +post_file_list.append([os.path.join(anl_dir, 'Data', f'{cice_rst_date}.cice_model.res.nc'), os.path.join(com_ice_restart, f'{cdate}.cice_model_anl.res.nc')]) FileHandler({'copy': post_file_list}).sync() diff --git a/scripts/exgdas_global_marine_analysis_prep.py b/scripts/exgdas_global_marine_analysis_prep.py index d80bfc834..686b4c8e8 100755 --- a/scripts/exgdas_global_marine_analysis_prep.py +++ b/scripts/exgdas_global_marine_analysis_prep.py @@ -70,10 +70,6 @@ def find_clim_ens(input_date): staticsoca_dir = os.getenv('SOCA_INPUT_FIX_DIR') nmem_ens = 0 nmem_ens = int(os.getenv('NMEM_ENS')) -if os.getenv('DOHYBVAR') == "YES": - dohybvar = True -else: - dohybvar = False # create analysis directories diags = os.path.join(anl_dir, 'diags') # output dir for soca DA obs space @@ -89,12 +85,25 @@ def find_clim_ens(input_date): window_begin = datetime.strptime(os.getenv('PDY')+os.getenv('cyc'), '%Y%m%d%H') - half_assim_freq window_begin_iso = window_begin.strftime('%Y-%m-%dT%H:%M:%SZ') window_middle_iso = window_middle.strftime('%Y-%m-%dT%H:%M:%SZ') -fcst_begin = datetime.strptime(os.getenv('PDY')+os.getenv('cyc'), '%Y%m%d%H') RUN = os.getenv('RUN') cyc = os.getenv('cyc') gcyc = os.getenv('gcyc') PDY = os.getenv('PDY') +# hybrid-envar switch +if os.getenv('DOHYBVAR') == "YES": + dohybvar = True +else: + dohybvar = False + +# switch for the cycling type +if os.getenv('DOIAU') == "YES": + # forecast initialized at the begining of the DA window + fcst_begin = window_begin +else: + # forecast initialized at the middle of the DA window + fcst_begin = datetime.strptime(os.getenv('PDY')+os.getenv('cyc'), '%Y%m%d%H') + ################################################################################ # fetch observations From a8daf6a1be634cabce5c923ff6aa442fd206d357 Mon Sep 17 00:00:00 2001 From: Cory Martin Date: Fri, 26 Apr 2024 19:05:16 +0000 Subject: [PATCH 04/10] Add modulefile for Dogwood/Cactus (#1073) While we are technically not supposed to build our own CRTM or FMS, we are asking for forgiveness and forging ahead. This PR should allow GDASApp to compile on WCOSS2 (Dogwood and Cactus) after a minor issue is resolved. That issue being a compile issue in IODA that Ron is already aware of and has a fix that works (https://github.com/JCSDA-internal/ioda/compare/develop...bugfix/read_script_convert_fix). --------- Co-authored-by: cory martin Co-authored-by: cory martin Co-authored-by: cory martin --- build.sh | 2 +- modulefiles/GDAS/wcoss2.intel.lua | 51 +++++++++++++++++++++++++++++++ sorc/ioda | 2 +- 3 files changed, 53 insertions(+), 2 deletions(-) create mode 100644 modulefiles/GDAS/wcoss2.intel.lua diff --git a/build.sh b/build.sh index a27da74e5..c81b47ae3 100755 --- a/build.sh +++ b/build.sh @@ -71,7 +71,7 @@ while getopts "p:t:c:hvdfa" opt; do done case ${BUILD_TARGET} in - hera | orion | hercules) + hera | orion | hercules | wcoss2) echo "Building GDASApp on $BUILD_TARGET" source $dir_root/ush/module-setup.sh module use $dir_root/modulefiles diff --git a/modulefiles/GDAS/wcoss2.intel.lua b/modulefiles/GDAS/wcoss2.intel.lua new file mode 100644 index 000000000..d21099b2e --- /dev/null +++ b/modulefiles/GDAS/wcoss2.intel.lua @@ -0,0 +1,51 @@ +help([[ +Load environment for running the GDAS application with Intel compilers and MPI. +]]) + +local pkgName = myModuleName() +local pkgVersion = myModuleVersion() +local pkgNameVer = myModuleFullName() + +prepend_path("MODULEPATH", "/apps/dev/lmodules/core") + +load("PrgEnv-intel/8.2.0") +load("cmake/3.20.2") +load("craype") +load("cray-pals") +load("git/2.29.0") +load("intel/19.1.3.304") +load("cray-mpich/8.1.12") +load("hdf5/1.12.2") +load("netcdf/4.7.4") +load("udunits/2.2.28") +load("eigen/3.4.0") +load("boost/1.79.0") +load("gsl-lite/v0.40.0") +load("sp/2.4.0") +load("python/3.8.6") +load("ecbuild/3.7.0") +load("qhull/2020.2") +load("eckit/1.24.4") +load("fckit/0.11.0") +load("atlas/0.35.0") +load("nccmp") + +-- hack for pybind11 +setenv("pybind11_ROOT", "/apps/spack/python/3.8.6/intel/19.1.3.304/pjn2nzkjvqgmjw4hmyz43v5x4jbxjzpk/lib/python3.8/site-packages/pybind11/share/cmake/pybind11") + +-- hack for wxflow +--prepend_path("PYTHONPATH", "/scratch1/NCEPDEV/da/python/gdasapp/wxflow/20240307/src") + +local mpiexec = '/pe/intel/compilers_and_libraries_2020.4.304/linux/mpi/intel64/bin/mpirun' +local mpinproc = '-n' +setenv('MPIEXEC_EXEC', mpiexec) +setenv('MPIEXEC_NPROC', mpinproc) + +setenv("CRTM_FIX","/lfs/h2/emc/da/noscrub/emc.da/GDASApp/fix/crtm/2.4.0") +setenv("GDASAPP_TESTDATA","/lfs/h2/emc/da/noscrub/emc.da/GDASApp/data") +setenv("GDASAPP_UNIT_TEST_DATA_PATH", "/lfs/h2/emc/da/noscrub/emc.da/GDASApp/data/test") + +whatis("Name: ".. pkgName) +whatis("Version: ".. pkgVersion) +whatis("Category: GDASApp") +whatis("Description: Load all libraries needed for GDASApp") diff --git a/sorc/ioda b/sorc/ioda index 91eb2c764..206eba708 160000 --- a/sorc/ioda +++ b/sorc/ioda @@ -1 +1 @@ -Subproject commit 91eb2c7643b33c1af2470d289eefb7aec2667387 +Subproject commit 206eba7084c08a3bc9a1c6ccbbff3e63d7cb602e From 7ecfc17f9df39859b340c8ba81c3c27bfb5f18ac Mon Sep 17 00:00:00 2001 From: Mindo Choi <141867620+apchoiCMD@users.noreply.github.com> Date: Fri, 26 Apr 2024 18:25:52 -0400 Subject: [PATCH 05/10] Using ioda util to convert the datetime in AMSR2 converter (#1077) #### This PR includes: - Replace the method to get datetime from epoch time by using `ioda utils` - `julianDate - 2440588` will make half day off - The above method will avoid dealing with decimal precision. - Closes #1062 --------- Co-authored-by: Guillaume Vernieres --- utils/obsproc/IcecAmsr2Ioda.h | 38 ++++++++++++++++++++++++----------- 1 file changed, 26 insertions(+), 12 deletions(-) diff --git a/utils/obsproc/IcecAmsr2Ioda.h b/utils/obsproc/IcecAmsr2Ioda.h index 681bbcd8d..3241d3ae2 100644 --- a/utils/obsproc/IcecAmsr2Ioda.h +++ b/utils/obsproc/IcecAmsr2Ioda.h @@ -13,8 +13,10 @@ #include // NOLINT +#include "ioda/../../../../core/IodaUtils.h" #include "ioda/Group.h" #include "ioda/ObsGroup.h" + #include "oops/util/dateFunctions.h" #include "NetCDFToIodaConverter.h" @@ -83,22 +85,34 @@ namespace gdasapp { int minute = oneTmpdateTimeVal[i+4]; int second = static_cast(oneTmpdateTimeVal[i+5]); - // Replace Fillvalue -9999 to 0 to avoid crash in dateToJulian + // Avoid crash util in ioda::convertDtimeToTimeOffsets if (year == -9999 || month == -9999 || day == -9999 || hour == -9999 || minute == -9999 || second == -9999) { year = month = day = hour = minute = second = 0; } - // Convert a date to Julian date - uint64_t julianDate = util::datefunctions::dateToJulian(year, month, day); - - // Subtract Julian day from January 1, 1970 (convert to epoch) - int daysSinceEpoch = julianDate - 2440588; - - // Calculate seconds only from HHMMSS - int secondsOffset = util::datefunctions::hmsToSeconds(hour, minute, second); - - iodaVars.datetime_(index) = static_cast(daysSinceEpoch*86400.0f) + secondsOffset; + // Construct iso8601 string format for each dateTime + std::stringstream ss; + ss << std::setfill('0') + << std::setw(4) << year << '-' + << std::setw(2) << month << '-' + << std::setw(2) << day << 'T' + << std::setw(2) << hour << ':' + << std::setw(2) << minute << ':' + << std::setw(2) << second << 'Z'; + std::string formattedDateTime = ss.str(); + util::DateTime dateTime(formattedDateTime); + + // Set epoch time for AMSR2_ICEC + util::DateTime epochDtime("1970-01-01T00:00:00Z"); + + // Convert Obs DateTime objects to epoch time offsets in seconds + // 0000-00-00T00:00:00Z will be converterd to negative seconds + int64_t timeOffsets + = ioda::convertDtimeToTimeOffsets(epochDtime, {dateTime})[0]; + + // Update datetime Eigen Arrays + iodaVars.datetime_(index) = timeOffsets; index++; } @@ -106,7 +120,7 @@ namespace gdasapp { for (int i = 0; i < iodaVars.location_; i++) { iodaVars.longitude_(i) = oneDimLonVal[i]; iodaVars.latitude_(i) = oneDimLatVal[i]; - iodaVars.obsVal_(i) = static_cast(oneDimObsVal[i]*0.01f); + iodaVars.obsVal_(i) = static_cast(oneDimObsVal[i]*0.01); iodaVars.obsError_(i) = 0.1; // Do something for obs error iodaVars.preQc_(i) = oneDimFlagsVal[i]; // Store optional metadata, set ocean basins to -999 for now From 69ef3fe4fefc05a6731dcafaefda210c219d4124 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Mon, 29 Apr 2024 15:48:13 -0400 Subject: [PATCH 06/10] The DA only uses the gdas bkg ... fixing again ... (#1079) - fixes #1078 --- scripts/exgdas_global_marine_analysis_post.py | 2 +- scripts/exgdas_global_marine_analysis_prep.py | 2 +- ush/soca/bkg_utils.py | 9 ++++----- ush/soca/marine_recenter.py | 2 +- 4 files changed, 7 insertions(+), 8 deletions(-) diff --git a/scripts/exgdas_global_marine_analysis_post.py b/scripts/exgdas_global_marine_analysis_post.py index 1376bdd52..d207840f0 100755 --- a/scripts/exgdas_global_marine_analysis_post.py +++ b/scripts/exgdas_global_marine_analysis_post.py @@ -176,7 +176,7 @@ def create_obs_space(data): # TODO(GorA): this should be setup properly in the g-w once gdassoca_obsstats is in develop gdassoca_obsstats_exec = os.path.join(os.getenv('HOMEgfs'), 'sorc', 'gdas.cd', 'build', 'bin', 'gdassoca_obsstats.x') -command = f"{os.getenv('launcher')} {gdassoca_obsstats_exec} {stats_yaml}" +command = f"{os.getenv('launcher')} -n 1 {gdassoca_obsstats_exec} {stats_yaml}" logger.info(f"{command}") result = subprocess.run(command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) diff --git a/scripts/exgdas_global_marine_analysis_prep.py b/scripts/exgdas_global_marine_analysis_prep.py index 686b4c8e8..265526eb0 100755 --- a/scripts/exgdas_global_marine_analysis_prep.py +++ b/scripts/exgdas_global_marine_analysis_prep.py @@ -336,7 +336,7 @@ def find_clim_ens(input_date): ################################################################################ # Copy initial condition -bkg_utils.stage_ic(bkg_dir, anl_dir, RUN, gcyc) +bkg_utils.stage_ic(bkg_dir, anl_dir, gcyc) ################################################################################ # prepare input.nml diff --git a/ush/soca/bkg_utils.py b/ush/soca/bkg_utils.py index 08be3dd89..488fcf697 100755 --- a/ush/soca/bkg_utils.py +++ b/ush/soca/bkg_utils.py @@ -101,13 +101,12 @@ def gen_bkg_list(bkg_path, out_path, window_begin=' ', yaml_name='bkg.yaml', ice bkg_date = window_begin # Construct list of background file names - RUN = os.getenv('RUN') cyc = str(os.getenv('cyc')).zfill(2) gcyc = str((int(cyc) - 6) % 24).zfill(2) # previous cycle fcst_hrs = list(range(3, 10, dt_pseudo)) files = [] for fcst_hr in fcst_hrs: - files.append(os.path.join(bkg_path, f'{RUN}.ocean.t'+gcyc+'z.inst.f'+str(fcst_hr).zfill(3)+'.nc')) + files.append(os.path.join(bkg_path, f'gdas.ocean.t'+gcyc+'z.inst.f'+str(fcst_hr).zfill(3)+'.nc')) # Identify the ocean background that will be used for the vertical coordinate remapping ocn_filename_ic = os.path.splitext(os.path.basename(files[0]))[0]+'.nc' @@ -159,17 +158,17 @@ def gen_bkg_list(bkg_path, out_path, window_begin=' ', yaml_name='bkg.yaml', ice FileHandler({'copy': bkg_list_src_dst}).sync() -def stage_ic(bkg_dir, anl_dir, RUN, gcyc): +def stage_ic(bkg_dir, anl_dir, gcyc): # Copy and rename initial condition ics_list = [] # ocean IC's - mom_ic_src = os.path.join(bkg_dir, f'{RUN}.ocean.t{gcyc}z.inst.f003.nc') + mom_ic_src = os.path.join(bkg_dir, f'gdas.ocean.t{gcyc}z.inst.f003.nc') mom_ic_dst = os.path.join(anl_dir, 'INPUT', 'MOM.res.nc') ics_list.append([mom_ic_src, mom_ic_dst]) # seaice IC's - cice_ic_src = os.path.join(bkg_dir, f'{RUN}.agg_ice.t{gcyc}z.inst.f003.nc') + cice_ic_src = os.path.join(bkg_dir, f'gdas.agg_ice.t{gcyc}z.inst.f003.nc') cice_ic_dst = os.path.join(anl_dir, 'INPUT', 'cice.res.nc') ics_list.append([cice_ic_src, cice_ic_dst]) FileHandler({'copy': ics_list}).sync() diff --git a/ush/soca/marine_recenter.py b/ush/soca/marine_recenter.py index 4c49bb530..5bb689e71 100644 --- a/ush/soca/marine_recenter.py +++ b/ush/soca/marine_recenter.py @@ -128,7 +128,7 @@ def initialize(self): ################################################################################ # Copy initial condition - bkg_utils.stage_ic(self.config.bkg_dir, self.runtime_config.DATA, RUN, gcyc) + bkg_utils.stage_ic(self.config.bkg_dir, self.runtime_config.DATA, gcyc) ################################################################################ # stage ensemble members From 8fdaf3d5229101364d38a6d86d75337871c3c374 Mon Sep 17 00:00:00 2001 From: Wei Huang Date: Tue, 30 Apr 2024 06:30:56 -0600 Subject: [PATCH 07/10] Add module files to compile on AWS (#1082) Add module files to compile on AWS Solve issue: #1081 https://github.com/NOAA-EMC/GDASApp/issues/1081 - Tested on AWS - Tested on Hera --- build.sh | 2 +- modulefiles/GDAS/noaacloud.intel.lua | 93 ++++++++++++++++++++++++++++ ush/module-setup.sh | 4 ++ 3 files changed, 98 insertions(+), 1 deletion(-) create mode 100644 modulefiles/GDAS/noaacloud.intel.lua diff --git a/build.sh b/build.sh index c81b47ae3..3a9e98fa8 100755 --- a/build.sh +++ b/build.sh @@ -71,7 +71,7 @@ while getopts "p:t:c:hvdfa" opt; do done case ${BUILD_TARGET} in - hera | orion | hercules | wcoss2) + hera | orion | hercules | wcoss2 | noaacloud) echo "Building GDASApp on $BUILD_TARGET" source $dir_root/ush/module-setup.sh module use $dir_root/modulefiles diff --git a/modulefiles/GDAS/noaacloud.intel.lua b/modulefiles/GDAS/noaacloud.intel.lua new file mode 100644 index 000000000..b976d75bc --- /dev/null +++ b/modulefiles/GDAS/noaacloud.intel.lua @@ -0,0 +1,93 @@ +help([[ +Load environment for running the GDAS application with Intel compilers and MPI. +]]) + +local pkgName = myModuleName() +local pkgVersion = myModuleVersion() +local pkgNameVer = myModuleFullName() + +prepend_path("MODULEPATH", "/contrib/spack-stack/spack-stack-1.6.0/envs/unified-env/install/modulefiles/Core") + +-- below two lines get us access to the spack-stack modules +load("stack-intel/2021.3.0") +load("stack-intel-oneapi-mpi/2021.3.0") +load("cmake/3.23.1") +load("gettext/0.19.8.1") +--load("libunistring/1.1") +--load("libidn2/2.3.4") +load("pcre2/10.42") +load("curl/8.4.0") +load("zlib/1.2.13") +load("git/1.8.3.1") +load("pkg-config/0.27.1") +load("hdf5/1.14.0") +load("parallel-netcdf/1.12.2") +load("netcdf-c/4.9.2") +load("nccmp/1.9.0.1") +load("netcdf-fortran/4.6.1") +load("nco/5.0.6") +load("parallelio/2.5.10") +load("wget/1.14") +load("boost/1.83.0") +load("bufr/12.0.1") +load("git-lfs/2.4.1") +load("ecbuild/3.7.2") +load("openjpeg/2.3.1") +load("eccodes/2.32.0") +load("eigen/3.4.0") +load("openblas/0.3.24") +load("eckit/1.24.5") +load("fftw/3.3.10") +load("fckit/0.11.0") +load("fiat/1.2.0") +load("ectrans/1.2.0") +load("atlas/0.35.1") +load("sp/2.5.0") +load("gsl-lite/0.37.0") +load("libjpeg/2.1.0") +load("krb5/1.15.1") +load("libtirpc/1.3.3") +load("hdf/4.2.15") +load("jedi-cmake/1.4.0") +load("libpng/1.6.37") +load("libxt/1.1.5") +load("libxmu/1.1.4") +load("libxpm/3.5.12") +load("libxaw/1.0.13") +load("udunits/2.2.28") +load("ncview/2.1.9") +load("netcdf-cxx4/4.3.1") +load("json/3.10.5") +load("crtm/2.4.0.1") +load("rocoto/1.3.3") +load("prod_util/2.1.1") + +load("py-jinja2/3.0.3") +load("py-netcdf4/1.5.8") +load("py-pybind11/2.11.0") +load("py-pycodestyle/2.11.0") +load("py-pyyaml/6.0") +load("py-scipy/1.11.3") +load("py-xarray/2023.7.0") +load("py-f90nml/1.4.3") +load("py-pip/23.1.2") + +-- hack for wxflow +--prepend_path("PYTHONPATH", "/contrib/Wei.Huang/data/glopara/data/gdasapp-wxflow-20240307/src/wxflow") + +setenv("CC","mpiicc") +setenv("FC","mpiifort") +setenv("CXX","mpiicpc") + +local mpiexec = '/apps/slurm/default/bin/srun' +local mpinproc = '-n' +setenv('MPIEXEC_EXEC', mpiexec) +setenv('MPIEXEC_NPROC', mpinproc) + +setenv("CRTM_FIX","/contrib/Wei.Huang/data/hack-orion/crtm/2.4.0_skylab_3.0") +setenv("GDASAPP_TESTDATA","/contrib/Wei.Huang/data/hack-orion/data") + +whatis("Name: ".. pkgName) +whatis("Version: ".. pkgVersion) +whatis("Category: GDASApp") +whatis("Description: Load all libraries needed for GDASApp") diff --git a/ush/module-setup.sh b/ush/module-setup.sh index ef7f73699..6a2aabf73 100755 --- a/ush/module-setup.sh +++ b/ush/module-setup.sh @@ -105,6 +105,10 @@ elif [[ $MACHINE_ID = discover* ]]; then export PATH=$PATH:$SPACK_ROOT/bin . $SPACK_ROOT/share/spack/setup-env.sh +elif [[ $MACHINE_ID = noaacloud* ]]; then + # We are on NOAA Cloud + module purge + else echo WARNING: UNKNOWN PLATFORM 1>&2 fi From 25e00c9dd5dbdafbe51a7bc67c75d2aa99e38b88 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Tue, 30 Apr 2024 12:36:31 -0400 Subject: [PATCH 08/10] Use the gdas bkg for the static B (#1084) - fixes #1083 --- parm/soca/berror/soca_diagb.yaml.j2 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/parm/soca/berror/soca_diagb.yaml.j2 b/parm/soca/berror/soca_diagb.yaml.j2 index 3e6b9d2ec..fcfd7442f 100644 --- a/parm/soca/berror/soca_diagb.yaml.j2 +++ b/parm/soca/berror/soca_diagb.yaml.j2 @@ -7,8 +7,8 @@ date: '{{ ATM_WINDOW_MIDDLE }}' background: date: '{{ ATM_WINDOW_MIDDLE }}' basename: ./bkg/ - ocn_filename: '{{ RUN }}.ocean.t{{ gcyc }}z.inst.f009.nc' - ice_filename: '{{ RUN }}.agg_ice.t{{ gcyc }}z.inst.f009.nc' + ocn_filename: 'gdas.ocean.t{{ gcyc }}z.inst.f009.nc' + ice_filename: 'gdas.agg_ice.t{{ gcyc }}z.inst.f009.nc' read_from_file: 1 background error: From 7436897ad63bbc1489c5bc1e09810728173c5e2b Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Wed, 1 May 2024 07:43:39 -0400 Subject: [PATCH 09/10] Time series of csv stats (#1086) Fixed a few issues in the time series plot. Does things like this: ![csv-ex](https://github.com/NOAA-EMC/GDASApp/assets/14031856/f33d91db-72c3-4e3d-a5fc-288ecbe51d36) --- ush/soca/gdassoca_obsstats.py | 51 ++++++++++++++++++++++++----------- 1 file changed, 35 insertions(+), 16 deletions(-) diff --git a/ush/soca/gdassoca_obsstats.py b/ush/soca/gdassoca_obsstats.py index 4741687bb..ac6f42185 100644 --- a/ush/soca/gdassoca_obsstats.py +++ b/ush/soca/gdassoca_obsstats.py @@ -10,6 +10,21 @@ import matplotlib.pyplot as plt import matplotlib.dates as mdates +colors = [ + "lightsteelblue", + "lightgreen", + "lightpink", + "lightsalmon", + "lightcoral", + "lightgoldenrodyellow", + "paleturquoise", + "palegreen", + "palegoldenrod", + "peachpuff", + "mistyrose", + "lavender" +] + class ObsStats: def __init__(self): @@ -19,14 +34,15 @@ def read_csv(self, filepaths): # Iterate through the list of file paths and append their data for filepath in filepaths: new_data = pd.read_csv(filepath) + # Convert date to datetime for easier plotting new_data['date'] = pd.to_datetime(new_data['date'], format='%Y%m%d%H') self.data = pd.concat([self.data, new_data], ignore_index=True) + self.data.sort_values('date', inplace=True) - def plot_timeseries(self, ocean, variable): + def plot_timeseries(self, ocean, variable, inst="", dirout=""): # Filter data for the given ocean and variable filtered_data = self.data[(self.data['Ocean'] == ocean) & (self.data['Variable'] == variable)] - if filtered_data.empty: print("No data available for the given ocean and variable combination.") return @@ -36,13 +52,16 @@ def plot_timeseries(self, ocean, variable): # Plot settings fig, axs = plt.subplots(3, 1, figsize=(10, 15), sharex=True) - fig.suptitle(f'{variable} statistics, {ocean} ocean', fontsize=16, fontweight='bold') + fig.suptitle(f'{inst} {variable} statistics, {ocean} ocean', fontsize=16, fontweight='bold') + exp_counter = 0 for exp in experiments: - exp_data = filtered_data[filtered_data['Exp'] == exp] + exp_data = self.data[(self.data['Ocean'] == ocean) & + (self.data['Variable'] == variable) & + (self.data['Exp'] == exp)] # Plot RMSE - axs[0].plot(exp_data['date'], exp_data['RMSE'], marker='o', linestyle='-', label=exp) + axs[0].plot(exp_data['date'], exp_data['RMSE'], marker='o', linestyle='-', color=colors[exp_counter], label=exp) axs[0].xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H')) axs[0].xaxis.set_major_locator(mdates.DayLocator()) axs[0].tick_params(labelbottom=False) @@ -51,42 +70,42 @@ def plot_timeseries(self, ocean, variable): axs[0].grid(True) # Plot Bias - axs[1].plot(exp_data['date'], exp_data['Bias'], marker='o', linestyle='-', label=exp) + axs[1].plot(exp_data['date'], exp_data['Bias'], marker='o', linestyle='-', color=colors[exp_counter], label=exp) axs[1].xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H')) axs[1].xaxis.set_major_locator(mdates.DayLocator()) axs[1].tick_params(labelbottom=False) axs[1].set_ylabel('Bias') - axs[1].legend() axs[1].grid(True) # Plot Count - axs[2].plot(exp_data['date'], exp_data['Count'], marker='o', linestyle='-', label=exp) + axs[2].plot(exp_data['date'], exp_data['Count'], marker='o', linestyle='-', color=colors[exp_counter], label=exp) axs[2].xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H')) axs[2].xaxis.set_major_locator(mdates.DayLocator()) axs[2].set_ylabel('Count') - axs[2].legend() axs[2].grid(True) # Improve layout and show plot plt.tight_layout(rect=[0, 0.03, 1, 0.95]) - plt.savefig(f'{ocean}_test.png') + plt.savefig(f'{dirout}/{inst}_{variable}_{ocean}_test.png') if __name__ == "__main__": - epilog = ["Usage examples: ./gdassoca_obsstats.py --exp gdas.*.csv"] + epilog = ["Usage examples: ./gdassoca_obsstats.py --exps cp1/COMROOT/cp1 cp2/COMROOT/cp2 --inst sst_abi_g16_l3c --dirout cp1vscp2"] parser = argparse.ArgumentParser(description="Observation space RMSE's and BIAS's", formatter_class=argparse.RawDescriptionHelpFormatter, epilog=os.linesep.join(epilog)) parser.add_argument("--exps", nargs='+', required=True, help="Path to the experiment's COMROOT") - parser.add_argument("--variable", required=True, help="ombg_qc or ombg_noqc") - parser.add_argument("--figname", required=True, help="The name of the output figure") + parser.add_argument("--inst", required=True, help="The name of the instrument/platform (ex: sst_abi_g16_l3c)") + parser.add_argument("--dirout", required=True, help="Output directory") args = parser.parse_args() flist = [] + inst = args.inst + os.makedirs(args.dirout, exist_ok=True) + for exp in args.exps: - print(exp) - wc = exp + '/*.*/??/analysis/ocean/*.t??z.stats.csv' + wc = exp + f'/*.*/??/analysis/ocean/*{inst}*.stats.csv' flist.append(glob.glob(wc)) flist = sum(flist, []) @@ -94,4 +113,4 @@ def plot_timeseries(self, ocean, variable): obsStats.read_csv(flist) for var, ocean in product(['ombg_noqc', 'ombg_qc'], ['Atlantic', 'Pacific', 'Indian', 'Arctic', 'Southern']): - obsStats.plot_timeseries(ocean, var) + obsStats.plot_timeseries(ocean, var, inst=inst, dirout=args.dirout) From 706e11a777b3fd63566ce00bbcaeff12f6135c11 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Wed, 1 May 2024 19:38:54 -0400 Subject: [PATCH 10/10] Visualize stats in simple html document (#1089) I got tiered of looking for figures and added a basic html file that displays what we need. Ugly and basic, but it does the trick. The page looks like this: ![html-stats](https://github.com/NOAA-EMC/GDASApp/assets/14031856/e6608697-c004-4d6c-82a0-d3169856d591) _ps: I'd like the same thing for the 100's of figure that used to come out of the verify job._ --- .../soca/fig_gallery}/gdassoca_obsstats.py | 3 +- utils/soca/fig_gallery/time_series_omb.html | 120 ++++++++++++++++++ 2 files changed, 121 insertions(+), 2 deletions(-) rename {ush/soca => utils/soca/fig_gallery}/gdassoca_obsstats.py (98%) create mode 100644 utils/soca/fig_gallery/time_series_omb.html diff --git a/ush/soca/gdassoca_obsstats.py b/utils/soca/fig_gallery/gdassoca_obsstats.py similarity index 98% rename from ush/soca/gdassoca_obsstats.py rename to utils/soca/fig_gallery/gdassoca_obsstats.py index ac6f42185..3f92f799d 100644 --- a/ush/soca/gdassoca_obsstats.py +++ b/utils/soca/fig_gallery/gdassoca_obsstats.py @@ -86,8 +86,7 @@ def plot_timeseries(self, ocean, variable, inst="", dirout=""): # Improve layout and show plot plt.tight_layout(rect=[0, 0.03, 1, 0.95]) - plt.savefig(f'{dirout}/{inst}_{variable}_{ocean}_test.png') - + plt.savefig(f'{dirout}/{inst}_{variable}_{ocean}.png') if __name__ == "__main__": epilog = ["Usage examples: ./gdassoca_obsstats.py --exps cp1/COMROOT/cp1 cp2/COMROOT/cp2 --inst sst_abi_g16_l3c --dirout cp1vscp2"] diff --git a/utils/soca/fig_gallery/time_series_omb.html b/utils/soca/fig_gallery/time_series_omb.html new file mode 100644 index 000000000..4e06b5e56 --- /dev/null +++ b/utils/soca/fig_gallery/time_series_omb.html @@ -0,0 +1,120 @@ + + + + + + + + +
+ + + +
+
+ + + + + + +
+
+ + +
+ +
+ + + +