diff --git a/parm/atm/obs/config/conv_ps.yaml b/parm/atm/obs/config/conv_ps.yaml new file mode 100644 index 000000000..8fa4ab185 --- /dev/null +++ b/parm/atm/obs/config/conv_ps.yaml @@ -0,0 +1,295 @@ +obs space: + name: surface_ps + obsdatain: + engine: + type: H5File + obsfile: $(DATA)/obs/$(OPREFIX)sondes.{{ current_cycle | to_YMDH }}.nc4 + obsdataout: + engine: + type: H5File + obsfile: $(DATA)/diags/diag_sondes_{{ current_cycle | to_YMDH }}.nc4 + io pool: + max pool size: 1 + simulated variables: [stationPressure] + +obs operator: + name: SfcPCorrected + variables: + - name: stationPressure + da_psfc_scheme: GSI + station_altitude: height + geovar_sfc_geomz: surface_altitude + geovar_geomz: geopotential_height + +obs prior filters: +# Initial Error Assignments for SFC Observations +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: ObsType/stationPressure + is_in: [181] + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + round_to_the_nearest_integer: true + xvar: + name: ObsValue/stationPressure + xvals: [80000, 75000, 70000, 65000, 60000, 55000 ] + errors: [110, 120, 120, 120, 120, 1.0e+11] + +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: ObsType/stationPressure + is_in: [187] + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + round_to_the_nearest_integer: true + xvar: + name: ObsValue/stationPressure + xvals: [85000, 80000, 75000, 70000, 65000, 60000, 55000 ] + errors: [ 120, 140, 140, 140, 140, 140, 1.0e+11] + +# Initial Error Assignments for SFCSHIP Observations +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: ObsType/stationPressure + is_in: [180] + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + round_to_the_nearest_integer: true + xvar: + name: ObsValue/stationPressure + xvals: [60000, 55000 ] + errors: [ 130, 1.0e+11] + +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: ObsType/stationPressure + is_in: [183] + action: + name: assign error + error parameter: 1.0e+11 + +# Initial Error Assignments for Radiosonde +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: ObsType/stationPressure + is_in: [120] + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + round_to_the_nearest_integer: true + xvar: + name: ObsValue/stationPressure + xvals: [80000, 75000, 70000, 65000, 60000, 55000 ] + errors: [ 110, 120, 120, 120, 120, 1.0e+11] + +obs post filters: +# Observation range sanity check +- filter: Bounds Check + filter variables: + - name: stationPressure + minvalue: 37499.0 + maxvalue: 106999.0 + action: + name: reject + +# Reject all ObsType 183 +- filter: RejectList + where: + - variable: + name: ObsType/stationPressure + is_in: 183 + +# Reject surface pressure below 500 hPa +- filter: Bounds Check + filter variables: + - name: stationPressure + minvalue: 50000.00 + action: + name: reject + +- filter: RejectList + where: + - variable: + name: PreQC/stationPressure + is_in: 4-15 + +# Inflate obs error based on obs type +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: PreQC/stationPressure + is_in: 3, 7 + action: + name: inflate error + inflation factor: 1.2 + +# Calculate obs error inflation factors for duplicated observations at the same location +- filter: Variable Assignment + assignments: + - name: ObsErrorFactorDuplicateCheck/stationPressure + type: float + function: + name: ObsFunction/ObsErrorFactorDuplicateCheck + options: + use_air_pressure: false + variable: stationPressure + +# Reduce effective observation error based on obs type and subtype +# In this case: reduce effective obs error for buoy +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: + name: ObsType/stationPressure + is_in: 180 + - variable: + name: ObsSubType/stationPressure + is_in: 0 + action: + name: inflate error + inflation factor: 0.7 + +# Reduce original observation error based on obs type and subtype +# In this case: reduce original obs error for buoy +- filter: Variable Assignment + where: + - variable: + name: ObsType/stationPressure + is_in: 180 + - variable: + name: ObsSubType/stationPressure + is_in: 0 + assignments: + - name: ObsError/stationPressure + type: float + function: + name: ObsFunction/Arithmetic + options: + variables: + - name: ObsError/stationPressure + coefs: [0.7] + +# Calculate obs error inflation factors for large discrepancies between model and observations +- filter: Variable Assignment + assignments: + - name: ObsErrorFactorSfcPressure/stationPressure + type: float + function: + name: ObsFunction/ObsErrorFactorSfcPressure + options: + geovar_sfc_geomz: surface_altitude + geovar_geomz: geopotential_height + station_altitude: height + +# Inflate surface pressure observation based on discrepancies between +# model and observations due to terrian +- filter: Perform Action + filter variables: + - name: stationPressure + action: + name: inflate error + inflation variable: + name: ObsErrorFactorSfcPressure/stationPressure + +- filter: Variable Assignment + assignments: + - name: DerivedMetaData/Innovation + type: float + function: + name: ObsFunction/Arithmetic + options: + variables: + - name: ObsValue/stationPressure + - name: HofX/stationPressure + coefs: [1, -1] + +- filter: Variable Assignment + assignments: + - name: DerivedMetaData/ObsErrorBoundSfcPressure1 + type: float + function: + name: ObsFunction/ObsErrorBoundConventional + options: + obsvar: stationPressure + obserr_bound_min: 100 + obserr_bound_max: 300 + obserr_bound_factor: 5.0 + +- filter: Background Check + filter variables: + - name: stationPressure + where: + - variable: PreQC/stationPressure + is_not_in: 3 + function absolute threshold: + - name: DerivedMetaData/ObsErrorBoundSfcPressure1 + action: + name: reject + +- filter: Variable Assignment + assignments: + - name: DerivedMetaData/ObsErrorBoundSfcPressure2 + type: float + function: + name: ObsFunction/ObsErrorBoundConventional + options: + obsvar: stationPressure + obserr_bound_min: 100 + obserr_bound_max: 300 + obserr_bound_factor: 3.5 + +- filter: Background Check + filter variables: + - name: stationPressure + where: + - variable: PreQC/stationPressure + is_in: 3 + function absolute threshold: + - name: DerivedMetaData/ObsErrorBoundSfcPressure2 + action: + name: reject + +# Inflate obs error based on duplicate check +- filter: Perform Action + filter variables: + - name: stationPressure + action: + name: inflate error + inflation variable: + name: ObsErrorFactorDuplicateCheck/stationPressure + +# Reject data based on PreUseFlag (usage in GSI) +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: PreUseFlag/stationPressure + is_not_in: 0, 1 + action: + name: reject + +# End of Filters + diff --git a/parm/atm/obs/config/scatwind_metop-a.yaml b/parm/atm/obs/config/scatwind_metop-a.yaml new file mode 100644 index 000000000..4735c31fe --- /dev/null +++ b/parm/atm/obs/config/scatwind_metop-a.yaml @@ -0,0 +1,280 @@ +obs space: + name: ascatw_ascat_metop-a + obsdatain: + engine: + type: H5File + obsfile: $(DATA)/obs/$(OPREFIX)ascatw.ascat_metop-a.{{ current_cycle | to_YMDH }}.nc4 + obsdataout: + engine: + type: H5File + obsfile: $(DATA)/diags/diag_ascatw_ascat_metop-a_{{ current_cycle | to_YMDH }}.nc4 + io pool: + max pool size: 1 + simulated variables: [windEastward, windNorthward] + +obs operator: + name: VertInterp + # Use height vertical coordinate first +# vertical coordinate: geometric_height + vertical coordinate: geopotential_height + observation vertical coordinate group: DerivedVariables + observation vertical coordinate: adjustedHeight + interpolation method: linear + hofx scaling field: SurfaceWindScalingHeight + hofx scaling field group: DerivedVariables + +obs prior filters: +# Apply variable changes needed for rescaled height coordinate +- filter: Variable Transforms + Transform: AdjustedHeightCoordinate + SkipWhenNoObs: False + +# Apply variable changes needed for wind scaling +- filter: Variable Transforms + Transform: SurfaceWindScalingHeight + SkipWhenNoObs: False + +# Assign the initial observation error (constant value, 1.5 m/s right now). +- filter: Perform Action + filter variables: + - name: windEastward + - name: windNorthward + action: + name: assign error + error parameter: 1.5 + +# 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 + +# Reject all obs with PreQC mark already set above 3 +# NOTE: All scatwinds have an automatic PreQC mark of 2 (hard-wired default from GSI) +# - filter: PreQC +# maxvalue: 3 +# action: +# name: reject + +obs post filters: +# Reject all ASCAT (Type 290) winds with tsavg <= 273.0 (surface temperature) +- filter: Perform Action + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 290 +# - variable: GeoVaLs/surface_temperature + - variable: GeoVaLs/surface_temperature_where_land + maxvalue: 273. + action: + name: reject + +# Reject all ASCAT (Type 290) winds with isflg /= 0 (non-water surface) +- filter: Perform Action + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 290 + - variable: GeoVaLs/water_area_fraction + maxvalue: 0.99 + action: + name: reject + +# Reject ASCAT (Type 290) when observed component deviates from background by more than 5.0 m/s +# NOTE: This check can reject a u- or v-component of the same observation independently, which +# is fundamentally different from how GSI rejects obs (both components are rejected if +# either component fails a check). +- filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/Arithmetic + options: + variables: + - name: ObsValue/windEastward + - name: HofX/windEastward + coefs: [1.0, -1.0] + minvalue: -5.0 + maxvalue: 5.0 + +- filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/Arithmetic + options: + variables: + - name: ObsValue/windNorthward + - name: HofX/windNorthward + coefs: [1.0, -1.0] + minvalue: -5.0 + maxvalue: 5.0 + +# Reject OSCAT (Type 291) when observed component deviates from background by more than 6.0 m/s +# NOTE: This check can reject a u- or v-component of the same observation independently, which +# is fundamentally different from how GSI rejects obs (both components are rejected if +# either component fails a check). +- filter: Background Check + filter variables: + - name: windEastward + - name: windNorthward + threshold: 6. + absolute threshold: 6. + where: + - variable: ObsType/windEastward + is_in: 291 + action: + name: reject + +# Reject ASCAT (Type 290) when ambiguity check fails (returned value is negative) +- filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 290 + test variables: + - name: ObsFunction/ScatWindsAmbiguityCheck + options: + minimum_uv: 0.0001 # hard-coding a minimum-uv for transparancy, want this to basically be zero + maxvalue: 0. + action: + name: reject + +# All scatwinds must adjust errors based on ObsErrorFactorPressureCheck. +# This check will inflate errors for obs that are too close to either +# the model top or bottom. +- filter: Perform Action + filter variables: + - name: windEastward + where: + - variable: + name: ObsType/windEastward + is_in: 290-291 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck + options: + variable: windEastward + inflation factor: 4.0 + +- filter: Perform Action + filter variables: + - name: windNorthward + where: + - variable: + name: ObsType/windNorthward + is_in: 290-291 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck + options: + variable: windNorthward + inflation factor: 4.0 + +# All scatwinds subject to a gross error check. This is contained within +# the WindsSPDBCheck, although it is not exclusive to satwinds. +- filter: Background Check + filter variables: + - name: windEastward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 290, 291] + cgross: [ 5.0, 5.0] + error_min: [1.4, 1.4] + error_max: [6.1, 6.1] + variable: windEastward + action: + name: reject + +- filter: Background Check + filter variables: + - name: windNorthward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 290, 291] + cgross: [ 5.0, 5.0] + error_min: [1.4, 1.4] + error_max: [6.1, 6.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. +#- 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 + +# There is no across-the-board inflation for nvqc=.true. for scatwinds, presumably because for +# this inflation to take place both nvqc must be .true. AND ibeta must be >0, see: +# https://github.com/NOAA-EMC/GSI/blob/14ae595af1b03471287d322596d35c0665336e95/src/gsi/setupw.f90#L1229 +# GSI settings must have ibeta>0 for satwinds, but not for scatwinds. +# +# If the ibeta settings for scatwinds were to change while nvqc remained .true., we would extend YAML to +# an additional filter that inflates final ob-errors across-the-board by 1/0.8 = 1.25. NOTE: the nvqc setting +# is defaulted to .false. in GSI code, but is overridden in global operational configuration. See: +# 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 functionality were to be activated for scatwinds, you would want to include this last inflation filter. +#- filter: Perform Action +# filter variables: +# - name: windEastward +# where: +# - variable: ObsType/windEastward +# is_in: 290-291 +# action: +# name: inflate error +# inflation factor: 1.25 +#- filter: Perform Action +# filter variables: +# - name: windNorthward +# where: +# - variable: ObsType/windNorthward +# is_in: 290-291 +# action: +# name: inflate error +# inflation factor: 1.25 + +# END OF FILTERS# diff --git a/parm/atm/obs/config/scatwind_metop-b.yaml b/parm/atm/obs/config/scatwind_metop-b.yaml new file mode 100644 index 000000000..5ed7c38bc --- /dev/null +++ b/parm/atm/obs/config/scatwind_metop-b.yaml @@ -0,0 +1,280 @@ +obs space: + name: ascatw_ascat_metop-b + obsdatain: + engine: + type: H5File + obsfile: $(DATA)/obs/$(OPREFIX)ascatw.ascat_metop-b.{{ current_cycle | to_YMDH }}.nc4 + obsdataout: + engine: + type: H5File + obsfile: $(DATA)/diags/diag_ascatw_ascat_metop-b_{{ current_cycle | to_YMDH }}.nc4 + io pool: + max pool size: 1 + simulated variables: [windEastward, windNorthward] + +obs operator: + name: VertInterp + # Use height vertical coordinate first +# vertical coordinate: geometric_height + vertical coordinate: geopotential_height + observation vertical coordinate group: DerivedVariables + observation vertical coordinate: adjustedHeight + interpolation method: linear + hofx scaling field: SurfaceWindScalingHeight + hofx scaling field group: DerivedVariables + +obs prior filters: +# Apply variable changes needed for rescaled height coordinate +- filter: Variable Transforms + Transform: AdjustedHeightCoordinate + SkipWhenNoObs: False + +# Apply variable changes needed for wind scaling +- filter: Variable Transforms + Transform: SurfaceWindScalingHeight + SkipWhenNoObs: False + +# Assign the initial observation error (constant value, 1.5 m/s right now). +- filter: Perform Action + filter variables: + - name: windEastward + - name: windNorthward + action: + name: assign error + error parameter: 1.5 + +# 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 + +# Reject all obs with PreQC mark already set above 3 +# NOTE: All scatwinds have an automatic PreQC mark of 2 (hard-wired default from GSI) +# - filter: PreQC +# maxvalue: 3 +# action: +# name: reject + +obs post filters: +# Reject all ASCAT (Type 290) winds with tsavg <= 273.0 (surface temperature) +- filter: Perform Action + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 290 +# - variable: GeoVaLs/surface_temperature + - variable: GeoVaLs/surface_temperature_where_land + maxvalue: 273. + action: + name: reject + +# Reject all ASCAT (Type 290) winds with isflg /= 0 (non-water surface) +- filter: Perform Action + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 290 + - variable: GeoVaLs/water_area_fraction + maxvalue: 0.99 + action: + name: reject + +# Reject ASCAT (Type 290) when observed component deviates from background by more than 5.0 m/s +# NOTE: This check can reject a u- or v-component of the same observation independently, which +# is fundamentally different from how GSI rejects obs (both components are rejected if +# either component fails a check). +- filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/Arithmetic + options: + variables: + - name: ObsValue/windEastward + - name: HofX/windEastward + coefs: [1.0, -1.0] + minvalue: -5.0 + maxvalue: 5.0 + +- filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/Arithmetic + options: + variables: + - name: ObsValue/windNorthward + - name: HofX/windNorthward + coefs: [1.0, -1.0] + minvalue: -5.0 + maxvalue: 5.0 + +# Reject OSCAT (Type 291) when observed component deviates from background by more than 6.0 m/s +# NOTE: This check can reject a u- or v-component of the same observation independently, which +# is fundamentally different from how GSI rejects obs (both components are rejected if +# either component fails a check). +- filter: Background Check + filter variables: + - name: windEastward + - name: windNorthward + threshold: 6. + absolute threshold: 6. + where: + - variable: ObsType/windEastward + is_in: 291 + action: + name: reject + +# Reject ASCAT (Type 290) when ambiguity check fails (returned value is negative) +- filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 290 + test variables: + - name: ObsFunction/ScatWindsAmbiguityCheck + options: + minimum_uv: 0.0001 # hard-coding a minimum-uv for transparancy, want this to basically be zero + maxvalue: 0. + action: + name: reject + +# All scatwinds must adjust errors based on ObsErrorFactorPressureCheck. +# This check will inflate errors for obs that are too close to either +# the model top or bottom. +- filter: Perform Action + filter variables: + - name: windEastward + where: + - variable: + name: ObsType/windEastward + is_in: 290-291 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck + options: + variable: windEastward + inflation factor: 4.0 + +- filter: Perform Action + filter variables: + - name: windNorthward + where: + - variable: + name: ObsType/windNorthward + is_in: 290-291 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck + options: + variable: windNorthward + inflation factor: 4.0 + +# All scatwinds subject to a gross error check. This is contained within +# the WindsSPDBCheck, although it is not exclusive to satwinds. +- filter: Background Check + filter variables: + - name: windEastward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 290, 291] + cgross: [ 5.0, 5.0] + error_min: [1.4, 1.4] + error_max: [6.1, 6.1] + variable: windEastward + action: + name: reject + +- filter: Background Check + filter variables: + - name: windNorthward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 290, 291] + cgross: [ 5.0, 5.0] + error_min: [1.4, 1.4] + error_max: [6.1, 6.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. +#- 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 + +# There is no across-the-board inflation for nvqc=.true. for scatwinds, presumably because for +# this inflation to take place both nvqc must be .true. AND ibeta must be >0, see: +# https://github.com/NOAA-EMC/GSI/blob/14ae595af1b03471287d322596d35c0665336e95/src/gsi/setupw.f90#L1229 +# GSI settings must have ibeta>0 for satwinds, but not for scatwinds. +# +# If the ibeta settings for scatwinds were to change while nvqc remained .true., we would extend YAML to +# an additional filter that inflates final ob-errors across-the-board by 1/0.8 = 1.25. NOTE: the nvqc setting +# is defaulted to .false. in GSI code, but is overridden in global operational configuration. See: +# 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 functionality were to be activated for scatwinds, you would want to include this last inflation filter. +#- filter: Perform Action +# filter variables: +# - name: windEastward +# where: +# - variable: ObsType/windEastward +# is_in: 290-291 +# action: +# name: inflate error +# inflation factor: 1.25 +#- filter: Perform Action +# filter variables: +# - name: windNorthward +# where: +# - variable: ObsType/windNorthward +# is_in: 290-291 +# action: +# name: inflate error +# inflation factor: 1.25 + +# END OF FILTERS# diff --git a/parm/atm/obs/testing/conv_ps.yaml b/parm/atm/obs/testing/conv_ps.yaml index f43bc9f41..20a449a8a 100644 --- a/parm/atm/obs/testing/conv_ps.yaml +++ b/parm/atm/obs/testing/conv_ps.yaml @@ -8,8 +8,8 @@ obs space: engine: type: H5File obsfile: !ENV conv_ps_diag_${CDATE}.nc4 - overwrite: true simulated variables: [stationPressure] + geovals: filename: !ENV conv_ps_geoval_${CDATE}.nc4 @@ -22,191 +22,276 @@ obs operator: geovar_sfc_geomz: surface_altitude geovar_geomz: geopotential_height +obs prior filters: +# Initial Error Assignments for SFC Observations +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: ObsType/stationPressure + is_in: [181] + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + round_to_the_nearest_integer: true + xvar: + name: ObsValue/stationPressure + xvals: [80000, 75000, 70000, 65000, 60000, 55000 ] + errors: [110, 120, 120, 120, 120, 1.0e+11] + +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: ObsType/stationPressure + is_in: [187] + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + round_to_the_nearest_integer: true + xvar: + name: ObsValue/stationPressure + xvals: [85000, 80000, 75000, 70000, 65000, 60000, 55000 ] + errors: [ 120, 140, 140, 140, 140, 140, 1.0e+11] + +# Initial Error Assignments for SFCSHIP Observations +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: ObsType/stationPressure + is_in: [180] + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + round_to_the_nearest_integer: true + xvar: + name: ObsValue/stationPressure + xvals: [60000, 55000 ] + errors: [ 130, 1.0e+11] + +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: ObsType/stationPressure + is_in: [183] + action: + name: assign error + error parameter: 1.0e+11 + +# Initial Error Assignments for Radiosonde +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: ObsType/stationPressure + is_in: [120] + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + round_to_the_nearest_integer: true + xvar: + name: ObsValue/stationPressure + xvals: [80000, 75000, 70000, 65000, 60000, 55000 ] + errors: [ 110, 120, 120, 120, 120, 1.0e+11] + obs post filters: - # Observation range sanity check - - filter: Bounds Check - filter variables: - - name: stationPressure - minvalue: 37499.0 - maxvalue: 106999.0 - action: - name: reject - - # Reject all ObsType 183 - - filter: RejectList - where: - - variable: - name: ObsType/stationPressure - is_in: 183 - - # Reject surface pressure below 500 hPa - - filter: Bounds Check - filter variables: - - name: stationPressure - minvalue: 50000.00 - action: - name: reject - - - filter: RejectList - where: - - variable: - name: PreQC/stationPressure - is_in: 4-15 - - # Inflate obs error based on obs type - - filter: Perform Action - filter variables: - - name: stationPressure - where: - - variable: PreQC/stationPressure - is_in: 3, 7 - action: - name: inflate error - inflation factor: 1.2 - - # Calculate obs error inflation factors for duplicated observations at the same location - - filter: Variable Assignment - assignments: - - name: ObsErrorFactorDuplicateCheck/stationPressure - type: float - function: - name: ObsFunction/ObsErrorFactorDuplicateCheck - options: - use_air_pressure: false - variable: stationPressure - - # Reduce effective observation error based on obs type and subtype - # In this case: reduce effective obs error for buoy - - filter: Perform Action - filter variables: - - name: stationPressure - where: - - variable: - name: ObsType/stationPressure - is_in: 180 - - variable: - name: ObsSubType/stationPressure - is_in: 0 - action: - name: inflate error - inflation factor: 0.7 - - # Reduce original observation error based on obs type and subtype - # In this case: reduce original obs error for buoy - - filter: Variable Assignment - where: - - variable: - name: ObsType/stationPressure - is_in: 180 - - variable: - name: ObsSubType/stationPressure - is_in: 0 - assignments: - - name: ObsError/stationPressure - type: float - function: - name: ObsFunction/Arithmetic - options: - variables: - - name: ObsError/stationPressure - coefs: [0.7] - - # Calculate obs error inflation factors for large discrepancies between model and observations - - filter: Variable Assignment - assignments: - - name: ObsErrorFactorSfcPressure/stationPressure - type: float - function: - name: ObsFunction/ObsErrorFactorSfcPressure - options: - geovar_sfc_geomz: surface_altitude - - # Inflate surface pressure observation based on discrepancies between - # model and observations due to terrian - - filter: Perform Action - filter variables: - - name: stationPressure - action: - name: inflate error - inflation variable: - name: ObsErrorFactorSfcPressure/stationPressure - - - filter: Variable Assignment - assignments: - - name: DerivedMetaData/Innovation - type: float - function: - name: ObsFunction/Arithmetic - options: - variables: - - name: ObsValue/stationPressure - - name: HofX/stationPressure - coefs: [1, -1] - - - filter: Variable Assignment - assignments: - - name: DerivedMetaData/ObsErrorBoundSfcPressure1 - type: float - function: - name: ObsFunction/ObsErrorBoundConventional - options: - obsvar: stationPressure - obserr_bound_min: 100 - obserr_bound_max: 300 - obserr_bound_factor: 5.0 - - - filter: Background Check - filter variables: - - name: stationPressure - where: - - variable: PreQC/stationPressure - is_not_in: 3 - function absolute threshold: - - name: DerivedMetaData/ObsErrorBoundSfcPressure1 - action: - name: reject - - - filter: Variable Assignment - assignments: - - name: DerivedMetaData/ObsErrorBoundSfcPressure2 - type: float - function: - name: ObsFunction/ObsErrorBoundConventional - options: - obsvar: stationPressure - obserr_bound_min: 100 - obserr_bound_max: 300 - obserr_bound_factor: 3.5 - - - filter: Background Check - filter variables: - - name: stationPressure - where: - - variable: PreQC/stationPressure - is_in: 3 - function absolute threshold: - - name: DerivedMetaData/ObsErrorBoundSfcPressure2 - action: - name: reject - - # Inflate obs error based on duplicate check - - filter: Perform Action - filter variables: - - name: stationPressure - action: - name: inflate error - inflation variable: - name: ObsErrorFactorDuplicateCheck/stationPressure - - # Reject data based on PreUseFlag (usage in GSI) - - filter: Perform Action - filter variables: - - name: stationPressure - where: - - variable: PreUseFlag/stationPressure - is_not_in: 0, 1 - action: - name: reject +# Observation range sanity check +- filter: Bounds Check + filter variables: + - name: stationPressure + minvalue: 37499.0 + maxvalue: 106999.0 + action: + name: reject + +# Reject all ObsType 183 +- filter: RejectList + where: + - variable: + name: ObsType/stationPressure + is_in: 183 + +# Reject surface pressure below 500 hPa +- filter: Bounds Check + filter variables: + - name: stationPressure + minvalue: 50000.00 + action: + name: reject + +- filter: RejectList + where: + - variable: + name: PreQC/stationPressure + is_in: 4-15 + +# Inflate obs error based on obs type +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: PreQC/stationPressure + is_in: 3, 7 + action: + name: inflate error + inflation factor: 1.2 + +# Calculate obs error inflation factors for duplicated observations at the same location +- filter: Variable Assignment + assignments: + - name: ObsErrorFactorDuplicateCheck/stationPressure + type: float + function: + name: ObsFunction/ObsErrorFactorDuplicateCheck + options: + use_air_pressure: false + variable: stationPressure + +# Reduce effective observation error based on obs type and subtype +# In this case: reduce effective obs error for buoy +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: + name: ObsType/stationPressure + is_in: 180 + - variable: + name: ObsSubType/stationPressure + is_in: 0 + action: + name: inflate error + inflation factor: 0.7 + +# Reduce original observation error based on obs type and subtype +# In this case: reduce original obs error for buoy +- filter: Variable Assignment + where: + - variable: + name: ObsType/stationPressure + is_in: 180 + - variable: + name: ObsSubType/stationPressure + is_in: 0 + assignments: + - name: ObsError/stationPressure + type: float + function: + name: ObsFunction/Arithmetic + options: + variables: + - name: ObsError/stationPressure + coefs: [0.7] + +# Calculate obs error inflation factors for large discrepancies between model and observations +- filter: Variable Assignment + assignments: + - name: ObsErrorFactorSfcPressure/stationPressure + type: float + function: + name: ObsFunction/ObsErrorFactorSfcPressure + options: + geovar_sfc_geomz: surface_altitude + geovar_geomz: geopotential_height + station_altitude: height + +# Inflate surface pressure observation based on discrepancies between +# model and observations due to terrian +- filter: Perform Action + filter variables: + - name: stationPressure + action: + name: inflate error + inflation variable: + name: ObsErrorFactorSfcPressure/stationPressure + +- filter: Variable Assignment + assignments: + - name: DerivedMetaData/Innovation + type: float + function: + name: ObsFunction/Arithmetic + options: + variables: + - name: ObsValue/stationPressure + - name: HofX/stationPressure + coefs: [1, -1] + +- filter: Variable Assignment + assignments: + - name: DerivedMetaData/ObsErrorBoundSfcPressure1 + type: float + function: + name: ObsFunction/ObsErrorBoundConventional + options: + obsvar: stationPressure + obserr_bound_min: 100 + obserr_bound_max: 300 + obserr_bound_factor: 5.0 + +- filter: Background Check + filter variables: + - name: stationPressure + where: + - variable: PreQC/stationPressure + is_not_in: 3 + function absolute threshold: + - name: DerivedMetaData/ObsErrorBoundSfcPressure1 + action: + name: reject + +- filter: Variable Assignment + assignments: + - name: DerivedMetaData/ObsErrorBoundSfcPressure2 + type: float + function: + name: ObsFunction/ObsErrorBoundConventional + options: + obsvar: stationPressure + obserr_bound_min: 100 + obserr_bound_max: 300 + obserr_bound_factor: 3.5 + +- filter: Background Check + filter variables: + - name: stationPressure + where: + - variable: PreQC/stationPressure + is_in: 3 + function absolute threshold: + - name: DerivedMetaData/ObsErrorBoundSfcPressure2 + action: + name: reject + +# Inflate obs error based on duplicate check +- filter: Perform Action + filter variables: + - name: stationPressure + action: + name: inflate error + inflation variable: + name: ObsErrorFactorDuplicateCheck/stationPressure + +# Reject data based on PreUseFlag (usage in GSI) +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: PreUseFlag/stationPressure + is_not_in: 0, 1 + action: + name: reject +# End of Filters passedBenchmark: 85378 diff --git a/parm/atm/obs/testing/conv_ps_noqc.yaml b/parm/atm/obs/testing/conv_ps_noqc.yaml index 44a6ef4c8..8d469729f 100644 --- a/parm/atm/obs/testing/conv_ps_noqc.yaml +++ b/parm/atm/obs/testing/conv_ps_noqc.yaml @@ -8,8 +8,8 @@ obs space: engine: type: H5File obsfile: !ENV conv_ps_diag_${CDATE}.nc4 - overwrite: true simulated variables: [stationPressure] + geovals: filename: !ENV conv_ps_geoval_${CDATE}.nc4 @@ -22,10 +22,89 @@ obs operator: geovar_sfc_geomz: surface_altitude geovar_geomz: geopotential_height -vector ref: GsiHofXBc -tolerance: 1.e-4 -#linear obs operator test: -# coef TL: 0.1 -# tolerance TL: 1.0e-13 -# tolerance AD: 1.0e-11 +obs prior filters: +# Initial Error Assignments for SFC Observations +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: ObsType/stationPressure + is_in: [181] + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + round_to_the_nearest_integer: true + xvar: + name: ObsValue/stationPressure + xvals: [80000, 75000, 70000, 65000, 60000, 55000 ] + errors: [110, 120, 120, 120, 120, 1.0e+11] + +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: ObsType/stationPressure + is_in: [187] + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + round_to_the_nearest_integer: true + xvar: + name: ObsValue/stationPressure + xvals: [85000, 80000, 75000, 70000, 65000, 60000, 55000 ] + errors: [ 120, 140, 140, 140, 140, 140, 1.0e+11] + +# Initial Error Assignments for SFCSHIP Observations +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: ObsType/stationPressure + is_in: [180] + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + round_to_the_nearest_integer: true + xvar: + name: ObsValue/stationPressure + xvals: [60000, 55000 ] + errors: [ 130, 1.0e+11] + +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: ObsType/stationPressure + is_in: [183] + action: + name: assign error + error parameter: 1.0e+11 + + # Initial Error Assignments for Radiosonde +- filter: Perform Action + filter variables: + - name: stationPressure + where: + - variable: ObsType/stationPressure + is_in: [120] + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + round_to_the_nearest_integer: true + xvar: + name: ObsValue/stationPressure + xvals: [80000, 75000, 70000, 65000, 60000, 55000 ] + errors: [ 110, 120, 120, 120, 120, 1.0e+11] + +passedBenchmark: 92824 # total: 92842; missing: 18 +#vector ref: GsiHofXBc +#tolerance: 1.e-4 diff --git a/parm/atm/obs/testing/scatwind.yaml b/parm/atm/obs/testing/scatwind.yaml index 30a1af6f7..9befdda81 100644 --- a/parm/atm/obs/testing/scatwind.yaml +++ b/parm/atm/obs/testing/scatwind.yaml @@ -9,228 +9,283 @@ obs space: type: H5File obsfile: !ENV scatwind_diag_${CDATE}.nc4 simulated variables: [windEastward, windNorthward] + geovals: filename: !ENV scatwind_geoval_${CDATE}.nc4 -vector ref: GsiHofXBc -tolerance: 0.01 + obs operator: name: VertInterp - apply near surface wind scaling: true -# -obs pre filters: - # Assign the initial observation error (constant value, 1.5 m/s right now). - - filter: Perform Action - filter variables: - - name: windEastward - - name: windNorthward - action: - name: assign error - error parameter: 1.5 -#obs prior filters: - # - # Reject all obs with PreQC mark already set above 3 - # NOTE: All scatwinds have an automatic PreQC mark of 2 (hard-wired default from GSI) - # - filter: PreQC - # maxvalue: 3 - # action: - # name: reject - # + # Use height vertical coordinate first +# vertical coordinate: geometric_height + vertical coordinate: geopotential_height + observation vertical coordinate group: DerivedVariables + observation vertical coordinate: adjustedHeight + interpolation method: linear + hofx scaling field: SurfaceWindScalingHeight + hofx scaling field group: DerivedVariables + +obs prior filters: +# Apply variable changes needed for rescaled height coordinate +- filter: Variable Transforms + Transform: AdjustedHeightCoordinate + SkipWhenNoObs: False + +# Apply variable changes needed for wind scaling +- filter: Variable Transforms + Transform: SurfaceWindScalingHeight + SkipWhenNoObs: False + +# Assign the initial observation error (constant value, 1.5 m/s right now). +- filter: Perform Action + filter variables: + - name: windEastward + - name: windNorthward + action: + name: assign error + error parameter: 1.5 + +# 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 + +# Reject all obs with PreQC mark already set above 3 +# NOTE: All scatwinds have an automatic PreQC mark of 2 (hard-wired default from GSI) +# - filter: PreQC +# maxvalue: 3 +# action: +# name: reject + obs post filters: - # - # Reject all ASCAT (Type 290) winds with tsavg <= 273.0 (surface temperature) - # - - filter: Perform Action - filter variables: - - name: windEastward - - name: windNorthward - where: - - variable: ObsType/windEastward - is_in: 290 - - variable: GeoVaLs/surface_temperature - maxvalue: 273. - action: - name: reject - # - # Reject all ASCAT (Type 290) winds with isflg /= 0 (non-water surface) - # - - filter: Perform Action - filter variables: - - name: windEastward - - name: windNorthward - where: - - variable: ObsType/windEastward - is_in: 290 - - variable: GeoVaLs/water_area_fraction - maxvalue: 0.001 - action: - name: reject - # Reject ASCAT (Type 290) when observed component deviates from background by more than 5.0 m/s - # NOTE: This check can reject a u- or v-component of the same observation independently, which - # is fundamentally different from how GSI rejects obs (both components are rejected if - # either component fails a check). - - filter: Background Check - filter variables: - - name: windEastward - - name: windNorthward - threshold: 5. - absolute threshold: 5. - where: - - variable: ObsType/windEastward - is_in: 290 - action: - name: reject - # Reject OSCAT (Type 291) when observed component deviates from background by more than 6.0 m/s - # NOTE: This check can reject a u- or v-component of the same observation independently, which - # is fundamentally different from how GSI rejects obs (both components are rejected if - # either component fails a check). - - filter: Background Check - filter variables: - - name: windEastward - - name: windNorthward - threshold: 6. - absolute threshold: 6. - where: - - variable: ObsType/windEastward - is_in: 291 - action: - name: reject - # Reject ASCAT (Type 290) when ambiguity check fails (returned value is negative) - - filter: Bounds Check - filter variables: - - name: windEastward - - name: windNorthward - where: - - variable: ObsType/windEastward - is_in: 290 - test variables: - - name: ObsFunction/ScatWindsAmbiguityCheck - options: - test_hofx: GsiHofX - minimum_uv: 0.0001 # hard-coding a minimum-uv for transparancy, want this to basically be zero - maxvalue: 0. - action: - name: reject -# - # All scatwinds must adjust errors based on ObsErrorFactorPressureCheck. - # This check will inflate errors for obs that are too close to either - # the model top or bottom. - - filter: Perform Action - filter variables: - - name: windEastward - where: - - variable: - name: ObsType/windEastward - is_in: 290-291 - action: - name: inflate error - inflation variable: - name: ObsFunction/ObsErrorFactorPressureCheck - options: - variable: windEastward - inflation factor: 4.0 -# - - filter: Perform Action - filter variables: - - name: windNorthward - where: - - variable: - name: ObsType/windNorthward - is_in: 290-291 - action: - name: inflate error - inflation variable: - name: ObsFunction/ObsErrorFactorPressureCheck - options: - variable: windNorthward - inflation factor: 4.0 - # All scatwinds subject to a gross error check. This is contained within - # the WindsSPDBCheck, although it is not exclusive to satwinds. - - filter: Background Check - filter variables: - - name: windEastward - function absolute threshold: - - name: ObsFunction/WindsSPDBCheck +# Reject all ASCAT (Type 290) winds with tsavg <= 273.0 (surface temperature) +- filter: Perform Action + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 290 +# - variable: GeoVaLs/surface_temperature + - variable: GeoVaLs/surface_temperature_where_land + maxvalue: 273. + action: + name: reject + +# Reject all ASCAT (Type 290) winds with isflg /= 0 (non-water surface) +- filter: Perform Action + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 290 + - variable: GeoVaLs/water_area_fraction + maxvalue: 0.99 +# - variable: GeoVaLs/surface_type +# minvalue: 0.001 + action: + name: reject +#passedBenchmark: 53558 # 2 variables (u,v), u=26779 passing, v=26779 passing + +# Reject ASCAT (Type 290) when observed component deviates from background by more than 5.0 m/s +# NOTE: This check can reject a u- or v-component of the same observation independently, which +# is fundamentally different from how GSI rejects obs (both components are rejected if +# either component fails a check). +- filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/Arithmetic + options: + variables: + - name: ObsValue/windEastward + - name: HofX/windEastward + coefs: [1.0, -1.0] + minvalue: -5.0 + maxvalue: 5.0 + +- filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/Arithmetic + options: + variables: + - name: ObsValue/windNorthward + - name: HofX/windNorthward + coefs: [1.0, -1.0] + minvalue: -5.0 + maxvalue: 5.0 +#passedBenchmark: 52362 # 2 variables (u,v), u=26181 passing, v=26181 passing + +# Reject OSCAT (Type 291) when observed component deviates from background by more than 6.0 m/s +# NOTE: This check can reject a u- or v-component of the same observation independently, which +# is fundamentally different from how GSI rejects obs (both components are rejected if +# either component fails a check). +- filter: Background Check + filter variables: + - name: windEastward + - name: windNorthward + threshold: 6. + absolute threshold: 6. + where: + - variable: ObsType/windEastward + is_in: 291 + action: + name: reject + +# Reject ASCAT (Type 290) when ambiguity check fails (returned value is negative) +- filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 290 + test variables: + - name: ObsFunction/ScatWindsAmbiguityCheck + options: +# test_hofx: GsiHofX + minimum_uv: 0.0001 # hard-coding a minimum-uv for transparancy, want this to basically be zero + maxvalue: 0. + action: + name: reject +#passedBenchmark: 51776 # 2 variables (u,v), u=25888 passing, v=25888 passing + +# All scatwinds must adjust errors based on ObsErrorFactorPressureCheck. +# This check will inflate errors for obs that are too close to either +# the model top or bottom. +- filter: Perform Action + filter variables: + - name: windEastward + where: + - variable: + name: ObsType/windEastward + is_in: 290-291 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck options: - wndtype: [ 290, 291] - cgross: [ 5.0, 5.0] - error_min: [1.4, 1.4] - error_max: [6.1, 6.1] variable: windEastward - action: - name: reject - - - filter: Background Check - filter variables: - - name: windNorthward - function absolute threshold: - - name: ObsFunction/WindsSPDBCheck + inflation factor: 4.0 + +- filter: Perform Action + filter variables: + - name: windNorthward + where: + - variable: + name: ObsType/windNorthward + is_in: 290-291 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck options: - wndtype: [ 290, 291] - cgross: [ 5.0, 5.0] - error_min: [1.4, 1.4] - error_max: [6.1, 6.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. - - filter: RejectList - filter variables: - - name: windEastward - action: - name: inflate error - inflation variable: - name: ObsFunction/ObsErrorFactorDuplicateCheck - options: - use_air_pressure: true - variable: windEastward - - filter: RejectList - filter variables: - - name: windNorthward - action: - name: inflate error - inflation variable: - name: ObsFunction/ObsErrorFactorDuplicateCheck - options: - use_air_pressure: true - variable: windNorthward - # There is no across-the-board inflation for nvqc=.true. for scatwinds, presumably because for - # this inflation to take place both nvqc must be .true. AND ibeta must be >0, see: - # https://github.com/NOAA-EMC/GSI/blob/14ae595af1b03471287d322596d35c0665336e95/src/gsi/setupw.f90#L1229 - # GSI settings must have ibeta>0 for satwinds, but not for scatwinds. - # - # - # If the ibeta settings for scatwinds were to change while nvqc remained .true., we would extend YAML to - # an additional filter that inflates final ob-errors across-the-board by 1/0.8 = 1.25. NOTE: the nvqc setting - # is defaulted to .false. in GSI code, but is overridden in global operational configuration. See: - # 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 functionality were to be activated for scatwinds, you would want to include this last inflation filter. -# - filter: Perform Action -# filter variables: -# - name: windEastward -# where: -# - variable: ObsType/windEastward -# is_in: 290-291 -# action: -# name: inflate error -# inflation factor: 1.25 -# - filter: Perform Action -# filter variables: -# - name: windNorthward -# where: -# - variable: ObsType/windNorthward -# is_in: 290-291 -# action: -# name: inflate error -# inflation factor: 1.25 -linear obs operator: - name: Identity -#passedBenchmark: 52097 # 2 variables (u,v), u=26087 passing, v=26010 passing -passedBenchmark: 52109 # 2 variables (u,v), u=26093 passing, v=26016 passing + inflation factor: 4.0 +#passedBenchmark: 51776 # 2 variables (u,v), u=25888 passing, v=25888 passing + +# All scatwinds subject to a gross error check. This is contained within +# the WindsSPDBCheck, although it is not exclusive to satwinds. +- filter: Background Check + filter variables: + - name: windEastward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 290, 291] + cgross: [ 5.0, 5.0] + error_min: [1.4, 1.4] + error_max: [6.1, 6.1] + variable: windEastward + action: + name: reject + +- filter: Background Check + filter variables: + - name: windNorthward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 290, 291] + cgross: [ 5.0, 5.0] + error_min: [1.4, 1.4] + error_max: [6.1, 6.1] + variable: windNorthward + action: + name: reject +#passedBenchmark: 51776 # 2 variables (u,v), u=25888 passing, v=25888 passing + +# 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. +#- 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 + +# There is no across-the-board inflation for nvqc=.true. for scatwinds, presumably because for +# this inflation to take place both nvqc must be .true. AND ibeta must be >0, see: +# https://github.com/NOAA-EMC/GSI/blob/14ae595af1b03471287d322596d35c0665336e95/src/gsi/setupw.f90#L1229 +# GSI settings must have ibeta>0 for satwinds, but not for scatwinds. +# +# If the ibeta settings for scatwinds were to change while nvqc remained .true., we would extend YAML to +# an additional filter that inflates final ob-errors across-the-board by 1/0.8 = 1.25. NOTE: the nvqc setting +# is defaulted to .false. in GSI code, but is overridden in global operational configuration. See: +# 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 functionality were to be activated for scatwinds, you would want to include this last inflation filter. +#- filter: Perform Action +# filter variables: +# - name: windEastward +# where: +# - variable: ObsType/windEastward +# is_in: 290-291 +# action: +# name: inflate error +# inflation factor: 1.25 +#- filter: Perform Action +# filter variables: +# - name: windNorthward +# where: +# - variable: ObsType/windNorthward +# is_in: 290-291 +# action: +# name: inflate error +# inflation factor: 1.25 +passedBenchmark: 51764 # 2 variables (u,v), u=25882 passing, v=25882 passing # GSI rejects both u- and v-component in first-guess check, UFO does not, but when UFO rej is reconciled btwn u and v these match GSI rej # u: 207 obs pass with corresponding v being rejected, 25880 obs in UFO/GSI agreement # v: 130 obs pass with corresponding u being rejected, 25880 obs in UFO/GSI agreement diff --git a/parm/atm/obs/testing/scatwind_noqc.yaml b/parm/atm/obs/testing/scatwind_noqc.yaml index f8f36efeb..ec4aef452 100644 --- a/parm/atm/obs/testing/scatwind_noqc.yaml +++ b/parm/atm/obs/testing/scatwind_noqc.yaml @@ -9,13 +9,42 @@ obs space: type: H5File obsfile: !ENV scatwind_diag_${CDATE}.nc4 simulated variables: [windEastward, windNorthward] + geovals: filename: !ENV scatwind_geoval_${CDATE}.nc4 -vector ref: GsiHofXBc -tolerance: 0.01 + obs operator: name: VertInterp - apply near surface wind scaling: true + # Use height vertical coordinate first +# vertical coordinate: geometric_height + vertical coordinate: geopotential_height + observation vertical coordinate group: DerivedVariables + observation vertical coordinate: adjustedHeight + interpolation method: linear + hofx scaling field: SurfaceWindScalingHeight + hofx scaling field group: DerivedVariables + +obs prior filters: +# Apply variable changes needed for rescaled height coordinate +- filter: Variable Transforms + Transform: AdjustedHeightCoordinate + SkipWhenNoObs: False + +# Apply variable changes needed for wind scaling +- filter: Variable Transforms + Transform: SurfaceWindScalingHeight + SkipWhenNoObs: False + +# Assign the initial observation error (constant value, 1.5 m/s) +- filter: Perform Action + filter variables: + - name: windEastward + - name: windNorthward + action: + name: assign error + error parameter: 1.5 -vector ref: GsiHofXBc -tolerance: 1.e-5 +passedBenchmark: 53546 +#passedBenchmark: 53558 +#vector ref: GsiHofXBc +#tolerance: 1.e-5 diff --git a/parm/ioda/bufr2ioda/bufr2ioda_adpsfc_prepbufr.json b/parm/ioda/bufr2ioda/bufr2ioda_adpsfc_prepbufr.json new file mode 100644 index 000000000..bdde19135 --- /dev/null +++ b/parm/ioda/bufr2ioda/bufr2ioda_adpsfc_prepbufr.json @@ -0,0 +1,12 @@ +{ + "data_format" : "prepbufr", + "subsets" : [ "ADPSFC" ], + "source" : "prepBUFR", + "data_type" : "ADPSFC", + "cycle_type" : "{{ RUN }}", + "cycle_datetime" : "{{ current_cycle | to_YMDH }}", + "dump_directory" : "{{ DMPDIR }}", + "ioda_directory" : "{{ COM_OBS }}", + "data_description" : "ADPSFC_prepbufr", + "data_provider" : "U.S. NOAA" +} diff --git a/parm/ioda/bufr2ioda/bufr2ioda_sfcshp_prepbufr.json b/parm/ioda/bufr2ioda/bufr2ioda_sfcshp_prepbufr.json new file mode 100644 index 000000000..4bdc8f781 --- /dev/null +++ b/parm/ioda/bufr2ioda/bufr2ioda_sfcshp_prepbufr.json @@ -0,0 +1,12 @@ +{ + "data_format" : "prepbufr", + "subsets" : [ "SFCSHP" ], + "source" : "prepBUFR", + "data_type" : "SFCSHP", + "cycle_type" : "{{ RUN }}", + "cycle_datetime" : "{{ current_cycle | to_YMDH }}", + "dump_directory" : "{{ DMPDIR }}", + "ioda_directory" : "{{ COM_OBS }}", + "data_description" : "SFCSHP_prepbufr", + "data_provider" : "U.S. NOAA" +} diff --git a/parm/soca/berror/soca_ensb.yaml b/parm/soca/berror/soca_ensb.yaml index 814dfc7ac..447d5b5fd 100644 --- a/parm/soca/berror/soca_ensb.yaml +++ b/parm/soca/berror/soca_ensb.yaml @@ -30,6 +30,12 @@ steric height: linear variable changes: - linear variable change name: BalanceSOCA # Only the steric balance is applied +ensemble mean output: + datadir: ./static_ens + date: '{{ATM_WINDOW_BEGIN}}' + exp: ens_mean + type: incr + ssh output: unbalanced: datadir: ./static_ens diff --git a/parm/soca/berror/soca_ensweights.yaml b/parm/soca/berror/soca_ensweights.yaml index c7ea40a98..094377f95 100644 --- a/parm/soca/berror/soca_ensweights.yaml +++ b/parm/soca/berror/soca_ensweights.yaml @@ -19,6 +19,16 @@ weights: # Need to provide weights^2 when reading from file ice: 0.0025 # 5% of original variance ocean: 0.0625 # 25% " " + # Apply localized weights to the ocean ens. B + ocean local weights: + - lon: -172.0 + lat: 11.0 + amplitude: -1.0 + length scale: 700.0 + - lon: -160.0 + lat: 12.0 + amplitude: -1.0 + length scale: 700.0 output: datadir: ./ diff --git a/parm/soca/obs/config/salt_profile_argo.yaml b/parm/soca/obs/config/salt_profile_argo.yaml new file mode 100644 index 000000000..9f3e2521a --- /dev/null +++ b/parm/soca/obs/config/salt_profile_argo.yaml @@ -0,0 +1,31 @@ +obs space: + name: salt_profile_argo + obsdatain: + engine: + type: H5File + obsfile: !ENV ${DATA}/obs/${OPREFIX}salt_profile_argo.${PDY}${cyc}.nc4 + obsdataout: + engine: + type: H5File + obsfile: !ENV ${DATA}/diags/salt_profile_argo.${PDY}${cyc}.nc4 + simulated variables: [salinity] + io pool: + max pool size: 1 +obs operator: + name: VertInterp + observation alias file: ./obsop_name_map.yaml + vertical coordinate: sea_water_depth + observation vertical coordinate: depth + interpolation method: linear +obs error: + covariance model: diagonal + +obs filter: + # Passivate obs where ocean fraction is > 90% + - filter: Domain Check + action: + name: passivate + where: + - variable: {name: GeoVaLs/sea_area_fraction} + maxvalue: 0.9 + diff --git a/parm/soca/obs/config/salt_profile_glider.yaml b/parm/soca/obs/config/salt_profile_glider.yaml new file mode 100644 index 000000000..a1504886b --- /dev/null +++ b/parm/soca/obs/config/salt_profile_glider.yaml @@ -0,0 +1,31 @@ +obs space: + name: salt_profile_glider + obsdatain: + engine: + type: H5File + obsfile: !ENV ${DATA}/obs/${OPREFIX}salt_profile_glider.${PDY}${cyc}.nc4 + obsdataout: + engine: + type: H5File + obsfile: !ENV ${DATA}/diags/salt_profile_glider.${PDY}${cyc}.nc4 + simulated variables: [salinity] + io pool: + max pool size: 1 +obs operator: + name: VertInterp + observation alias file: ./obsop_name_map.yaml + vertical coordinate: sea_water_depth + observation vertical coordinate: depth + interpolation method: linear +obs error: + covariance model: diagonal + +obs filter: + # Passivate obs where ocean fraction is > 90% + - filter: Domain Check + action: + name: passivate + where: + - variable: {name: GeoVaLs/sea_area_fraction} + maxvalue: 0.9 + diff --git a/parm/soca/obs/config/salt_profile_tesac.yaml b/parm/soca/obs/config/salt_profile_tesac.yaml new file mode 100644 index 000000000..53cc4496e --- /dev/null +++ b/parm/soca/obs/config/salt_profile_tesac.yaml @@ -0,0 +1,31 @@ +obs space: + name: salt_profile_tesac + obsdatain: + engine: + type: H5File + obsfile: !ENV ${DATA}/obs/${OPREFIX}salt_profile_tesac.${PDY}${cyc}.nc4 + obsdataout: + engine: + type: H5File + obsfile: !ENV ${DATA}/diags/salt_profile_tesac.${PDY}${cyc}.nc4 + simulated variables: [salinity] + io pool: + max pool size: 1 +obs operator: + name: VertInterp + observation alias file: ./obsop_name_map.yaml + vertical coordinate: sea_water_depth + observation vertical coordinate: depth + interpolation method: linear +obs error: + covariance model: diagonal + +obs filter: + # Passivate obs where ocean fraction is > 90% + - filter: Domain Check + action: + name: passivate + where: + - variable: {name: GeoVaLs/sea_area_fraction} + maxvalue: 0.9 + diff --git a/parm/soca/obs/config/temp_profile_argo.yaml b/parm/soca/obs/config/temp_profile_argo.yaml new file mode 100644 index 000000000..ff1b8ce9d --- /dev/null +++ b/parm/soca/obs/config/temp_profile_argo.yaml @@ -0,0 +1,28 @@ +obs space: + name: temp_profile_argo + obsdatain: + engine: + type: H5File + obsfile: !ENV ${DATA}/obs/${OPREFIX}temp_profile_argo.${PDY}${cyc}.nc4 + obsdataout: + engine: + type: H5File + obsfile: !ENV ${DATA}/diags/temp_profile_argo.${PDY}${cyc}.nc4 + simulated variables: [waterTemperature] + io pool: + max pool size: 1 +obs operator: + name: InsituTemperature + +obs error: + covariance model: diagonal + +obs filters: + # Passivate obs where ocean fraction is > 90% + - filter: Domain Check + action: + name: passivate + where: + - variable: {name: GeoVaLs/sea_area_fraction} + maxvalue: 0.9 + diff --git a/parm/soca/obs/config/temp_profile_bathy.yaml b/parm/soca/obs/config/temp_profile_bathy.yaml new file mode 100644 index 000000000..90c0569ea --- /dev/null +++ b/parm/soca/obs/config/temp_profile_bathy.yaml @@ -0,0 +1,28 @@ +obs space: + name: temp_profile_bathy + obsdatain: + engine: + type: H5File + obsfile: !ENV ${DATA}/obs/${OPREFIX}temp_profile_bathy.${PDY}${cyc}.nc4 + obsdataout: + engine: + type: H5File + obsfile: !ENV ${DATA}/diags/temp_profile_bathy.${PDY}${cyc}.nc4 + simulated variables: [waterTemperature] + io pool: + max pool size: 1 +obs operator: + name: InsituTemperature + +obs error: + covariance model: diagonal + +obs filters: + # Passivate obs where ocean fraction is > 90% + - filter: Domain Check + action: + name: passivate + where: + - variable: {name: GeoVaLs/sea_area_fraction} + maxvalue: 0.9 + diff --git a/parm/soca/obs/config/temp_profile_glider.yaml b/parm/soca/obs/config/temp_profile_glider.yaml new file mode 100644 index 000000000..721fccafb --- /dev/null +++ b/parm/soca/obs/config/temp_profile_glider.yaml @@ -0,0 +1,28 @@ +obs space: + name: temp_profile_glider + obsdatain: + engine: + type: H5File + obsfile: !ENV ${DATA}/obs/${OPREFIX}temp_profile_glider.${PDY}${cyc}.nc4 + obsdataout: + engine: + type: H5File + obsfile: !ENV ${DATA}/diags/temp_profile_glider.${PDY}${cyc}.nc4 + simulated variables: [waterTemperature] + io pool: + max pool size: 1 +obs operator: + name: InsituTemperature + +obs error: + covariance model: diagonal + +obs filters: + # Passivate obs where ocean fraction is > 90% + - filter: Domain Check + action: + name: passivate + where: + - variable: {name: GeoVaLs/sea_area_fraction} + maxvalue: 0.9 + diff --git a/parm/soca/obs/config/temp_profile_tesac.yaml b/parm/soca/obs/config/temp_profile_tesac.yaml new file mode 100644 index 000000000..8cbe84bd6 --- /dev/null +++ b/parm/soca/obs/config/temp_profile_tesac.yaml @@ -0,0 +1,28 @@ +obs space: + name: temp_profile_tesac + obsdatain: + engine: + type: H5File + obsfile: !ENV ${DATA}/obs/${OPREFIX}temp_profile_tesac.${PDY}${cyc}.nc4 + obsdataout: + engine: + type: H5File + obsfile: !ENV ${DATA}/diags/temp_profile_tesac.${PDY}${cyc}.nc4 + simulated variables: [waterTemperature] + io pool: + max pool size: 1 +obs operator: + name: InsituTemperature + +obs error: + covariance model: diagonal + +obs filters: + # Passivate obs where ocean fraction is > 90% + - filter: Domain Check + action: + name: passivate + where: + - variable: {name: GeoVaLs/sea_area_fraction} + maxvalue: 0.9 + diff --git a/parm/soca/obs/obs_list.yaml b/parm/soca/obs/obs_list.yaml index d710f3461..89b2659a1 100644 --- a/parm/soca/obs/obs_list.yaml +++ b/parm/soca/obs/obs_list.yaml @@ -13,3 +13,10 @@ observers: #- !INC ${OBS_YAML_DIR}/icec_ssmis_f17_south.yaml #- !INC ${OBS_YAML_DIR}/icec_ssmis_f18_north.yaml - !INC ${OBS_YAML_DIR}/icec_ssmis_f18_south.yaml +#- !INC ${OBS_YAML_DIR}/temp_profile_argo.yaml +#- !INC ${OBS_YAML_DIR}/temp_profile_bathy.yaml +#- !INC ${OBS_YAML_DIR}/temp_profile_glider.yaml +#- !INC ${OBS_YAML_DIR}/temp_profile_tesac.yaml +#- !INC ${OBS_YAML_DIR}/salt_profile_argo.yaml +#- !INC ${OBS_YAML_DIR}/salt_profile_glider.yaml +#- !INC ${OBS_YAML_DIR}/salt_profile_tesac.yaml \ No newline at end of file diff --git a/test/land/letkfoi_land.yaml b/test/land/letkfoi_land.yaml index c9069f0e8..98cd3d171 100644 --- a/test/land/letkfoi_land.yaml +++ b/test/land/letkfoi_land.yaml @@ -64,7 +64,7 @@ observations: where: - minvalue: '-999.0' variable: - name: MetaData/height + name: MetaData/stationElevation - filter: Domain Check where: - maxvalue: '1.5' diff --git a/test/soca/gw/run_jjobs.yaml.test b/test/soca/gw/run_jjobs.yaml.test index f773a8c79..f908c841f 100644 --- a/test/soca/gw/run_jjobs.yaml.test +++ b/test/soca/gw/run_jjobs.yaml.test @@ -33,6 +33,9 @@ gw environement: OOPS_DEBUG: 1 OMP_NUM_THREADS: 1 + run scripts: + GDASPREPOCNOBSPY: @HOMEgfs@/sorc/gdas.cd/scripts/exglobal_prep_ocean_obs.py + setup_expt config: base: DO_JEDIATMVAR: "NO" @@ -51,6 +54,8 @@ setup_expt config: SABER_BLOCKS_YAML: @HOMEgfs@/sorc/gdas.cd/parm/soca/berror/saber_blocks.yaml NICAS_RESOL: 1 NICAS_GRID_SIZE: 150 + prepoceanobs: + SOCA_OBS_LIST: @HOMEgfs@/sorc/gdas.cd/parm/soca/obs/obs_list_small.yaml job options: account: da-cpu diff --git a/test/soca/testinput/socahybridweights.yaml b/test/soca/testinput/socahybridweights.yaml index cdad1ab27..855dd6d60 100644 --- a/test/soca/testinput/socahybridweights.yaml +++ b/test/soca/testinput/socahybridweights.yaml @@ -16,8 +16,19 @@ background: read_from_file: 1 weights: - ice: 0.1 - ocean: 0.5 + # Need to provide weights^2 when reading from file + ice: 0.0025 # 5% of original variance + ocean: 0.0625 # 25% " " + # Apply localized weights to the ocean ens. B + ocean local weights: + - lon: -172.0 + lat: 11.0 + amplitude: -1.0 + length scale: 700.0 + - lon: -160.0 + lat: 12.0 + amplitude: -1.0 + length scale: 700.0 output: datadir: ./ diff --git a/ush/eva/jedi_gsi_compare_conv.yaml b/ush/eva/jedi_gsi_compare_conv.yaml index 02293a7b0..c82a0dbd9 100644 --- a/ush/eva/jedi_gsi_compare_conv.yaml +++ b/ush/eva/jedi_gsi_compare_conv.yaml @@ -107,6 +107,24 @@ transforms: for: variable: *variables + # Obs Error that passed QC for JEDI + - transform: accept where + new name: experiment::EffectiveErrorQc::${variable} + starting field: experiment::EffectiveError::${variable} + where: + - experiment::EffectiveQC::${variable} == 0 + for: + variable: *variables + + # ObsError that passed QC for GSI + - transform: accept where + new name: experiment::GsiFinalObsErrorQc::${variable} + starting field: experiment::GsiFinalObsError::${variable} + where: + - experiment::GsiEffectiveQC::${variable} == 0 + for: + variable: *variables + graphics: # ---------- Scatter Plots ---------- @@ -276,6 +294,7 @@ graphics: markersize: 5 color: 'red' label: 'hofxdiff vs pressure' + do_linear_regression: False statistics: fields: - field_name: experiment::HofXDiff::${variable} @@ -354,6 +373,42 @@ graphics: markersize: 5 color: 'red' label: 'errordiff vs pressure' + do_linear_regression: False + # Error comparison as a function of pressure + - batch figure: + variables: *variables + @CHANNELSKEY@ + figure: + layout: [1,1] + title: 'JEDI-GSI Comparison (Passed QC) | @NAME@ @CYCLE@ | ${variable_title}' + output name: observation_scatter_plots/errorcomp_vs_pressure_@CYCLE@_@NAME@_${variable}@CHANNELVAR@.png + plots: + - add_xlabel: 'Observation Error' + add_ylabel: 'Observation Pressure' + add_grid: + add_legend: + loc: 'lower left' + layers: + - type: Scatter + x: + variable: experiment::GsiFinalObsErrorQc::${variable} + y: + variable: experiment::MetaData::pressure + @CHANNELKEY@ + markersize: 3 + color: 'red' + label: 'GSI Error' + do_linear_regression: False + - type: Scatter + x: + variable: experiment::EffectiveErrorQc::${variable} + y: + variable: experiment::MetaData::pressure + @CHANNELKEY@ + markersize: 5 + color: 'green' + label: 'JEDI Error' + do_linear_regression: False # ---------- Histograms ---------- # Histogram of h(x) difference diff --git a/ush/ioda/bufr2ioda/bufr2ioda_adpsfc_prepbufr.py b/ush/ioda/bufr2ioda/bufr2ioda_adpsfc_prepbufr.py index 0f2b60311..a5700574a 100755 --- a/ush/ioda/bufr2ioda/bufr2ioda_adpsfc_prepbufr.py +++ b/ush/ioda/bufr2ioda/bufr2ioda_adpsfc_prepbufr.py @@ -52,7 +52,7 @@ def bufr_to_ioda(config, logger): reference_time = reference_time.strftime("%Y-%m-%dT%H:%M:%SZ") reference_time_full = f"{yyyymmdd}{hh}00" - logger.info(f"reference_time = {reference_time}") + logger.debug(f"reference_time = {reference_time}") # General informaton converter = 'BUFR to IODA Converter' @@ -62,34 +62,38 @@ def bufr_to_ioda(config, logger): DATA_PATH = os.path.join(dump_dir, f"{cycle_type}.{yyyymmdd}", str(hh), bufrfile) - logger.info(f"The DATA_PATH is: {DATA_PATH}") + 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 ...') + logger.debug('Making QuerySet ...') q = bufr.QuerySet(subsets) + # ObsType + q.add('observationType', '*/TYP') + # MetaData q.add('stationIdentification', '*/SID') q.add('latitude', '*/YOB') q.add('longitude', '*/XOB') q.add('obsTimeMinusCycleTime', '*/DHR') - q.add('stationElevation', '*/ELV') - q.add('observationType', '*/TYP') + q.add('heightOfStation', '*/Z___INFO/Z__EVENT{1}/ZOB') q.add('pressure', '*/P___INFO/P__EVENT{1}/POB') # # Quality Infomation (Quality Indicator) q.add('qualityMarkerStationPressure', '*/P___INFO/P__EVENT{1}/PQM') + q.add('qualityMarkerStationElevation', '*/Z___INFO/Z__EVENT{1}/ZQM') # ObsValue q.add('stationPressure', '*/P___INFO/P__EVENT{1}/POB') + q.add('stationElevation', '*/ELV') end_time = time.time() running_time = end_time - start_time - logger.info(f"Running time for making QuerySet: {running_time} seconds") + logger.debug(f"Running time for making QuerySet: {running_time} seconds") # ============================================================== # Open the BUFR file and execute the QuerySet to get ResultSet @@ -97,65 +101,80 @@ def bufr_to_ioda(config, logger): # ============================================================== start_time = time.time() - logger.info(f"Executing QuerySet to get ResultSet ...") + logger.debug(f"Executing QuerySet to get ResultSet ...") with bufr.File(DATA_PATH) as f: r = f.execute(q) - logger.info(" ... Executing QuerySet: get MetaData: basic ...") + logger.debug(" ... Executing QuerySet: get ObsType ...") + # ObsType + typ = r.get('observationType') + + logger.debug(" ... Executing QuerySet: get MetaData ...") # MetaData sid = r.get('stationIdentification') lat = r.get('latitude') lon = r.get('longitude') lon[lon > 180] -= 360 - elv = r.get('stationElevation') - typ = r.get('observationType') + zob = r.get('heightOfStation', type='float') pressure = r.get('pressure') pressure *= 100 - logger.info(f" ... Executing QuerySet: get QualityMarker information ...") + logger.debug(f" ... Executing QuerySet: get QualityMarker information ...") # Quality Information pobqm = r.get('qualityMarkerStationPressure') + zobqm = r.get('qualityMarkerStationElevation') - logger.info(f" ... Executing QuerySet: get obsvalue: stationPressure ...") + logger.debug(f" ... Executing QuerySet: get ObsValue ...") # ObsValue + elv = r.get('stationElevation', type='float') pob = r.get('stationPressure') pob *= 100 - logger.info(f" ... Executing QuerySet: get datatime: observation time ...") + logger.debug(f" ... Executing QuerySet: get dateTime ...") # DateTime: seconds since Epoch time # IODA has no support for numpy datetime arrays dtype=datetime64[s] dhr = r.get('obsTimeMinusCycleTime', type='int64') - logger.info(f" ... Executing QuerySet: Done!") + logger.debug(f" ... Executing QuerySet: Done!") - logger.info(f" ... Executing QuerySet: Check BUFR variable generic \ + logger.debug(f" ... Executing QuerySet: Check BUFR variable generic \ dimension and type ...") # Check BUFR variable generic dimension and type - logger.info(f" sid shape = {sid.shape}") - logger.info(f" dhr shape = {dhr.shape}") - logger.info(f" lat shape = {lat.shape}") - logger.info(f" lon shape = {lon.shape}") - logger.info(f" elv shape = {elv.shape}") - logger.info(f" typ shape = {typ.shape}") - logger.info(f" pressure shape = {pressure.shape}") - - logger.info(f" pobqm shape = {pobqm.shape}") - logger.info(f" pob shape = {pob.shape}") - - logger.info(f" sid type = {sid.dtype}") - logger.info(f" dhr type = {dhr.dtype}") - logger.info(f" lat type = {lat.dtype}") - logger.info(f" lon type = {lon.dtype}") - logger.info(f" elv type = {elv.dtype}") - logger.info(f" typ type = {typ.dtype}") - logger.info(f" pressure type = {pressure.dtype}") - - logger.info(f" pobqm type = {pobqm.dtype}") - logger.info(f" pob type = {pob.dtype}") + logger.debug(f" typ shape = {typ.shape}") + logger.debug(f" sid shape = {sid.shape}") + logger.debug(f" dhr shape = {dhr.shape}") + logger.debug(f" lat shape = {lat.shape}") + logger.debug(f" lon shape = {lon.shape}") + logger.debug(f" zob shape = {zob.shape}") + logger.debug(f" pressure shape = {pressure.shape}") + + logger.debug(f" pobqm shape = {pobqm.shape}") + logger.debug(f" zobqm shape = {zobqm.shape}") + + logger.debug(f" elv shape = {elv.shape}") + logger.debug(f" pob shape = {pob.shape}") + + logger.debug(f" dhr type = {dhr.shape}") + + logger.debug(f" sid type = {sid.dtype}") + logger.debug(f" dhr type = {dhr.dtype}") + logger.debug(f" lat type = {lat.dtype}") + logger.debug(f" lon type = {lon.dtype}") + logger.debug(f" zob type = {zob.dtype}") + logger.debug(f" typ type = {typ.dtype}") + logger.debug(f" pressure type = {pressure.dtype}") + + logger.debug(f" pobqm type = {pobqm.dtype}") + logger.debug(f" zobqm type = {zobqm.dtype}") + + logger.debug(f" elv type = {elv.dtype}") + logger.debug(f" pob type = {pob.dtype}") + + logger.debug(f" dhr type = {dhr.dtype}") end_time = time.time() running_time = end_time - start_time - logger.info(f"Running time for executing QuerySet to get ResultSet: \ + logger.debug(f"Running time for executing QuerySet to get ResultSet: \ {running_time} seconds") # ========================= @@ -163,19 +182,19 @@ def bufr_to_ioda(config, logger): # ========================= start_time = time.time() - logger.info(f"Creating derived variables - dateTime ...") + logger.debug(f"Creating derived variables - dateTime ...") cycleTimeSinceEpoch = np.int64(calendar.timegm(time.strptime( reference_time_full, '%Y%m%d%H%M'))) dateTime = Compute_dateTime(cycleTimeSinceEpoch, dhr) - logger.info(f" Check derived variables type ... ") - logger.info(f" dateTime shape = {dateTime.shape}") - logger.info(f" dateTime type = {dateTime.dtype}") + logger.debug(f" Check derived variables type ... ") + logger.debug(f" dateTime shape = {dateTime.shape}") + logger.debug(f" dateTime type = {dateTime.dtype}") end_time = time.time() running_time = end_time - start_time - logger.info(f"Running time for creating derived variables: \ + logger.debug(f"Running time for creating derived variables: \ {running_time} seconds") # ===================================== @@ -188,7 +207,7 @@ def bufr_to_ioda(config, logger): iodafile = f"{cycle_type}.t{hh}z.{data_type}.{data_format}.nc" OUTPUT_PATH = os.path.join(ioda_dir, iodafile) - logger.info(f" ... ... Create OUTPUT file: {OUTPUT_PATH}") + logger.debug(f" ... ... Create OUTPUT file: {OUTPUT_PATH}") path, fname = os.path.split(OUTPUT_PATH) if path and not os.path.exists(path): @@ -197,7 +216,7 @@ def bufr_to_ioda(config, logger): obsspace = ioda_ospace.ObsSpace(OUTPUT_PATH, mode='w', dim_dict=dims) # Create Global attributes - logger.info(f" ... ... Create global attributes") + logger.debug(f" ... ... Create global attributes") obsspace.write_attr('Converter', converter) obsspace.write_attr('source', source) @@ -210,7 +229,20 @@ def bufr_to_ioda(config, logger): obsspace.write_attr('platformLongDescription', platform_description) # Create IODA variables - logger.info(f" ... ... Create variables: name, type, units, & attributes") + logger.debug(f" ... ... Create variables: name, type, units, & attributes") + + # Observation Type - Station Elevation + obsspace.create_var('ObsType/stationElevation', dtype=typ.dtype, + fillval=typ.fill_value) \ + .write_attr('long_name', 'Station Elevation Observation Type') \ + .write_data(typ) + + # Observation Type - Station Pressure + obsspace.create_var('ObsType/stationPressure', dtype=typ.dtype, + fillval=typ.fill_value) \ + .write_attr('long_name', 'Station Pressure Observation Type') \ + .write_data(typ) + # Longitude obsspace.create_var('MetaData/longitude', dtype=lon.dtype, fillval=lon.fill_value) \ @@ -240,18 +272,12 @@ def bufr_to_ioda(config, logger): .write_attr('long_name', 'Station Identification') \ .write_data(sid) - # Station Elevation - obsspace.create_var('MetaData/stationElevation', dtype=elv.dtype, - fillval=elv.fill_value) \ + # Height Of Station + obsspace.create_var('MetaData/heightOfStation', dtype=zob.dtype, + fillval=zob.fill_value) \ .write_attr('units', 'm') \ - .write_attr('long_name', 'Station Elevation') \ - .write_data(elv) - - # Observation Type - obsspace.create_var('MetaData/observationType', dtype=typ.dtype, - fillval=typ.fill_value) \ - .write_attr('long_name', 'Observation Type') \ - .write_data(typ) + .write_attr('long_name', 'Height Of Station') \ + .write_data(zob) # Pressure obsspace.create_var('MetaData/pressure', dtype=pressure.dtype, @@ -260,14 +286,27 @@ def bufr_to_ioda(config, logger): .write_attr('long_name', 'Pressure') \ .write_data(pressure) - # Quality: Percent Confidence - Quality Information Without Forecast + # QualityMarker - Station Elevation + obsspace.create_var('QualityMarker/stationElevation', dtype=zobqm.dtype, + fillval=zobqm.fill_value) \ + .write_attr('long_name', 'Station Elevation Quality Marker') \ + .write_data(zobqm) + + # QualityMarker - Station Pressure obsspace.create_var('QualityMarker/stationPressure', dtype=pobqm.dtype, fillval=pobqm.fill_value) \ .write_attr('long_name', 'Station Pressure Quality Marker') \ .write_data(pobqm) + # Station Elevation + obsspace.create_var('ObsValue/stationElevation', dtype=elv.dtype, + fillval=elv.fill_value) \ + .write_attr('units', 'm') \ + .write_attr('long_name', 'Station Elevation') \ + .write_data(elv) + # Station Pressure - obsspace.create_var('ObsValue/pressure', dtype=pob.dtype, + obsspace.create_var('ObsValue/stationPressure', dtype=pob.dtype, fillval=pob.fill_value) \ .write_attr('units', 'Pa') \ .write_attr('long_name', 'Station Pressure') \ @@ -275,10 +314,10 @@ def bufr_to_ioda(config, logger): end_time = time.time() running_time = end_time - start_time - logger.info(f"Running time for splitting and output IODA: \ + logger.debug(f"Running time for splitting and output IODA: \ {running_time} seconds") - logger.info("All Done!") + logger.debug("All Done!") if __name__ == '__main__': @@ -305,4 +344,4 @@ def bufr_to_ioda(config, logger): end_time = time.time() running_time = end_time - start_time - logger.info(f"Total running time: {running_time} seconds") + logger.debug(f"Total running time: {running_time} seconds") diff --git a/ush/ioda/bufr2ioda/bufr2ioda_sfcshp_prepbufr.py b/ush/ioda/bufr2ioda/bufr2ioda_sfcshp_prepbufr.py old mode 100644 new mode 100755 index 48439c774..47207a3e9 --- a/ush/ioda/bufr2ioda/bufr2ioda_sfcshp_prepbufr.py +++ b/ush/ioda/bufr2ioda/bufr2ioda_sfcshp_prepbufr.py @@ -52,37 +52,40 @@ def bufr_to_ioda(config, logger): reference_time = reference_time.strftime("%Y-%m-%dT%H:%M:%SZ") reference_time_full = f"{yyyymmdd}{hh}00" - logger.info(f"Reference time = {reference_time}") + logger.debug(f"Reference time = {reference_time}") # General informaton converter = 'BUFR to IODA Converter' platform_description = 'SFCSHP data from prepBUFR format' bufrfile = f"{cycle_type}.t{hh}z.{data_format}" - DATA_PATH = os.path.join(dump_dir, f"{cycle_type}.{yyyymmdd}", str(hh), - bufrfile) + DATA_PATH = os.path.join(dump_dir, f"{cycle_type}.{yyyymmdd}", + str(hh), bufrfile) - logger.info("The DATA_PATH is: {DATA_PATH}") + logger.debug(f"The DATA_PATH is: {DATA_PATH}") # ============================================ # Make the QuerySet for all the data we want # ============================================ start_time = time.time() - logger.info(f"Making QuerySet ...") + logger.debug("Making QuerySet ...") q = bufr.QuerySet(subsets) + # ObsType + q.add('observationType', '*/TYP') + # MetaData q.add('stationIdentification', '*/SID') q.add('latitude', '*/YOB') q.add('longitude', '*/XOB') q.add('obsTimeMinusCycleTime', '*/DHR') - q.add('stationElevation', '*/ELV') - q.add('observationType', '*/TYP') + q.add('heightOfStation', '*/Z___INFO/Z__EVENT{1}/ZOB') q.add('pressure', '*/P___INFO/P__EVENT{1}/POB') q.add('temperatureEventCode', '*/T___INFO/T__EVENT{1}/TPC') # # Quality Infomation (Quality Indicator) + q.add('qualityMarkerStationElevation', '*/Z___INFO/Z__EVENT{1}/ZQM') q.add('qualityMarkerStationPressure', '*/P___INFO/P__EVENT{1}/PQM') q.add('qualityMarkerAirTemperature', '*/T___INFO/T__EVENT{1}/TQM') q.add('qualityMarkerSpecificHumidity', '*/Q___INFO/Q__EVENT{1}/QQM') @@ -90,9 +93,9 @@ def bufr_to_ioda(config, logger): q.add('qualityMarkerSeaSurfaceTemperature', '*/SST_INFO/SSTEVENT{1}/SSTQM') # ObsValue + q.add('stationElevation', '*/ELV') q.add('stationPressure', '*/P___INFO/P__EVENT{1}/POB') q.add('airTemperature', '*/T___INFO/T__EVENT{1}/TOB') - q.add('virtualTemperature', '*/T___INFO/TVO') q.add('specificHumidity', '*/Q___INFO/Q__EVENT{1}/QOB') q.add('windNorthward', '*/W___INFO/W__EVENT{1}/VOB') q.add('windEastward', '*/W___INFO/W__EVENT{1}/UOB') @@ -100,7 +103,7 @@ def bufr_to_ioda(config, logger): end_time = time.time() running_time = end_time - start_time - logger.info(f"Running time for making QuerySet : {running_time} seconds") + logger.debug(f"Running time for making QuerySet: {running_time} seconds") # ============================================================== # Open the BUFR file and execute the QuerySet to get ResultSet @@ -108,103 +111,124 @@ def bufr_to_ioda(config, logger): # ============================================================== start_time = time.time() - logger.info(f"Executing QuerySet to get ResultSet ...") + logger.debug(f"Executing QuerySet to get ResultSet ...") with bufr.File(DATA_PATH) as f: r = f.execute(q) - logger.info(f" ... Executing QuerySet: get metadata: basic ...") + logger.debug(f" ... Executing QuerySet: get metadata: basic ...") + # ObsType + logger.debug(" ... Executing QuerySet: get ObsType ...") + typ = r.get('observationType') + # MetaData + logger.debug(" ... Executing QuerySet: MetaData ...") sid = r.get('stationIdentification') lat = r.get('latitude') lon = r.get('longitude') lon[lon > 180] -= 360 - elv = r.get('stationElevation', type='float') - typ = r.get('observationType') + zob = r.get('heightOfStation', type='float') pressure = r.get('pressure') pressure *= 100 tpc = r.get('temperatureEventCode') - logger.info(f" ... Executing QuerySet: get QualityMarker information ...") # Quality Information - pqm = r.get('qualityMarkerStationPressure') - tqm = r.get('qualityMarkerAirTemperature') - qqm = r.get('qualityMarkerSpecificHumidity') - wqm = r.get('qualityMarkerWindNorthward') + logger.debug(f" ... Executing QuerySet: QualityMarker ...") + zobqm = r.get('qualityMarkerStationElevation') + pobqm = r.get('qualityMarkerStationPressure') + tobqm = r.get('qualityMarkerAirTemperature') + tsenqm = np.full(tobqm.shape[0], tobqm.fill_value) + tsenqm = np.where(((tpc >= 1) & (tpc < 8)), tobqm, tsenqm) + tvoqm = np.full(tobqm.shape[0], tobqm.fill_value) + tvoqm = np.where((tpc == 8), tobqm, tvoqm) + qobqm = r.get('qualityMarkerSpecificHumidity') + wobqm = r.get('qualityMarkerWindNorthward') sstqm = r.get('qualityMarkerSeaSurfaceTemperature') - logger.info(f" ... Executing QuerySet: get obsvalue: stationPressure ...") + logger.debug(f" ... Executing QuerySet: ObsValue ...") # ObsValue + elv = r.get('stationElevation', type='float') pob = r.get('stationPressure') pob *= 100 tob = r.get('airTemperature') tob += 273.15 - tvo = r.get('virtualTemperature') - tvo += 273.15 + tsen = np.full(tob.shape[0], tob.fill_value) + tsen = np.where(((tpc >= 1) & (tpc < 8)), tob, tsen) + tvo = np.full(tob.shape[0], tob.fill_value) + tvo = np.where((tpc == 8), tob, tvo) qob = r.get('specificHumidity', type='float') qob *= 0.000001 uob = r.get('windEastward') vob = r.get('windNorthward') sst1 = r.get('seaSurfaceTemperature') - logger.info(f" ... Executing QuerySet: get datatime: observation time ...") + logger.debug(f" ... Executing QuerySet: get datatime: observation time ...") # DateTime: seconds since Epoch time # IODA has no support for numpy datetime arrays dtype=datetime64[s] dhr = r.get('obsTimeMinusCycleTime', type='int64') - logger.info(f" ... Executing QuerySet: Done!") + logger.debug(f" ... Executing QuerySet: Done!") - logger.info(f" ... Executing QuerySet: Check BUFR variable generic \ + logger.debug(f" ... Executing QuerySet: Check BUFR variable generic \ dimension and type ...") # Check BUFR variable generic dimension and type - logger.info(f" sid shape = {sid.shape}") - logger.info(f" dhr shape = {dhr.shape}") - logger.info(f" lat shape = {lat.shape}") - logger.info(f" lon shape = {lon.shape}") - logger.info(f" elv shape = {elv.shape}") - logger.info(f" typ shape = {typ.shape}") - logger.info(f" pressure shape = {pressure.shape}") - logger.info(f" tpc shape = {tpc.shape}") - - logger.info(f" pqm shape = {pqm.shape}") - logger.info(f" tqm shape = {tqm.shape}") - logger.info(f" qqm shape = {qqm.shape}") - logger.info(f" wqm shape = {wqm.shape}") - logger.info(f" sstqm shape = {sstqm.shape}") - - logger.info(f" pob shape = {pob.shape}") - logger.info(f" tob shape = {pob.shape}") - logger.info(f" tvo shape = {tvo.shape}") - logger.info(f" qob shape = {qob.shape}") - logger.info(f" uob shape = {uob.shape}") - logger.info(f" vob shape = {vob.shape}") - logger.info(f" sst1 shape = {sst1.shape}") - - logger.info(f" sid type = {sid.dtype}") - logger.info(f" dhr type = {dhr.dtype}") - logger.info(f" lat type = {lat.dtype}") - logger.info(f" lon type = {lon.dtype}") - logger.info(f" elv type = {elv.dtype}") - logger.info(f" typ type = {typ.dtype}") - logger.info(f" pressure type = {pressure.dtype}") - logger.info(f" tpc type = {tpc.dtype}") - - logger.info(f" pqm type = {pqm.dtype}") - logger.info(f" tqm type = {tqm.dtype}") - logger.info(f" qqm type = {qqm.dtype}") - logger.info(f" wqm type = {wqm.dtype}") - logger.info(f" sstqm type = {sstqm.dtype}") - - logger.info(f" pob type = {pob.dtype}") - logger.info(f" tob type = {tob.dtype}") - logger.info(f" tvo type = {tvo.dtype}") - logger.info(f" qob type = {qob.dtype}") - logger.info(f" uob type = {uob.dtype}") - logger.info(f" vob type = {vob.dtype}") - logger.info(f" sst1 type = {sst1.dtype}") + logger.debug(f" typ shape = {typ.shape}") + logger.debug(f" sid shape = {sid.shape}") + logger.debug(f" dhr shape = {dhr.shape}") + logger.debug(f" lat shape = {lat.shape}") + logger.debug(f" lon shape = {lon.shape}") + logger.debug(f" zob shape = {zob.shape}") + logger.debug(f" pressure shape = {pressure.shape}") + logger.debug(f" tpc shape = {tpc.shape}") + + logger.debug(f" zobqm shape = {zobqm.shape}") + logger.debug(f" pobqm shape = {pobqm.shape}") + logger.debug(f" tobqm shape = {tobqm.shape}") + logger.debug(f" tsenqm shape = {tsenqm.shape}") + logger.debug(f" tvoqm shape = {tvoqm.shape}") + logger.debug(f" qobqm shape = {qobqm.shape}") + logger.debug(f" wobqm shape = {wobqm.shape}") + logger.debug(f" sstqm shape = {sstqm.shape}") + + logger.debug(f" elv shape = {elv.shape}") + logger.debug(f" pob shape = {pob.shape}") + logger.debug(f" tob shape = {tob.shape}") + logger.debug(f" tsen shape = {tsen.shape}") + logger.debug(f" tvo shape = {tvo.shape}") + logger.debug(f" qob shape = {qob.shape}") + logger.debug(f" uob shape = {uob.shape}") + logger.debug(f" vob shape = {vob.shape}") + logger.debug(f" sst1 shape = {sst1.shape}") + + logger.debug(f" typ type = {typ.dtype}") + logger.debug(f" sid type = {sid.dtype}") + logger.debug(f" dhr type = {dhr.dtype}") + logger.debug(f" lat type = {lat.dtype}") + logger.debug(f" lon type = {lon.dtype}") + logger.debug(f" zob type = {zob.dtype}") + logger.debug(f" pressure type = {pressure.dtype}") + logger.debug(f" tpc type = {tpc.dtype}") + + logger.debug(f" pobqm type = {pobqm.dtype}") + logger.debug(f" tobqm type = {tobqm.dtype}") + logger.debug(f" tsenqm type = {tsenqm.dtype}") + logger.debug(f" tvoqm type = {tvoqm.dtype}") + logger.debug(f" qobqm type = {qobqm.dtype}") + logger.debug(f" wobqm type = {wobqm.dtype}") + logger.debug(f" sstqm type = {sstqm.dtype}") + + logger.debug(f" elv type = {elv.dtype}") + logger.debug(f" pob type = {pob.dtype}") + logger.debug(f" tob type = {tob.dtype}") + logger.debug(f" tsen type = {tsen.dtype}") + logger.debug(f" tvo type = {tvo.dtype}") + logger.debug(f" qob type = {qob.dtype}") + logger.debug(f" uob type = {uob.dtype}") + logger.debug(f" vob type = {vob.dtype}") + logger.debug(f" sst1 type = {sst1.dtype}") end_time = time.time() running_time = end_time - start_time - logger.info(f"Running time for executing QuerySet to get ResultSet : \ + logger.debug(f"Running time for executing QuerySet to get ResultSet: \ {running_time} seconds") # ========================= @@ -212,19 +236,19 @@ def bufr_to_ioda(config, logger): # ========================= start_time = time.time() - logger.info(f"Creating derived variables - dateTime ...") + logger.debug(f"Creating derived variables - dateTime ...") cycleTimeSinceEpoch = np.int64(calendar.timegm(time.strptime( reference_time_full, '%Y%m%d%H%M'))) dateTime = Compute_dateTime(cycleTimeSinceEpoch, dhr) - logger.info(f" Check derived variables type ... ") - logger.info(f" dateTime shape = {dateTime.shape}") - logger.info(f" dateTime type = {dateTime.dtype}") + logger.debug(f" Check derived variables type ... ") + logger.debug(f" dateTime shape = {dateTime.shape}") + logger.debug(f" dateTime type = {dateTime.dtype}") end_time = time.time() running_time = end_time - start_time - logger.info(f"Running time for creating derived variables: {running_time} \ + logger.debug(f"Running time for creating derived variables: {running_time} \ seconds") # ===================================== @@ -237,7 +261,7 @@ def bufr_to_ioda(config, logger): iodafile = f"{cycle_type}.t{hh}z.{data_type}.{data_format}.nc" OUTPUT_PATH = os.path.join(ioda_dir, iodafile) - logger.info(f" ... ... Create OUTPUT file: {OUTPUT_PATH}") + logger.debug(f" ... ... Create OUTPUT file: {OUTPUT_PATH}") path, fname = os.path.split(OUTPUT_PATH) if path and not os.path.exists(path): @@ -246,7 +270,7 @@ def bufr_to_ioda(config, logger): obsspace = ioda_ospace.ObsSpace(OUTPUT_PATH, mode='w', dim_dict=dims) # Create Global attributes - logger.info(f" ... ... Create global attributes") + logger.debug(f" ... ... Create global attributes") obsspace.write_attr('Converter', converter) obsspace.write_attr('source', source) @@ -254,12 +278,61 @@ def bufr_to_ioda(config, logger): obsspace.write_attr('dataProviderOrigin', data_provider) obsspace.write_attr('description', data_description) obsspace.write_attr('datetimeReference', reference_time) - obsspace.write_attr('datetimeRange', [str(min(dateTime)), - str(max(dateTime))]) + obsspace.write_attr('datetimeRange', + [str(min(dateTime)), str(max(dateTime))]) obsspace.write_attr('platformLongDescription', platform_description) # Create IODA variables - logger.info(f" ... ... Create variables: name, type, units, & attributes") + logger.debug(f" ... ... Create variables: name, type, units, & attributes") + + # ObsType: stationElevation + obsspace.create_var('ObsType/stationElevation', dtype=typ.dtype, + fillval=typ.fill_value) \ + .write_attr('long_name', 'Station Elevation Observation Type') \ + .write_data(typ) + + # ObsType: stationPressure + obsspace.create_var('ObsType/stationPressure', dtype=typ.dtype, + fillval=typ.fill_value) \ + .write_attr('long_name', 'Station Pressure Observation Type') \ + .write_data(typ) + + # ObsType: air Temperature + obsspace.create_var('ObsType/airTemperature', dtype=typ.dtype, + fillval=typ.fill_value) \ + .write_attr('long_name', 'Air Temperature Observation Type') \ + .write_data(typ) + + # ObsType: virtual Temperature + obsspace.create_var('ObsType/virtualTemperature', dtype=typ.dtype, + fillval=typ.fill_value) \ + .write_attr('long_name', 'Virtual Temperature Observation Type') \ + .write_data(typ) + + # ObsType: Specific Humidity + obsspace.create_var('ObsType/specificHumidity', dtype=typ.dtype, + fillval=typ.fill_value) \ + .write_attr('long_name', 'Specific Humidity Observation Type') \ + .write_data(typ) + + # ObsType: wind Eastward + obsspace.create_var('ObsType/windEastward', dtype=typ.dtype, + fillval=typ.fill_value) \ + .write_attr('long_name', 'Wind Eastward Observation Type') \ + .write_data(typ) + + # ObsType: wind Northward + obsspace.create_var('ObsType/windNorthward', dtype=typ.dtype, + fillval=typ.fill_value) \ + .write_attr('long_name', 'Wind Northward Observation Type') \ + .write_data(typ) + + # ObsType: Sea Surface Temperature + obsspace.create_var('ObsType/seaSurfaceTemperature', dtype=typ.dtype, + fillval=typ.fill_value) \ + .write_attr('long_name', 'Sea Surface Temperature Observation Type') \ + .write_data(typ) + # Longitude obsspace.create_var('MetaData/longitude', dtype=lon.dtype, fillval=lon.fill_value) \ @@ -286,22 +359,15 @@ def bufr_to_ioda(config, logger): # Station Identification obsspace.create_var('MetaData/stationIdentification', dtype=sid.dtype, fillval=sid.fill_value) \ - .write_attr('units', '1') \ .write_attr('long_name', 'Station Identification') \ .write_data(sid) # Station Elevation - obsspace.create_var('MetaData/stationElevation', dtype=elv.dtype, - fillval=elv.fill_value) \ + obsspace.create_var('MetaData/heightOfStation', dtype=zob.dtype, + fillval=zob.fill_value) \ .write_attr('units', 'm') \ - .write_attr('long_name', 'Station Elevation') \ - .write_data(elv) - - # Observation Type - obsspace.create_var('MetaData/observationType', dtype=typ.dtype, - fillval=typ.fill_value) \ - .write_attr('long_name', 'Observation Type') \ - .write_data(typ) + .write_attr('long_name', 'Height Of Station') \ + .write_data(zob) # Pressure obsspace.create_var('MetaData/pressure', dtype=pressure.dtype, @@ -314,31 +380,43 @@ def bufr_to_ioda(config, logger): obsspace.create_var('MetaData/temperatureEventCode', dtype=tpc.dtype, fillval=tpc.fill_value) \ .write_attr('long_name', 'Temperature Event Code') \ - .write_data(typ) + .write_data(tpc) # Quality Marker: Station Pressure - obsspace.create_var('QualityMarker/stationPressure', dtype=pqm.dtype, - fillval=pqm.fill_value) \ + obsspace.create_var('QualityMarker/stationPressure', dtype=pobqm.dtype, + fillval=pobqm.fill_value) \ .write_attr('long_name', 'Station Pressure Quality Marker') \ - .write_data(pqm) + .write_data(pobqm) # Quality Marker: Air Temperature - obsspace.create_var('QualityMarker/airTemperature', dtype=tqm.dtype, - fillval=tqm.fill_value) \ + obsspace.create_var('QualityMarker/airTemperature', dtype=tobqm.dtype, + fillval=tobqm.fill_value) \ .write_attr('long_name', 'Air Temperature Quality Marker') \ - .write_data(tqm) + .write_data(tsenqm) + + # Quality Marker: Virtual Temperature + obsspace.create_var('QualityMarker/virtualTemperature', dtype=tobqm.dtype, + fillval=tobqm.fill_value) \ + .write_attr('long_name', 'Virtual Temperature Quality Marker') \ + .write_data(tvoqm) # Quality Marker: Specific Humidity - obsspace.create_var('QualityMarker/specificHumidity', dtype=qqm.dtype, - fillval=qqm.fill_value) \ + obsspace.create_var('QualityMarker/specificHumidity', dtype=qobqm.dtype, + fillval=qobqm.fill_value) \ .write_attr('long_name', 'Specific Humidity Quality Marker') \ - .write_data(qqm) + .write_data(qobqm) + + # Quality Marker: Eastward Wind + obsspace.create_var('QualityMarker/windEastward', dtype=wobqm.dtype, + fillval=wobqm.fill_value) \ + .write_attr('long_name', 'Eastward Wind Quality Marker') \ + .write_data(wobqm) # Quality Marker: Northward Wind - obsspace.create_var('QualityMarker/windNorthward', dtype=wqm.dtype, - fillval=wqm.fill_value) \ + obsspace.create_var('QualityMarker/windNorthward', dtype=wobqm.dtype, + fillval=wobqm.fill_value) \ .write_attr('long_name', 'Northward Wind Quality Marker') \ - .write_data(wqm) + .write_data(wobqm) # Quality Marker: Sea Surface Temperature obsspace.create_var('QualityMarker/seaSurfaceTemperature', @@ -346,8 +424,15 @@ def bufr_to_ioda(config, logger): .write_attr('long_name', 'Sea Surface Temperature Quality Marker') \ .write_data(sstqm) + # Station Elevation + obsspace.create_var('ObsValue/stationElevation', dtype=elv.dtype, + fillval=elv.fill_value) \ + .write_attr('units', 'm') \ + .write_attr('long_name', 'Station Elevation') \ + .write_data(elv) + # Station Pressure - obsspace.create_var('ObsValue/pressure', dtype=pob.dtype, + obsspace.create_var('ObsValue/stationPressure', dtype=pob.dtype, fillval=pob.fill_value) \ .write_attr('units', 'Pa') \ .write_attr('long_name', 'Station Pressure') \ @@ -358,11 +443,11 @@ def bufr_to_ioda(config, logger): fillval=tob.fill_value) \ .write_attr('units', 'K') \ .write_attr('long_name', 'Air Temperature') \ - .write_data(tob) + .write_data(tsen) # Virtual Temperature - obsspace.create_var('ObsValue/virtualTemperature', dtype=tvo.dtype, - fillval=tvo.fill_value) \ + obsspace.create_var('ObsValue/virtualTemperature', dtype=tob.dtype, + fillval=tob.fill_value) \ .write_attr('units', 'K') \ .write_attr('long_name', 'Virtual Temperature') \ .write_data(tvo) @@ -397,10 +482,10 @@ def bufr_to_ioda(config, logger): end_time = time.time() running_time = end_time - start_time - logger.info(f"Running time for splitting and output IODA: {running_time} \ - seconds") + logger.debug(f"Running time for splitting and output IODA: \ + {running_time} seconds") - logger.info("All Done!") + logger.debug("All Done!") if __name__ == '__main__': @@ -409,7 +494,8 @@ def bufr_to_ioda(config, logger): parser = argparse.ArgumentParser() parser.add_argument('-c', '--config', type=str, - help='Input JSON configuration', required=True) + help='Input JSON configuration', + required=True) parser.add_argument('-v', '--verbose', help='print debug logging information', action='store_true') @@ -426,4 +512,4 @@ def bufr_to_ioda(config, logger): end_time = time.time() running_time = end_time - start_time - logger.info(f"Total running time: {running_time} seconds") + logger.debug(f"Total running time: {running_time} seconds") diff --git a/ush/ioda/bufr2ioda/run_bufr2ioda.py b/ush/ioda/bufr2ioda/run_bufr2ioda.py index d7a7bb130..da01dd131 100755 --- a/ush/ioda/bufr2ioda/run_bufr2ioda.py +++ b/ush/ioda/bufr2ioda/run_bufr2ioda.py @@ -30,7 +30,8 @@ def bufr2ioda(current_cycle, RUN, DMPDIR, config_template_dir, COM_OBS): } # Specify observation types to be processed by a script - BUFR_py = ["satwind_amv_goes", "satwind_scat", "adpupa_prepbufr"] + + BUFR_py = ["satwind_amv_goes", "satwind_scat", "adpupa_prepbufr", "adpsfc_prepbufr", "sfcshp_prepbufr"] for obtype in BUFR_py: logger.info(f"Convert {obtype}...") diff --git a/ush/ufoeval/run_ufo_hofx_test.sh b/ush/ufoeval/run_ufo_hofx_test.sh index ccc043bd4..230966284 100755 --- a/ush/ufoeval/run_ufo_hofx_test.sh +++ b/ush/ufoeval/run_ufo_hofx_test.sh @@ -33,6 +33,7 @@ run_filtering=YES run_eva=YES eva_stats_only=NO keep_output=NO +window_shift=False while getopts "c:hsxq" opt; do case $opt in @@ -99,7 +100,6 @@ if [ $run_filtering == NO ]; then else yamlpath=$GDASApp/parm/atm/obs/testing/${obtype}.yaml fi - exename=test_ObsFilters.x #-------------- Do not modify below this line ---------------- @@ -178,11 +178,13 @@ export OPREFIX=gdas.t${cyc}z export APREFIX=gdas.t${cyc}z export GPREFIX=gdas.t${gcyc}z +if [ $obtype == "conv_ps" ]; then window_shift=True; fi cat > $workdir/temp.yaml << EOF window begin: '{{ WINDOW_BEGIN | to_isotime }}' window end: '{{ WINDOW_END | to_isotime }}' observations: - !INC $yamlpath +window shift: ${window_shift} EOF $GDASApp/ush/genYAML --input $workdir/temp.yaml --output $workdir/${obtype}_${cycle}.yaml diff --git a/ush/ufoeval/test_yamls.sh b/ush/ufoeval/test_yamls.sh index 7e98d80e2..d46056bd8 100755 --- a/ush/ufoeval/test_yamls.sh +++ b/ush/ufoeval/test_yamls.sh @@ -13,10 +13,23 @@ if [ $? -ne 0 ]; then exit 1 fi -for file in `ls ../../parm/atm/obs/testing/*yaml`; do +# Process tests wiht QC +for file in `find ../../parm/atm/obs/testing/*.yaml -type f -not -name "*noqc*"`; do basefile=${file##*/} - inst="${basefile%.*}" - ./run_ufo_hofx_test.sh -x $inst > $WORKDIR/$inst.log 2> $WORKDIR/$inst.err + obtype="${basefile%.*}" + ./run_ufo_hofx_test.sh -x $obtype > $WORKDIR/$obtype.log 2> $WORKDIR/$obtype.err + if [ $? == 0 ]; then + echo $basefile Passes \(yay!\) + else + echo $basefile Fails \(boo!\) + fi +done + +# Process tests without QC (HofX + Observation error assignment) +for file in `ls ../../parm/atm/obs/testing/*_noqc.yaml`; do + basefile=${file##*/} + obtype="${basefile%_noqc.*}" + ./run_ufo_hofx_test.sh -x -q $obtype > $WORKDIR/${obtype}_noqc.log 2> $WORKDIR/${obtype}_noqc.err if [ $? == 0 ]; then echo $basefile Passes \(yay!\) else diff --git a/utils/soca/gdas_ens_handler.h b/utils/soca/gdas_ens_handler.h index 44e994a0b..5cf862651 100644 --- a/utils/soca/gdas_ens_handler.h +++ b/utils/soca/gdas_ens_handler.h @@ -104,6 +104,10 @@ namespace gdasapp { gdasapp_ens_utils::ensMoments(ensMembers, ensMean, ensStd, ensVariance); oops::Log::info() << "mean: " << ensMean << std::endl; oops::Log::info() << "std: " << ensStd << std::endl; + if ( fullConfig.has("ensemble mean output") ) { + const eckit::LocalConfiguration ensMeanOutputConfig(fullConfig, "ensemble mean output"); + ensMean.write(ensMeanOutputConfig); + } // Remove mean from ensemble members for (size_t i = 0; i < postProcIncr.ensSize_; ++i) { @@ -115,6 +119,7 @@ namespace gdasapp { eckit::LocalConfiguration stericVarChangeConfig; fullConfig.get("steric height", stericVarChangeConfig); oops::Log::info() << "steric config 0000: " << stericVarChangeConfig << std::endl; + // Initialize trajectories const eckit::LocalConfiguration trajConfig(fullConfig, "trajectory"); soca::State cycleTraj(geom, trajConfig); // trajectory of the cycle @@ -148,7 +153,7 @@ namespace gdasapp { soca::Increment incr = postProcIncr.appendLayer(ensMembers[i]); // Save total ssh - oops::Log::info() << "ssh ensemble memnber " << i << std::endl; + oops::Log::info() << "ssh ensemble member " << i << std::endl; soca::Increment ssh_tmp(geom, socaSshVar, postProcIncr.dt_); ssh_tmp = ensMembers[i]; sshTotal.push_back(ssh_tmp); @@ -186,14 +191,14 @@ namespace gdasapp { // Add the unbalanced ssh to the recentered perturbation // this assumes ssh_u is independent of the trajectory - oops::Log::info() << "&&&&& before adding ssh_u " << incr << std::endl; + oops::Log::debug() << "&&&&& before adding ssh_u " << incr << std::endl; atlas::FieldSet incrFs; incr.toFieldSet(incrFs); atlas::FieldSet sshNonStericFs; sshNonSteric[i].toFieldSet(sshNonStericFs); util::addFieldSets(incrFs, sshNonStericFs); incr.fromFieldSet(incrFs); - oops::Log::info() << "&&&&& after adding ssh_u " << incr << std::endl; + oops::Log::debug() << "&&&&& after adding ssh_u " << incr << std::endl; // Save final perturbation, used in the offline EnVAR result = postProcIncr.save(incr, i+1); diff --git a/utils/soca/gdas_socahybridweights.h b/utils/soca/gdas_socahybridweights.h index adf20188a..219450a8e 100644 --- a/utils/soca/gdas_socahybridweights.h +++ b/utils/soca/gdas_socahybridweights.h @@ -3,12 +3,17 @@ #include #include #include +#include #include "eckit/config/LocalConfiguration.h" #include "atlas/field.h" +#include "atlas/util/Earth.h" +#include "atlas/util/Geometry.h" +#include "atlas/util/Point.h" #include "oops/base/PostProcessor.h" +#include "oops/generic/gc99.h" #include "oops/mpi/mpi.h" #include "oops/runs/Application.h" #include "oops/util/DateTime.h" @@ -20,6 +25,38 @@ #include "soca/State/State.h" namespace gdasapp { + // Create a simple mask based on a Gaussian function + void gaussianMask(const soca::Geometry & geom, soca::Increment & gaussIncr, + eckit::LocalConfiguration conf) { + // Get the 2D grid + std::vector lats; + std::vector lons; + bool halo = true; + geom.latlon(lats, lons, halo); + + // Prepare fieldset from increment + atlas::FieldSet gaussIncrFs; + gaussIncr.toFieldSet(gaussIncrFs); + + // Get the GC99 parameters from config + double amp = conf.getDouble("amplitude"); + double scale = conf.getDouble("length scale"); + const atlas::PointLonLat p0(conf.getDouble("lon"), conf.getDouble("lat")); + + // Recompute weights + for (auto & field : gaussIncrFs) { + oops::Log::info() << "---------- Field name: " << field.name() << std::endl; + auto view = atlas::array::make_view(field); + for (int jnode = 0; jnode < field.shape(0); ++jnode) { + atlas::PointLonLat p1(lons[jnode], lats[jnode]); + double d = atlas::util::Earth::distance(p0, p1)/1000.0; + for (int jlevel = 0; jlevel < field.shape(1); ++jlevel) { + view(jnode, jlevel) += amp * oops::gc99(d/scale); + } + } + } + gaussIncr.fromFieldSet(gaussIncrFs); + } class SocaHybridWeights : public oops::Application { public: @@ -45,8 +82,7 @@ namespace gdasapp { socaVars += socaOcnVars; /// Read the background - // TODO(guillaume): Use the ice extent to set the weights ... no clue if this is - // possible at this level + // TODO(guillaume): Use the ice extent to set the weights soca::State socaBkg(geom, socaVars, dt); const eckit::LocalConfiguration socaBkgConfig(fullConfig, "background"); socaBkg.read(socaBkgConfig); @@ -70,6 +106,17 @@ namespace gdasapp { /// Create fields of weights for the ocean soca::Increment socaOcnHW(geom, socaOcnVars, dt); socaOcnHW.ones(); + + /// Apply localized gaussians to the weights + eckit::LocalConfiguration localWeightsConfigs(fullConfig, "weights.ocean local weights"); + std::vector localWeightsList = + localWeightsConfigs.getSubConfigurations(); + for (auto & conf : localWeightsList) { + gaussianMask(geom, socaOcnHW, conf); + oops::Log::info() << "Local weights for socaOcnHW: " << std::endl << conf << std::endl; + oops::Log::info() << socaOcnHW << std::endl; + } + socaOcnHW *= wOcean; oops::Log::info() << "socaOcnHW: " << std::endl << socaOcnHW << std::endl; socaOcnHW.write(socaHWOutConfig);