diff --git a/parm/atm/obs/config/satwind_viirs_n20.yaml.j2 b/parm/atm/obs/config/satwind_viirs_n20.yaml.j2 new file mode 100644 index 000000000..2a0c2258a --- /dev/null +++ b/parm/atm/obs/config/satwind_viirs_n20.yaml.j2 @@ -0,0 +1,289 @@ +- obs space: + name: satwind_viirs_n20 + obsdatain: + engine: + type: H5File + obsfile: '{{ DATA }}/obs/{{ OPREFIX }}satwnd.viirs_n20.tm00.nc' + obsdataout: + engine: + type: H5File + obsfile: '{{ DATA }}/diags/diag_satwind_viirs_n20_{{ current_cycle | to_YMDH }}.nc' + io pool: + max pool size: 1 + simulated variables: [windEastward, windNorthward] + + obs operator: + name: VertInterp + hofx scaling field: SurfaceWindScalingPressure + hofx scaling field group: DerivedVariables + + linear obs operator: + name: VertInterp + + # NOTE: Tests using the Gaussian Thinning filter (below) to duplicate GSI's thinning of AHI/Himawari-8 satwinds + # results in more JEDI satwinds in the diag file than in GSI, but far fewer JEDI satwinds assimilated than + # GSI. JEDI under-counts assimilated winds by roughly 25-40%, relative to GSI, and this under-count is not + # even including the temporal thinning which is applied in GSI but not JEDI (by this filter below). See + # GDASApp Issue #741 for details: https://github.com/NOAA-EMC/GDASApp/issues/741 + #obs pre filters: + #- filter: Gaussian Thinning + # horizontal_mesh: 200 + # vertical_mesh: 10000 + # use_reduced_horizontal_grid: true + # round_horizontal_bin_count_to_nearest: true + # partition_longitude_bins_using_mesh: true + + obs prior filters: + # Apply variable changes needed for wind scaling + # For wind observations with pressure provided + - filter: Variable Transforms + Transform: SurfaceWindScalingPressure + SkipWhenNoObs: False + + # Calculate error inflation factor for duplicate observations + #- filter: Variable Assignment + # assignments: + # - name: ObsErrorFactorDuplicateCheck/windEastward + # type: float + # function: + # name: ObsFunction/ObsErrorFactorDuplicateCheck + # options: + # use_air_pressure: true + # variable: windEastward + + #- filter: Variable Assignment + # assignments: + # - name: ObsErrorFactorDuplicateCheck/windNorthward + # type: float + # function: + # name: ObsFunction/ObsErrorFactorDuplicateCheck + # options: + # use_air_pressure: true + # variable: windNorthward + + obs post filters: + # Assign the initial observation error, based on height/pressure + # Hard-wiring to prepobs_errtable.global by Type + # ObsError is currently not updating in diag file, but passes directly to EffectiveError when no inflation is specified in YAML + + # Type 260 (VIIRS LWIR) + - filter: Perform Action + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 260 + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + xvar: + name: MetaData/pressure + xvals: [110000.,105000.,100000.,95000.,90000.,85000.,80000.,75000.,70000.,65000.,60000.,55000.,50000.,45000.,40000.,35000.,30000.,25000.,20000.,15000.,10000.,7500.,5000.,4000.,3000.,2000.,1000.,500.,400.,300.,200.,100.,0.] + errors: [3.8,3.8,3.8,3.8,3.8,3.8,3.8,3.8,3.9,3.9,4.,4.,4.1,5.,6.,6.3,6.6,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.] + + # sanity-check criteria + # Observation Range Sanity Check + # NOT EXPLICITLY CLEARED: No obs in this range in file, so 0 Bounds Check rejects (which is correct) but essentially untested + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + minvalue: -130. + maxvalue: 130. + action: + name: reject + + # Velocity Sanity Check + # NOT EXPLICITLY CLEARED: No obs in this range in file, so 0 Bounds Check rejects (which is correct) but essentially untested + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/Velocity + maxvalue: 130. + action: + name: reject + + # Reject any observation with a /=0 surface type (non-water surface) within + # 200 hPa of the surface pressure (as part of the LNVD check). + # CLEARED: maxvalue is rejecting >, not >= as per a Perform Action, so threshold is unchanged + - filter: Difference Check + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: + name: GeoVaLs/water_area_fraction + maxvalue: 0.99 + reference: GeoVaLs/surface_pressure + value: MetaData/pressure + maxvalue: -20000. # within 200 hPa above surface pressure, negative p-diff + action: + name: reject + + # LNVD check + # CLEARED: maxvalue is rejecting >, not >= as per a Perform Action, so threshold is unchanged + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/SatWindsLNVDCheck + maxvalue: 3. + action: + name: reject + + # GSI setupw routine QC + # Reject any ob Type [240–260] when pressure greater than 950 mb. + # CLEARED: minvalue/maxvalue are >=/<=, not >/<, so editing range by 1 Pa + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 240-260 + test variables: + - name: MetaData/pressure + maxvalue: 95001. + action: + name: reject + + # Multiple satellite platforms, reject when pressure is more than 50 mb above tropopause. + # CLEARED: minvalue is rejecting <, not <= as per a Perform Action, so threshold is unchanged + # Notes (eliu): This tropopause check reject too many obs; probably due to tropopause pressure estimation + # Turn this check off for now. + # Need to check if troposphere pressure was implemented correctly in fv3-jed + - filter: Difference Check + filter variables: + - name: windEastward + - name: windNorthward + reference: GeoVaLs/tropopause_pressure + value: MetaData/pressure + minvalue: -5000. # 50 hPa above tropopause level, negative p-diff + action: + name: reject + + # All satwinds must adjust errors based on ObsErrorFactorPressureCheck + # prior to the SPDB check (i.e. the gross-error check). The gross-error + # check uses the adjusted errors for error-bound tightening and rejection, + # so this check has to come first. This check will inflate errors for obs + # that are too close to either the model top or bottom. + # Notes (eliu): GMAO added a required parameter: adjusted_error_name. + - filter: Perform Action + filter variables: + - name: windEastward + where: + - variable: + name: ObsType/windEastward + is_in: 240-260 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck + options: + surface_obs: false + variable: windEastward + inflation factor: 4.0 + + - filter: Perform Action + filter variables: + - name: windNorthward + where: + - variable: + name: ObsType/windNorthward + is_in: 240-260 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck + options: + variable: windNorthward + inflation factor: 4.0 + + # All satwinds subject to a gross-error check that contains significant + # modifiers for satwinds with a negative speed-bias. ALL wind gross-error + # checks are currently being done by the SatWindsSPDBCheck. + # CLEARED + - filter: Background Check + filter variables: + - name: windEastward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260] + cgross: [ 2.5, 2.5, 2.5, 1.5, 2.5, 1.3, 1.3, 2.5, 2.5, 2.5, 2.5, 1.3, 2.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5] + error_min: [1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4] + error_max: [6.1, 6.1, 15.0, 15.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.1, 20.1, 20.1, 20.1, 20.1, 20.1] + variable: windEastward + action: + name: reject + + - filter: Background Check + filter variables: + - name: windNorthward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260] + cgross: [ 2.5, 2.5, 2.5, 1.5, 2.5, 1.3, 1.3, 2.5, 2.5, 2.5, 2.5, 1.3, 2.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5] + error_min: [1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4] + error_max: [6.1, 6.1, 15.0, 15.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.1, 20.1, 20.1, 20.1, 20.1, 20.1] + variable: windNorthward + action: + name: reject + + # The last error inflation check is for duplicate observations. This one needs + # to come last, because we don't want to inflate errors for duplication if one + # of the duplicates should be rejected. + # Notes (eliu): ObsErrorFactorDuplicateCheck obsfunction requires PreUseFlag (usage parameter from read_satwnd.f90). + # : Turn off duplicate check for now. + #- filter: Perform Action + # filter variables: + # - name: windEastward + # action: + # name: inflate error + # inflation variable: + # name: ObsErrorFactorDuplicateCheck/windEastward + + #- filter: Perform Action + # filter variables: + # - name: windNorthward + # action: + # name: inflate error + # inflation variable: + # name: ObsErrorFactorDuplicateCheck/windNorthward + + # We are extending this to an additional filter that inflates final ob-errors across-the-board by + # 1/0.8 = 1.25. This is caused by the GSI value of nvqc being set to .true. in the global operational + # configuration, see: https://github.com/NOAA-EMC/global-workflow/blob/d5ae3328fa4041b177357b1133f6b92e81c859d7/scripts/exglobal_atmos_analysis.sh#L750 + # This setting activates Line 1229 of setupw.f90 to scale ratio_errors by 0.8, which is applied in + # the denominator of the final ob-error, so 1/0.8 = 1.25 factor of ob-error inflation. + # + # If this nvqc functionality were to be switched off (i.e. if variational qc were to be turned off), + # you would want to remove this last inflation filter. + #- filter: Perform Action + # filter variables: + # - name: windEastward + # where: + # - variable: ObsType/windEastward + # is_in: 240-260 + # action: + # name: inflate error + # inflation factor: 1.25 + + #- filter: Perform Action + # filter variables: + # - name: windNorthward + # where: + # - variable: ObsType/windNorthward + # is_in: 240-260 + # action: + # name: inflate error + # inflation factor: 1.25 + + # End of Filters diff --git a/parm/atm/obs/config/satwind_viirs_npp.yaml.j2 b/parm/atm/obs/config/satwind_viirs_npp.yaml.j2 new file mode 100644 index 000000000..ce8a377a2 --- /dev/null +++ b/parm/atm/obs/config/satwind_viirs_npp.yaml.j2 @@ -0,0 +1,289 @@ +- obs space: + name: satwind_viirs_npp + obsdatain: + engine: + type: H5File + obsfile: '{{ DATA }}/obs/{{ OPREFIX }}satwnd.viirs_npp.tm00.nc' + obsdataout: + engine: + type: H5File + obsfile: '{{ DATA }}/diags/diag_satwind_viirs_npp_{{ current_cycle | to_YMDH }}.nc' + io pool: + max pool size: 1 + simulated variables: [windEastward, windNorthward] + + obs operator: + name: VertInterp + hofx scaling field: SurfaceWindScalingPressure + hofx scaling field group: DerivedVariables + + linear obs operator: + name: VertInterp + + # NOTE: Tests using the Gaussian Thinning filter (below) to duplicate GSI's thinning of AHI/Himawari-8 satwinds + # results in more JEDI satwinds in the diag file than in GSI, but far fewer JEDI satwinds assimilated than + # GSI. JEDI under-counts assimilated winds by roughly 25-40%, relative to GSI, and this under-count is not + # even including the temporal thinning which is applied in GSI but not JEDI (by this filter below). See + # GDASApp Issue #741 for details: https://github.com/NOAA-EMC/GDASApp/issues/741 + #obs pre filters: + #- filter: Gaussian Thinning + # horizontal_mesh: 200 + # vertical_mesh: 10000 + # use_reduced_horizontal_grid: true + # round_horizontal_bin_count_to_nearest: true + # partition_longitude_bins_using_mesh: true + + obs prior filters: + # Apply variable changes needed for wind scaling + # For wind observations with pressure provided + - filter: Variable Transforms + Transform: SurfaceWindScalingPressure + SkipWhenNoObs: False + + # Calculate error inflation factor for duplicate observations + #- filter: Variable Assignment + # assignments: + # - name: ObsErrorFactorDuplicateCheck/windEastward + # type: float + # function: + # name: ObsFunction/ObsErrorFactorDuplicateCheck + # options: + # use_air_pressure: true + # variable: windEastward + + #- filter: Variable Assignment + # assignments: + # - name: ObsErrorFactorDuplicateCheck/windNorthward + # type: float + # function: + # name: ObsFunction/ObsErrorFactorDuplicateCheck + # options: + # use_air_pressure: true + # variable: windNorthward + + obs post filters: + # Assign the initial observation error, based on height/pressure + # Hard-wiring to prepobs_errtable.global by Type + # ObsError is currently not updating in diag file, but passes directly to EffectiveError when no inflation is specified in YAML + + # Type 260 (VIIRS LWIR) + - filter: Perform Action + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 260 + action: + name: assign error + error function: + name: ObsFunction/ObsErrorModelStepwiseLinear + options: + xvar: + name: MetaData/pressure + xvals: [110000.,105000.,100000.,95000.,90000.,85000.,80000.,75000.,70000.,65000.,60000.,55000.,50000.,45000.,40000.,35000.,30000.,25000.,20000.,15000.,10000.,7500.,5000.,4000.,3000.,2000.,1000.,500.,400.,300.,200.,100.,0.] + errors: [3.8,3.8,3.8,3.8,3.8,3.8,3.8,3.8,3.9,3.9,4.,4.,4.1,5.,6.,6.3,6.6,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.,7.] + + # sanity-check criteria + # Observation Range Sanity Check + # NOT EXPLICITLY CLEARED: No obs in this range in file, so 0 Bounds Check rejects (which is correct) but essentially untested + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + minvalue: -130. + maxvalue: 130. + action: + name: reject + + # Velocity Sanity Check + # NOT EXPLICITLY CLEARED: No obs in this range in file, so 0 Bounds Check rejects (which is correct) but essentially untested + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/Velocity + maxvalue: 130. + action: + name: reject + + # Reject any observation with a /=0 surface type (non-water surface) within + # 200 hPa of the surface pressure (as part of the LNVD check). + # CLEARED: maxvalue is rejecting >, not >= as per a Perform Action, so threshold is unchanged + - filter: Difference Check + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: + name: GeoVaLs/water_area_fraction + maxvalue: 0.99 + reference: GeoVaLs/surface_pressure + value: MetaData/pressure + maxvalue: -20000. # within 200 hPa above surface pressure, negative p-diff + action: + name: reject + + # LNVD check + # CLEARED: maxvalue is rejecting >, not >= as per a Perform Action, so threshold is unchanged + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + test variables: + - name: ObsFunction/SatWindsLNVDCheck + maxvalue: 3. + action: + name: reject + + # GSI setupw routine QC + # Reject any ob Type [240–260] when pressure greater than 950 mb. + # CLEARED: minvalue/maxvalue are >=/<=, not >/<, so editing range by 1 Pa + - filter: Bounds Check + filter variables: + - name: windEastward + - name: windNorthward + where: + - variable: ObsType/windEastward + is_in: 240-260 + test variables: + - name: MetaData/pressure + maxvalue: 95001. + action: + name: reject + + # Multiple satellite platforms, reject when pressure is more than 50 mb above tropopause. + # CLEARED: minvalue is rejecting <, not <= as per a Perform Action, so threshold is unchanged + # Notes (eliu): This tropopause check reject too many obs; probably due to tropopause pressure estimation + # Turn this check off for now. + # Need to check if troposphere pressure was implemented correctly in fv3-jed + - filter: Difference Check + filter variables: + - name: windEastward + - name: windNorthward + reference: GeoVaLs/tropopause_pressure + value: MetaData/pressure + minvalue: -5000. # 50 hPa above tropopause level, negative p-diff + action: + name: reject + + # All satwinds must adjust errors based on ObsErrorFactorPressureCheck + # prior to the SPDB check (i.e. the gross-error check). The gross-error + # check uses the adjusted errors for error-bound tightening and rejection, + # so this check has to come first. This check will inflate errors for obs + # that are too close to either the model top or bottom. + # Notes (eliu): GMAO added a required parameter: adjusted_error_name. + - filter: Perform Action + filter variables: + - name: windEastward + where: + - variable: + name: ObsType/windEastward + is_in: 240-260 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck + options: + surface_obs: false + variable: windEastward + inflation factor: 4.0 + + - filter: Perform Action + filter variables: + - name: windNorthward + where: + - variable: + name: ObsType/windNorthward + is_in: 240-260 + action: + name: inflate error + inflation variable: + name: ObsFunction/ObsErrorFactorPressureCheck + options: + variable: windNorthward + inflation factor: 4.0 + + # All satwinds subject to a gross-error check that contains significant + # modifiers for satwinds with a negative speed-bias. ALL wind gross-error + # checks are currently being done by the SatWindsSPDBCheck. + # CLEARED + - filter: Background Check + filter variables: + - name: windEastward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260] + cgross: [ 2.5, 2.5, 2.5, 1.5, 2.5, 1.3, 1.3, 2.5, 2.5, 2.5, 2.5, 1.3, 2.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5] + error_min: [1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4] + error_max: [6.1, 6.1, 15.0, 15.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.1, 20.1, 20.1, 20.1, 20.1, 20.1] + variable: windEastward + action: + name: reject + + - filter: Background Check + filter variables: + - name: windNorthward + function absolute threshold: + - name: ObsFunction/WindsSPDBCheck + options: + wndtype: [ 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260] + cgross: [ 2.5, 2.5, 2.5, 1.5, 2.5, 1.3, 1.3, 2.5, 2.5, 2.5, 2.5, 1.3, 2.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5] + error_min: [1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4] + error_max: [6.1, 6.1, 15.0, 15.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.1, 20.1, 20.1, 20.1, 20.1, 20.1] + variable: windNorthward + action: + name: reject + + # The last error inflation check is for duplicate observations. This one needs + # to come last, because we don't want to inflate errors for duplication if one + # of the duplicates should be rejected. + # Notes (eliu): ObsErrorFactorDuplicateCheck obsfunction requires PreUseFlag (usage parameter from read_satwnd.f90). + # : Turn off duplicate check for now. + #- filter: Perform Action + # filter variables: + # - name: windEastward + # action: + # name: inflate error + # inflation variable: + # name: ObsErrorFactorDuplicateCheck/windEastward + + #- filter: Perform Action + # filter variables: + # - name: windNorthward + # action: + # name: inflate error + # inflation variable: + # name: ObsErrorFactorDuplicateCheck/windNorthward + + # We are extending this to an additional filter that inflates final ob-errors across-the-board by + # 1/0.8 = 1.25. This is caused by the GSI value of nvqc being set to .true. in the global operational + # configuration, see: https://github.com/NOAA-EMC/global-workflow/blob/d5ae3328fa4041b177357b1133f6b92e81c859d7/scripts/exglobal_atmos_analysis.sh#L750 + # This setting activates Line 1229 of setupw.f90 to scale ratio_errors by 0.8, which is applied in + # the denominator of the final ob-error, so 1/0.8 = 1.25 factor of ob-error inflation. + # + # If this nvqc functionality were to be switched off (i.e. if variational qc were to be turned off), + # you would want to remove this last inflation filter. + #- filter: Perform Action + # filter variables: + # - name: windEastward + # where: + # - variable: ObsType/windEastward + # is_in: 240-260 + # action: + # name: inflate error + # inflation factor: 1.25 + + #- filter: Perform Action + # filter variables: + # - name: windNorthward + # where: + # - variable: ObsType/windNorthward + # is_in: 240-260 + # action: + # name: inflate error + # inflation factor: 1.25 + + # End of Filters diff --git a/parm/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.json b/parm/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.json new file mode 100644 index 000000000..eab19c6b2 --- /dev/null +++ b/parm/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.json @@ -0,0 +1,16 @@ +{ + "data_format" : "bufr_d", + "data_type" : "satwnd", + "cycle_type" : "{{ RUN }}", + "cycle_datetime" : "{{ current_cycle | to_YMDH }}", + "dump_directory" : "{{ DMPDIR }}", + "ioda_directory" : "{{ COM_OBS }}", + "subsets" : [ "NC005091" ], + "data_description" : "NC005091 S-NPP/NOAA-20 SATWIND VIIRS IR(LW)", + "data_provider" : "U.S. NOAA/NESDIS", + "sensor_info" : { "sensor_name": "VIIRS", "sensor_full_name": "Visible Infrared Imaging Radiometer Suite", "sensor_id": 999 }, + "satellite_info" : [ + { "satellite_name": "npp", "satellite_full_name": "SUOMI-NPP", "satellite_id": 224, "launch time": "YYYYMMDD" }, + { "satellite_name": "n20", "satellite_full_name": "NOAA-20", "satellite_id": 225, "launch time": "YYYYMMDD" }, + ] +} diff --git a/scripts/exgdas_global_marine_analysis_post.py b/scripts/exgdas_global_marine_analysis_post.py index 54cce2363..1376bdd52 100755 --- a/scripts/exgdas_global_marine_analysis_post.py +++ b/scripts/exgdas_global_marine_analysis_post.py @@ -50,6 +50,7 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]): cyc = str(os.getenv('cyc')).zfill(2) bcyc = str((int(cyc) - 3) % 24).zfill(2) gcyc = str((int(cyc) - 6) % 24).zfill(2) # previous cycle +cdatedt = datetime.strptime(cdate, '%Y%m%d%H') bdatedt = datetime.strptime(cdate, '%Y%m%d%H') - timedelta(hours=3) bdate = datetime.strftime(bdatedt, '%Y-%m-%dT%H:00:00Z') mdate = datetime.strftime(datetime.strptime(cdate, '%Y%m%d%H'), '%Y-%m-%dT%H:00:00Z') @@ -93,8 +94,12 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]): os.path.join(com_ocean_analysis, f'{RUN}.t{bcyc}z.ocngrid.nc')]) # Copy the CICE analysis restart -cdateice = pdy + '.' + cyc + '0000' -post_file_list.append([os.path.join(anl_dir, 'Data', f'{cdateice}.cice_model.res.nc'), +if os.getenv('DOIAU') == "YES": + cice_rst_date = bdatedt.strftime('%Y%m%d.%H%M%S') +else: + cice_rst_date = cdatedt.strftime('%Y%m%d.%H%M%S') + +post_file_list.append([os.path.join(anl_dir, 'Data', f'{cice_rst_date}.cice_model.res.nc'), os.path.join(com_ice_restart, f'{cdate}.cice_model_anl.res.nc')]) FileHandler({'copy': post_file_list}).sync() diff --git a/scripts/exgdas_global_marine_analysis_prep.py b/scripts/exgdas_global_marine_analysis_prep.py index d80bfc834..686b4c8e8 100755 --- a/scripts/exgdas_global_marine_analysis_prep.py +++ b/scripts/exgdas_global_marine_analysis_prep.py @@ -70,10 +70,6 @@ def find_clim_ens(input_date): staticsoca_dir = os.getenv('SOCA_INPUT_FIX_DIR') nmem_ens = 0 nmem_ens = int(os.getenv('NMEM_ENS')) -if os.getenv('DOHYBVAR') == "YES": - dohybvar = True -else: - dohybvar = False # create analysis directories diags = os.path.join(anl_dir, 'diags') # output dir for soca DA obs space @@ -89,12 +85,25 @@ def find_clim_ens(input_date): window_begin = datetime.strptime(os.getenv('PDY')+os.getenv('cyc'), '%Y%m%d%H') - half_assim_freq window_begin_iso = window_begin.strftime('%Y-%m-%dT%H:%M:%SZ') window_middle_iso = window_middle.strftime('%Y-%m-%dT%H:%M:%SZ') -fcst_begin = datetime.strptime(os.getenv('PDY')+os.getenv('cyc'), '%Y%m%d%H') RUN = os.getenv('RUN') cyc = os.getenv('cyc') gcyc = os.getenv('gcyc') PDY = os.getenv('PDY') +# hybrid-envar switch +if os.getenv('DOHYBVAR') == "YES": + dohybvar = True +else: + dohybvar = False + +# switch for the cycling type +if os.getenv('DOIAU') == "YES": + # forecast initialized at the begining of the DA window + fcst_begin = window_begin +else: + # forecast initialized at the middle of the DA window + fcst_begin = datetime.strptime(os.getenv('PDY')+os.getenv('cyc'), '%Y%m%d%H') + ################################################################################ # fetch observations diff --git a/ush/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.py b/ush/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.py new file mode 100755 index 000000000..18f87db04 --- /dev/null +++ b/ush/ioda/bufr2ioda/bufr2ioda_satwind_amv_viirs.py @@ -0,0 +1,470 @@ +#!/usr/bin/env python3 +import argparse +import numpy as np +import numpy.ma as ma +from pyiodaconv import bufr +import calendar +import json +import time +import math +import datetime +import os +from datetime import datetime +from pyioda import ioda_obs_space as ioda_ospace +from wxflow import Logger + +# ==================================================================== +# Satellite Winds (AMV) BUFR dump file for VIIRS/S-NPP,NOAA-20 +# ==================================================================== +# Subset | Spectral Band | Code (002023) | ObsType +# -------------------------------------------------------------------- +# NC005091 | IRLW (Freq < 5E+13) | Method 1 | 260 +# ==================================================================== + +# Define and initialize global variables +global float32_fill_value +global int32_fill_value +global int64_fill_value + +float32_fill_value = np.float32(0) +int32_fill_value = np.int32(0) +int64_fill_value = np.int64(0) + + +def Compute_WindComponents_from_WindDirection_and_WindSpeed(wdir, wspd): + + uob = (-wspd * np.sin(np.radians(wdir))).astype(np.float32) + vob = (-wspd * np.cos(np.radians(wdir))).astype(np.float32) + + return uob, vob + + +def Get_ObsType(swcm, chanfreq): + + obstype = swcm.copy() + + # Use numpy vectorized operations + obstype = np.where(swcm == 1, 260, obstype) # IRLW + + if not np.any(np.isin(obstype, [260])): + raise ValueError("Error: Unassigned ObsType found ... ") + + return obstype + + +def bufr_to_ioda(config, logger): + + subsets = config["subsets"] + logger.debug(f"Checking subsets = {subsets}") + + # Get parameters from configuration + subsets = config["subsets"] + data_format = config["data_format"] + data_type = config["data_type"] + data_description = config["data_description"] + data_provider = config["data_provider"] + cycle_type = config["cycle_type"] + dump_dir = config["dump_directory"] + ioda_dir = config["ioda_directory"] + cycle = config["cycle_datetime"] + yyyymmdd = cycle[0:8] + hh = cycle[8:10] + + satellite_info_array = config["satellite_info"] + sensor_name = config["sensor_info"]["sensor_name"] + sensor_full_name = config["sensor_info"]["sensor_full_name"] + sensor_id = config["sensor_info"]["sensor_id"] + + # Get derived parameters + yyyymmdd = cycle[0:8] + hh = cycle[8:10] + reference_time = datetime.strptime(cycle, "%Y%m%d%H") + reference_time = reference_time.strftime("%Y-%m-%dT%H:%M:%SZ") + + # General informaton + converter = 'BUFR to IODA Converter' + process_level = 'Level-2' + platform_description = 'TERRA/AQUA' + sensor_description = 'Moderate Resolution Imaging Spectroradiometer' + + logger.info(f'sensor_name = {sensor_name}') + logger.info(f'sensor_full_name = {sensor_full_name}') + logger.info(f'sensor_id = {sensor_id}') + logger.info(f'reference_time = {reference_time}') + + bufrfile = f"{cycle_type}.t{hh}z.{data_type}.tm00.{data_format}" + DATA_PATH = os.path.join(dump_dir, f"{cycle_type}.{yyyymmdd}", str(hh), 'atmos', bufrfile) + if not os.path.isfile(DATA_PATH): + logger.info(f"DATA_PATH {DATA_PATH} does not exist") + return + logger.debug(f"The DATA_PATH is: {DATA_PATH}") + + # ============================================ + # Make the QuerySet for all the data we want + # ============================================ + start_time = time.time() + + logger.info('Making QuerySet') + q = bufr.QuerySet(subsets) + + # MetaData + q.add('latitude', '*/CLATH') + q.add('longitude', '*/CLONH') + q.add('satelliteId', '*/SAID') + q.add('year', '*/YEAR') + q.add('month', '*/MNTH') + q.add('day', '*/DAYS') + q.add('hour', '*/HOUR') + q.add('minute', '*/MINU') + q.add('second', '*/SECO') + q.add('satelliteZenithAngle', '*/SAZA') + q.add('sensorCentralFrequency', '*/SCCF') + q.add('pressure', '*/PRLC[1]') + + # Processing Center + q.add('dataProviderOrigin', '*/OGCE[1]') + q.add('windGeneratingApplication', '*/AMVQIC/GNAPS') + +# # Quality Infomation (Quality Indicator w/o forecast) + q.add('qualityInformationWithoutForecast', '*/AMVQIC/PCCF') + + # Wind Retrieval Method Information + q.add('windComputationMethod', '*/SWCM') + + # ObsValue + q.add('windDirection', '*/WDIR') + q.add('windSpeed', '*/WSPD') + + end_time = time.time() + running_time = end_time - start_time + logger.debug(f'Processing time for making QuerySet : {running_time} seconds') + + # ============================================================== + # Open the BUFR file and execute the QuerySet to get ResultSet + # Use the ResultSet returned to get numpy arrays of the data + # ============================================================== + start_time = time.time() + + logger.info('Executing QuerySet to get ResultSet') + with bufr.File(DATA_PATH) as f: + try: + r = f.execute(q) + except Exception as err: + logger.info(f'Return with {err}') + return + + # MetaData + satid = r.get('satelliteId') + year = r.get('year') + month = r.get('month') + day = r.get('day') + hour = r.get('hour') + minute = r.get('minute') + second = r.get('second') + lat = r.get('latitude') + lon = r.get('longitude') + satzenang = r.get('satelliteZenithAngle') + pressure = r.get('pressure', type='float') + chanfreq = r.get('sensorCentralFrequency', type='float') + + # Processing Center + ogce = r.get('dataProviderOrigin') + ga = r.get('windGeneratingApplication') + + # Quality Information + qi = r.get('qualityInformationWithoutForecast', type='float') + # For TERRA/AQUA MODIS data, qi w/o forecast (qifn) is packaged in same + # vector of qi with ga = 5 (QI without forecast), and EE is + # packaged in same vector of qi with ga=7 (Estimated Error (EE) in m/s + # converted to a percent confidence)shape (4,nobs). Must conduct a + # search and extract the correct vector for gnap and qi + # 1. Find dimension-sizes of ga and qi (should be the same!) + gDim1, gDim2 = np.shape(ga) + qDim1, qDim2 = np.shape(qi) + logger.info(f'Generating Application and Quality Information SEARCH:') + logger.info(f'Dimension size of GNAP ({gDim1},{gDim2})') + logger.info(f'Dimension size of PCCF ({qDim1},{qDim2})') + # 2. Initialize gnap and qifn as None, and search for dimension of + # ga with values of 5. If the same column exists for qi, assign + # gnap to ga[:,i] and qifn to qi[:,i], else raise warning that no + # appropriate GNAP/PCCF combination was found + gnap = None + qifn = None + for i in range(gDim2): + if np.unique(ga[:, i].squeeze()) == 5: + if i <= qDim2: + logger.info(f'GNAP/PCCF found for column {i}') + gnap = ga[:, i].squeeze() + qifn = qi[:, i].squeeze() + else: + logger.info(f'ERROR: GNAP column {i} outside of PCCF dimension {qDim2}') + if (gnap is None) & (qifn is None): + logger.info(f'ERROR: GNAP == 5 NOT FOUND OR OUT OF PCCF DIMENSION-RANGE, WILL FAIL!') + # If EE is needed, key search on np.unique(ga[:,i].squeeze()) == 7 instead + + # Wind Retrieval Method Information + swcm = r.get('windComputationMethod') + # ObsValue + # Wind direction and Speed + wdir = r.get('windDirection', type='float') + wspd = r.get('windSpeed') + + # DateTime: seconds since Epoch time + # IODA has no support for numpy datetime arrays dtype=datetime64[s] + timestamp = r.get_datetime('year', 'month', 'day', 'hour', 'minute', 'second').astype(np.int64) + + # Check BUFR variable generic dimension and type + + # Global variables declaration + # Set global fill values + float32_fill_value = satzenang.fill_value + int32_fill_value = satid.fill_value + int64_fill_value = timestamp.fill_value.astype(np.int64) + + end_time = time.time() + running_time = end_time - start_time + logger.info(f'Processing time for executing QuerySet to get ResultSet : {running_time} seconds') + + # ========================= + # Create derived variables + # ========================= + start_time = time.time() + + logger.info('Creating derived variables') + logger.debug('Creating derived variables - wind components (uob and vob)') + + uob, vob = Compute_WindComponents_from_WindDirection_and_WindSpeed(wdir, wspd) + + logger.debug(f' uob min/max = {uob.min()} {uob.max()}') + logger.debug(f' vob min/max = {vob.min()} {vob.max()}') + + obstype = Get_ObsType(swcm, chanfreq) + + height = np.full_like(pressure, fill_value=pressure.fill_value, dtype=np.float32) + stnelev = np.full_like(pressure, fill_value=pressure.fill_value, dtype=np.float32) + + end_time = time.time() + running_time = end_time - start_time + logger.info(f'Processing time for creating derived variables : {running_time} seconds') + + # ===================================== + # Split output based on satellite id + # Create IODA ObsSpace + # Write IODA output + # ===================================== + logger.info('Create IODA ObsSpace and Write IODA output based on satellite ID') + + # Find unique satellite identifiers in data to process + unique_satids = np.unique(satid) + logger.info(f'Number of Unique satellite identifiers: {len(unique_satids)}') + logger.info(f'Unique satellite identifiers: {unique_satids}') + + logger.debug(f'Loop through unique satellite identifier {unique_satids}') + total_ob_processed = 0 + for sat in unique_satids.tolist(): + start_time = time.time() + + matched = False + for satellite_info in satellite_info_array: + if (satellite_info["satellite_id"] == sat): + matched = True + satellite_id = satellite_info["satellite_id"] + satellite_name = satellite_info["satellite_name"] + satinst = sensor_name.lower()+'_'+satellite_name.lower() + logger.debug(f'Split data for {satinst} satid = {sat}') + + if matched: + + # Define a boolean mask to subset data from the original data object + mask = satid == sat + # MetaData + lon2 = lon[mask] + lat2 = lat[mask] + timestamp2 = timestamp[mask] + satid2 = satid[mask] + satzenang2 = satzenang[mask] + chanfreq2 = chanfreq[mask] + obstype2 = obstype[mask] + pressure2 = pressure[mask] + height2 = height[mask] + stnelev2 = stnelev[mask] + + # Processing Center + ogce2 = ogce[mask] + + # QC Info + qifn2 = qifn[mask] + + # Method + swcm2 = swcm[mask] + + # ObsValue + wdir2 = wdir[mask] + wspd2 = wspd[mask] + uob2 = uob[mask] + vob2 = vob[mask] + + # Timestamp Range + timestamp2_min = datetime.fromtimestamp(timestamp2.min()) + timestamp2_max = datetime.fromtimestamp(timestamp2.max()) + + # Check unique observation time + unique_timestamp2 = np.unique(timestamp2) + logger.debug(f'Processing output for satid {sat}') + + # Create the dimensions + dims = { + 'Location': np.arange(0, wdir2.shape[0]) + } + + # Create IODA ObsSpace + iodafile = f"{cycle_type}.t{hh}z.{data_type}.{satinst}.tm00.nc" + OUTPUT_PATH = os.path.join(ioda_dir, iodafile) + logger.info(f'Create output file : {OUTPUT_PATH}') + obsspace = ioda_ospace.ObsSpace(OUTPUT_PATH, mode='w', dim_dict=dims) + + # Create Global attributes + logger.debug('Write global attributes') + obsspace.write_attr('Converter', converter) + obsspace.write_attr('sourceFiles', bufrfile) + obsspace.write_attr('dataProviderOrigin', data_provider) + obsspace.write_attr('description', data_description) + obsspace.write_attr('datetimeReference', reference_time) + obsspace.write_attr('datetimeRange', [str(timestamp2_min), str(timestamp2_max)]) + obsspace.write_attr('sensor', sensor_id) + obsspace.write_attr('platform', satellite_id) + obsspace.write_attr('platformCommonName', satellite_name) + obsspace.write_attr('sensorCommonName', sensor_name) + obsspace.write_attr('processingLevel', process_level) + obsspace.write_attr('platformLongDescription', platform_description) + obsspace.write_attr('sensorLongDescription', sensor_description) + + # Create IODA variables + logger.debug('Write variables: name, type, units, and attributes') + # Longitude + obsspace.create_var('MetaData/longitude', dtype=lon2.dtype, fillval=lon2.fill_value) \ + .write_attr('units', 'degrees_east') \ + .write_attr('valid_range', np.array([-180, 180], dtype=np.float32)) \ + .write_attr('long_name', 'Longitude') \ + .write_data(lon2) + + # Latitude + obsspace.create_var('MetaData/latitude', dtype=lat.dtype, fillval=lat2.fill_value) \ + .write_attr('units', 'degrees_north') \ + .write_attr('valid_range', np.array([-90, 90], dtype=np.float32)) \ + .write_attr('long_name', 'Latitude') \ + .write_data(lat2) + + # Datetime + obsspace.create_var('MetaData/dateTime', dtype=np.int64, fillval=int64_fill_value) \ + .write_attr('units', 'seconds since 1970-01-01T00:00:00Z') \ + .write_attr('long_name', 'Datetime') \ + .write_data(timestamp2) + + # Satellite Identifier + obsspace.create_var('MetaData/satelliteIdentifier', dtype=satid2.dtype, fillval=satid2.fill_value) \ + .write_attr('long_name', 'Satellite Identifier') \ + .write_data(satid2) + + # Sensor Zenith Angle + obsspace.create_var('MetaData/satelliteZenithAngle', dtype=satzenang2.dtype, fillval=satzenang2.fill_value) \ + .write_attr('units', 'degree') \ + .write_attr('valid_range', np.array([0, 90], dtype=np.float32)) \ + .write_attr('long_name', 'Satellite Zenith Angle') \ + .write_data(satzenang2) + + # Sensor Centrall Frequency + obsspace.create_var('MetaData/sensorCentralFrequency', dtype=chanfreq2.dtype, fillval=chanfreq2.fill_value) \ + .write_attr('units', 'Hz') \ + .write_attr('long_name', 'Satellite Channel Center Frequency') \ + .write_data(chanfreq2) + + # Data Provider + obsspace.create_var('MetaData/dataProviderOrigin', dtype=ogce2.dtype, fillval=ogce2.fill_value) \ + .write_attr('long_name', 'Identification of Originating/Generating Center') \ + .write_data(ogce2) + + # Quality: Percent Confidence - Quality Information Without Forecast + obsspace.create_var('MetaData/qualityInformationWithoutForecast', dtype=qifn2.dtype, fillval=qifn2.fill_value) \ + .write_attr('long_name', 'Quality Information Without Forecast') \ + .write_data(qifn2) + + # Wind Computation Method + obsspace.create_var('MetaData/windComputationMethod', dtype=swcm2.dtype, fillval=swcm2.fill_value) \ + .write_attr('long_name', 'Satellite-derived Wind Computation Method') \ + .write_data(swcm2) + + # Pressure + obsspace.create_var('MetaData/pressure', dtype=pressure2.dtype, fillval=pressure2.fill_value) \ + .write_attr('units', 'pa') \ + .write_attr('long_name', 'Pressure') \ + .write_data(pressure2) + + # Height (mimic prepbufr) + obsspace.create_var('MetaData/height', dtype=height2.dtype, fillval=height2.fill_value) \ + .write_attr('units', 'm') \ + .write_attr('long_name', 'Height of Observation') \ + .write_data(height2) + + # Station Elevation (mimic prepbufr) + obsspace.create_var('MetaData/stationElevation', dtype=stnelev2.dtype, fillval=stnelev2.fill_value) \ + .write_attr('units', 'm') \ + .write_attr('long_name', 'Station Elevation') \ + .write_data(stnelev2) + + # ObsType based on computation method/spectral band + obsspace.create_var('ObsType/windEastward', dtype=obstype2.dtype, fillval=swcm2.fill_value) \ + .write_attr('long_name', 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band') \ + .write_data(obstype2) + + # ObsType based on computation method/spectral band + obsspace.create_var('ObsType/windNorthward', dtype=obstype2.dtype, fillval=swcm2.fill_value) \ + .write_attr('long_name', 'Observation Type based on Satellite-derived Wind Computation Method and Spectral Band') \ + .write_data(obstype2) + + # U-Wind Component + obsspace.create_var('ObsValue/windEastward', dtype=uob2.dtype, fillval=uob2.fill_value) \ + .write_attr('units', 'm s-1') \ + .write_attr('long_name', 'Eastward Wind Component') \ + .write_data(uob2) + + # V-Wind Component + obsspace.create_var('ObsValue/windNorthward', dtype=vob2.dtype, fillval=vob2.fill_value) \ + .write_attr('units', 'm s-1') \ + .write_attr('long_name', 'Northward Wind Component') \ + .write_data(vob2) + + end_time = time.time() + running_time = end_time - start_time + total_ob_processed += len(satid2) + logger.debug(f'Number of observation processed : {len(satid2)}') + logger.debug(f'Processing time for splitting and output IODA for {satinst} : {running_time} seconds') + + else: + logger.info(f"Do not find this satellite id in the configuration: satid = {sat}") + + logger.info("All Done!") + logger.info(f'Total number of observation processed : {total_ob_processed}') + + +if __name__ == '__main__': + + start_time = time.time() + + parser = argparse.ArgumentParser() + parser.add_argument('-c', '--config', type=str, help='Input JSON configuration', required=True) + parser.add_argument('-v', '--verbose', help='print debug logging information', + action='store_true') + args = parser.parse_args() + + log_level = 'DEBUG' if args.verbose else 'INFO' + logger = Logger('BUFR2IODA_satwind_amv_ahi.py', level=log_level, colored_log=True) + + with open(args.config, "r") as json_file: + config = json.load(json_file) + + bufr_to_ioda(config, logger) + + end_time = time.time() + running_time = end_time - start_time + logger.info(f"Total running time: {running_time} seconds")