convert tropical storm central pressure from tcvitals files to IODA format

src/conventional/CMakeLists.txt

@@ -10,6 +10,7 @@ list(APPEND programs
   sonde_bufr2ioda.py
   sonde_tac2ioda.py
   synop_bufr2ioda.py
+  tcv_txt2ioda.py
 )
 
 file( RELATIVE_PATH SCRIPT_LIB_PATH ${CMAKE_BINARY_DIR}/bin ${PYIODACONV_BUILD_LIBDIR} )

src/conventional/tcv_txt2ioda.py

@@ -0,0 +1,147 @@
+#!/usr/bin/env python3
+
+# converter for tcvitals format (tropical storm central pressure) data
+
+import numpy as np
+import netCDF4 as nc
+import argparse
+import pandas as pd
+from datetime import datetime
+import dateutil.parser
+from builtins import str
+from pyioda import ioda_obs_space as ioda_ospace
+
+
+def read_file(file_name):
+
+    names = ["center", "storm_id", "storm_name", "datestr",
+             "lat", "latns", "lon", "lonew", "stdir", "stspd", "pcen"]
+
+#                ctr      id     name    datestr     lat      latns
+    colspecs = [(0, 4), (5, 8), (9, 18), (19, 32), (33, 36), (36, 37),
+                (38, 42), (42, 43), (44, 47), (48, 51), (52, 56)]
+#                   lon      lonew      dir       spd       pres
+
+    data = pd.read_fwf(file_name, colspecs=colspecs, header=None, names=names)
+
+    return data
+
+
+def create_output(data, output_name):
+
+    float_missing_value = nc.default_fillvals['f4']
+
+    lon = data.lon / 10.0
+    lon[data.lonew == 'W'] = 360. - lon
+    lat = data.lat / 10.0
+    lat[data.latns == 'S'] = - lat
+# calculate 'obs error' as in GSI
+    tcp_refps = 1000.0
+    tcp_width = 50.0
+    tcp_ermin = 0.75
+    tcp_ermax = 5.0
+#      alpha=max(min(psdif/tcp_width,one),zero)
+    psdif = tcp_refps - data.pcen
+    alpha = psdif / tcp_width
+    alpha[alpha > 1.0] = 1.0
+    alpha[alpha < 0.0] = 0.0
+    oberr = tcp_ermin+(tcp_ermax-tcp_ermin)*alpha
+    nobs = len(lat)
+    height = np.zeros(nobs)
+    obstype = np.full(nobs, 112)
+    obssubtype = np.full(nobs, 0)
+    preqc = np.full(nobs, 2)
+    preuseflg = np.full(nobs, 1)
+    tempK = np.full(nobs, float_missing_value)
+    psminPa = data.pcen*100.
+    oberrPa = oberr*100.
+    epoch = datetime.fromisoformat("1970-01-01T00:00:00Z")
+    obdate = pd.to_datetime(data.datestr, format="%Y%m%d %H%M")
+    eparr = np.full_like(obdate, epoch)
+    obdiff = obdate - eparr
+    sec_since_epoch = obdiff.astype('int64') // 10**9
+    station_id = data['center'].str.cat(data['storm_id'], sep='_')
+    dims = {
+        'Location': np.arange(0, psminPa.shape[0]),
+    }
+    obsspace = ioda_ospace.ObsSpace(output_name, mode='w', dim_dict=dims)
+    obsspace.create_var('MetaData/dateTime', dtype='int64') \
+        .write_attr('units', 'seconds since 1970-01-01T00:00:00Z') \
+        .write_attr('long_name', 'dateTime') \
+        .write_data(sec_since_epoch)
+    obsspace.create_var('MetaData/stationElevation', dtype=np.float32) \
+        .write_attr('units', 'm') \
+        .write_attr('long_name', 'Height Of Station') \
+        .write_data(height)
+    obsspace.create_var('MetaData/height', dtype=np.float32) \
+        .write_attr('units', 'm') \
+        .write_attr('long_name', 'Height Of Station') \
+        .write_data(height)
+    obsspace.create_var('MetaData/latitude', dtype=np.float32) \
+        .write_attr('units', 'degrees_north') \
+        .write_attr('valid_range', np.array([-90, 90], dtype=np.float32)) \
+        .write_attr('long_name', 'Latitude') \
+        .write_data(lat)
+    obsspace.create_var('MetaData/longitude', dtype=np.float32) \
+        .write_attr('units', 'degrees_east') \
+        .write_attr('valid_range', np.array([0, 360], dtype=np.float32)) \
+        .write_attr('long_name', 'Longitude') \
+        .write_data(lon)
+    obsspace.create_var('MetaData/pressure', dtype=np.float32) \
+        .write_attr('units', 'Pa') \
+        .write_attr('valid_range', np.array([105000, 80000], dtype=np.float32)) \
+        .write_attr('long_name', 'Pressure') \
+        .write_data(psminPa)
+    obsspace.create_var('MetaData/stationIdentification', dtype=station_id.dtype) \
+        .write_attr('long_name', 'stationIdentification') \
+        .write_data(station_id)
+    obsspace.create_var('ObsType/stationPressure', dtype=np.int32) \
+        .write_attr('long_name', 'Station Pressure Observation Type') \
+        .write_data(obstype)
+
+    obsspace.create_var('ObsSubType/stationPressure', dtype=np.int32) \
+        .write_attr('long_name', 'Station Pressure Observation subType') \
+        .write_data(obssubtype)
+
+# ObsError: initial error values loaded from the input ioda file
+    obsspace.create_var('ObsError/stationPressure', dtype=np.float32) \
+        .write_attr('units', 'Pa') \
+        .write_attr('long_name', 'ObsError') \
+        .write_attr('coordinates', 'longitude latitude') \
+        .write_data(oberrPa)
+
+    obsspace.create_var('ObsValue/stationPressure', dtype=np.float32) \
+        .write_attr('units', 'Pa') \
+        .write_attr('valid_range', np.array([105000, 80000], dtype=np.float32)) \
+        .write_attr('long_name', 'ObsValue') \
+        .write_attr('coordinates', 'longitude latitude') \
+        .write_data(psminPa)
+
+# temporary fake array for temperature requested by SfcCorrected
+    obsspace.create_var('ObsValue/airTemperature', dtype=np.float32) \
+        .write_attr('units', 'K') \
+        .write_attr('valid_range', np.array([250, 380], dtype=np.float32)) \
+        .write_attr('long_name', 'ObsValue') \
+        .write_attr('coordinates', 'longitude latitude') \
+        .write_data(tempK)
+
+
+def main():
+
+    desc = "Reads a tcvitals (tropical storm center data) file and converts into IODA format"
+    parser = argparse.ArgumentParser(description=desc)
+    parser.add_argument(
+        '-i', '--input', help='Input tcvitals text file',
+        type=str, required=True, default=None)
+    parser.add_argument(
+        '-o', '--output', help='Output tcvitals ioda nc4 file',
+        type=str, required=True, default=None)
+    args = parser.parse_args()
+
+    data = read_file(args.input)
+
+    create_output(data, args.output)
+
+
+if __name__ == '__main__':
+    main()

src/gsi_ncdiag/gsi_ncdiag.py

@@ -53,6 +53,9 @@
     ],
     "conv_sst": [
         'sst',
+    ],
+    "conv_tcp": [
+        'tcp',
     ]
 }
 
@@ -81,6 +84,7 @@
             722, 723, 740, 741, 742, 743, 744, 745, \
             750, 751, 752, 753, 754, 755, 786, 803, 804, 820, 821, 825],
     "sst": [181, 182, 183, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202],
+    "tcp": [112],
     # 132 are dropsondes
 }
 
@@ -195,6 +199,7 @@
     "bend": ["bendingAngle"],
     "refract": ["atmosphericRefractivity"],
     "sst": ["seaSurfaceTemperature"],
+    "tcp": ["stationPressure"],
 }
 
 conv_gsivarnames = {
@@ -206,6 +211,7 @@
     "bend": ["Observation"],
     "refract": ["Observation"],
     "sst": ["Observation"],
+    "tcp": ["Observation"],
 }
 
 gsi_add_vars_allsky = {
@@ -782,6 +788,9 @@ def toGeovals(self, OutDir, clobber=True):
                 if (v == 'sst'):
                     outname = OutDir + '/' + v + '_geoval_' + \
                         self.validtime.strftime("%Y%m%d%H") + '.nc4'
+                if (v == 'tcp'):
+                    outname = OutDir + '/' + v + '_geoval_' + \
+                        self.validtime.strftime("%Y%m%d%H") + '.nc4'
                 if (p == 'windprof' or p == 'satwind' or p == 'scatwind' or p == 'vadwind' or p == 'pibal'):
                     outname = OutDir + '/' + p + '_geoval_' + \
                         self.validtime.strftime("%Y%m%d%H") + '.nc4'
@@ -905,6 +914,9 @@ def toIODAobs(self, OutDir, TotalBias=False, clobber=True, platforms=None):
                 if (v == 'sst'):
                     outname = OutDir + '/' + v + '_obs_' + \
                         self.validtime.strftime("%Y%m%d%H") + '.nc4'
+                if (v == 'tcp'):
+                    outname = OutDir + '/' + v + '_obs_' + \
+                        self.validtime.strftime("%Y%m%d%H") + '.nc4'
                 if (p == 'windprof' or p == 'satwind' or p == 'scatwind' or p == 'vadwind' or p == 'pibal'):
                     outname = OutDir + '/' + p + '_obs_' + \
                         self.validtime.strftime("%Y%m%d%H") + '.nc4'

test/CMakeLists.txt

@@ -112,6 +112,7 @@ list( APPEND test_input
   testinput/balloonsonde_sample_20241001T00Z_PT15M.json
   testinput/tec_groundBased_TENET_20201001T00Z.tec
   testinput/COWVR_TSDR.020421.20240430T215845.20240430T230345.V1001.sample.h5
+  testinput/gfs.231010.t00z.syndata.tcvitals.txt
 )
 
 list( APPEND test_output
@@ -215,6 +216,7 @@ list( APPEND test_output
   testoutput/obs.20241001T000000Z_PT15M_balloonsonde.nc4
   testoutput/obs.20201001T000000Z_PT1M_tec_groundbased.nc4
   testoutput/obs.20240430T220600Z_PT2M_cowvr_iss.nc4
+  testoutput/tcp_obs_2023101000.nc4
 )
 
 if( iodaconv_gnssro_ENABLED )
@@ -873,15 +875,15 @@ ecbuild_add_test( TARGET  test_${PROJECT_NAME}_pace_radiance_L1B
                           -d 2019032112"
                           pace_radiance_L1B.nc ${IODA_CONV_COMP_TOL})
 
-ecbuild_add_test( TARGET  test_${PROJECT_NAME}_modis_water_leaving_radiance_L2                                               
-                  TYPE    SCRIPT                                                                               
-                  ENVIRONMENT "PYTHONPATH=${IODACONV_PYTHONPATH}"                                              
-                  COMMAND bash                                                                                 
-                  ARGS    ${CMAKE_BINARY_DIR}/bin/iodaconv_comp.sh                                             
-                          netcdf                                                                               
-                          "${Python3_EXECUTABLE} ${CMAKE_BINARY_DIR}/bin/oc_l2_radiance2ioda.py                 
-                          -i testinput/modis_aqua_Rrs_OC_L2.nc                                                    
-                          -o testrun/modis_water_leaving_radiance_L2.nc                                                      
+ecbuild_add_test( TARGET  test_${PROJECT_NAME}_modis_water_leaving_radiance_L2
+                  TYPE    SCRIPT
+                  ENVIRONMENT "PYTHONPATH=${IODACONV_PYTHONPATH}"
+                  COMMAND bash
+                  ARGS    ${CMAKE_BINARY_DIR}/bin/iodaconv_comp.sh
+                          netcdf
+                          "${Python3_EXECUTABLE} ${CMAKE_BINARY_DIR}/bin/oc_l2_radiance2ioda.py
+                          -i testinput/modis_aqua_Rrs_OC_L2.nc
+                          -o testrun/modis_water_leaving_radiance_L2.nc
                           -d 2021080112"
                           modis_water_leaving_radiance_L2.nc ${IODA_CONV_COMP_TOL})
 
@@ -979,6 +981,21 @@ ecbuild_add_test( TARGET  test_${PROJECT_NAME}_sondetac
                           --date 2021-08-16T12:00:00Z --netcdf"
                           2021081612_sonde_small.nc ${IODA_CONV_COMP_TOL})
 
+#===============================================================================
+# Conventional data, tropical storm central pressure
+#===============================================================================
+
+ecbuild_add_test( TARGET  test_${PROJECT_NAME}_tcvitals
+                  TYPE    SCRIPT
+                  ENVIRONMENT "PYTHONPATH=${IODACONV_PYTHONPATH}"
+                  COMMAND bash
+                  ARGS    ${CMAKE_BINARY_DIR}/bin/iodaconv_comp.sh
+                          netcdf
+                          "${Python3_EXECUTABLE} ${CMAKE_BINARY_DIR}/bin/tcv_txt2ioda.py
+                          -i testinput/gfs.231010.t00z.syndata.tcvitals.txt
+                          -o testrun/tcp_obs_2023101000.nc4"
+                          tcp_obs_2023101000.nc4 ${IODA_CONV_COMP_TOL})
+
 #===============================================================================
 # MRMS (multi-radar, multi-sensor) radar reflectivity CONUS converter
 #===============================================================================
@@ -1021,7 +1038,7 @@ ecbuild_add_test( TARGET  test_iodaconv_tropics
                   -o testrun/20220216T00Z_PT3H_tms_tropics-01.nc4
                   -d 2022021600"
                   20220216T00Z_PT3H_tms_tropics-01.nc4 ${IODA_CONV_COMP_TOL})
-                  
+
 ecbuild_add_test( TARGET  test_iodaconv_amsr2_gcom-w1
                   TYPE    SCRIPT
                   ENVIRONMENT "PYTHONPATH=${IODACONV_PYTHONPATH}"
@@ -1033,7 +1050,7 @@ ecbuild_add_test( TARGET  test_iodaconv_amsr2_gcom-w1
                   -o testrun/20220216T00Z_PT1H_amsr2_sample.nc4
                   -d 2022021600"
                   20220216T00Z_PT1H_amsr2_sample.nc4 ${IODA_CONV_COMP_TOL})
-                  
+
 ecbuild_add_test( TARGET  test_iodaconv_gmi_gpm
                   TYPE    SCRIPT
                   ENVIRONMENT "PYTHONPATH=${IODACONV_PYTHONPATH}"
@@ -1479,7 +1496,7 @@ ecbuild_add_test( TARGET  test_${PROJECT_NAME}_modis_aod
                           -i testinput/modis_aod.hdf
                           --platform Terra
 			  --date_range 2021080100 2021080101
-                          -o testrun/modis_aod.nc"			
+                          -o testrun/modis_aod.nc"
                           modis_aod.nc ${IODA_CONV_COMP_TOL})
 
 ecbuild_add_test( TARGET  test_${PROJECT_NAME}_aeronet_aod
@@ -2307,8 +2324,8 @@ endif()
                     COMMAND bash
                     ARGS    ${CMAKE_BINARY_DIR}/bin/iodaconv_comp.sh
                             netcdf
-                            "${Python3_EXECUTABLE} ${CMAKE_BINARY_DIR}/bin/cowvr_hdf5_2ioda.py 
-                            -i testinput/COWVR_TSDR.020421.20240430T215845.20240430T230345.V1001.sample.h5 
+                            "${Python3_EXECUTABLE} ${CMAKE_BINARY_DIR}/bin/cowvr_hdf5_2ioda.py
+                            -i testinput/COWVR_TSDR.020421.20240430T215845.20240430T230345.V1001.sample.h5
                             -o testrun/obs.20240430T220600Z_PT2M_cowvr_iss.nc4
                             -d 2024050100"
                             obs.20240430T220600Z_PT2M_cowvr_iss.nc4 ${IODA_CONV_COMP_TOL})

test/testinput/gfs.231010.t00z.syndata.tcvitals.txt

@@ -0,0 +1,4 @@
+NHC  15E LIDIA     20231010 0000 182N 1106W 065 041 0983 1005 0408 33 037 0130 0130 0093 0093 D
+NHC  16E MAX       20231010 0000 179N 1011W 020 026 1002 1005 0111 18 037 0056 0056 0056 -999 S
+JTWC 15W BOLAVEN   20231010 0000 137N 1465E 310 067 0990 1006 0509 36 027 0213 0222 0158 0139 D
+JTWC 14W KOINU     20231010 0000 208N 1103E 250 047 1011 1013 0111 13 018 -999 -999 -999 -999 S