Fix function compute_scan_angle (QC flags), remove unused argument

src/compo/viirs_l1bnc2ioda.py

@@ -157,7 +157,7 @@ def _read(self):
             solar_aa = geo_ncd.groups['geolocation_data'].variables['solar_azimuth'][:].data.ravel()
             sensor_za = geo_ncd.groups['geolocation_data'].variables['sensor_zenith'][:].data.ravel()
             sensor_aa = geo_ncd.groups['geolocation_data'].variables['sensor_azimuth'][:].data.ravel()
-            sensor_va = compute_scan_angle(sensor_za, np.full_like(sensor_za, orbit_height), sensor_za)
+            sensor_va = compute_scan_angle(np.full_like(sensor_za, orbit_height), sensor_za)
             landwat_mask = geo_ncd.groups['geolocation_data'].variables['land_water_mask'][:].data.ravel()
 
             nlocs = lons.size

src/hdf5/amsr2_2ioda.py

@@ -170,7 +170,6 @@ def get_data(f, obs_data):
     sat_altitude = np.empty_like(instr_scan_ang)
     sat_altitude[:] = f.attrs['SatelliteAltitude'].item().strip('km')
     obs_data[('sensorViewAngle', metaDataName)] = compute_scan_angle(
-        instr_scan_ang,
         sat_altitude,
         instr_scan_ang)
 

src/hdf5/cowvr_hdf5_2ioda.py

@@ -179,7 +179,6 @@ def get_tempest_data(f, obs_data, add_qc=True):
     obs_data[('sensorZenithAngle', metaDataName)] = np.array(f['Geolocation']['earth_inc_ang'], dtype='float32')
     obs_data[('sensorAzimuthAngle', metaDataName)] = np.array(f['Geolocation']['earth_az_ang'], dtype='float32')
     obs_data[('sensorViewAngle', metaDataName)] = compute_scan_angle(
-        np.array(f['Geolocation']['instr_scan_ang'], dtype='float32'),
         sensor_altitude,
         np.array(f['Geolocation']['earth_inc_ang'], dtype='float32'),
         qc_flag=qc_flag)
@@ -246,7 +245,6 @@ def get_cowvr_data(f, obs_data, add_qc=True):
     obs_data[('sensorZenithAngle', metaDataName)] = np.array(f['GeolocationAndFlags']['earth_inc_ang'], dtype='float32')
     obs_data[('sensorAzimuthAngle', metaDataName)] = np.array(f['GeolocationAndFlags']['earth_az_ang'], dtype='float32')
     obs_data[('sensorViewAngle', metaDataName)] = compute_scan_angle(
-        np.array(f['GeolocationAndFlags']['instr_scan_ang'], dtype='float32'),
         sensor_altitude,
         np.array(f['GeolocationAndFlags']['earth_inc_ang'], dtype='float32'),
         qc_flag=sat_alt_flag)

src/hdf5/gmi_2ioda.py

@@ -155,7 +155,6 @@ def get_data(f, obs_data):
     sat_altitude = np.empty_like(instr_scan_ang)
     sat_altitude[:] = 407.0
     obs_data[('sensorViewAngle', metaDataName)] = compute_scan_angle(
-        instr_scan_ang,
         sat_altitude,
         instr_scan_ang)
 

src/hdf5/ssmis_upp_2ioda.py

@@ -213,7 +213,6 @@ def populate_obsValue(line, local_data, WMO_sat_ID=int_missing_value, ssmis_uas=
     # 60 km                          52.62
     local_data[('sensorViewAngle', metaDataName)].append(sensor_zenith)
 #   local_data[('sensorViewAngle', metaDataName)] = compute_scan_angle(
-#       sensor_zenith,
 #       sensor_altitude,
 #       sensor_zenith,
 #       qc_flag=[[int(irej)]])

src/hdf5/tropics_2ioda.py

@@ -331,7 +331,6 @@ def get_data_deprecated(f, obs_data, skip=1):
     sat_altitude = np.empty_like(instr_scan_ang)
     sat_altitude[:] = 550.
     obs_data[('sensorViewAngle', metaDataName)] = compute_scan_angle(
-        instr_scan_ang,
         sat_altitude,
         instr_scan_ang)
     obs_data[('dateTime', metaDataName)] = np.array(get_epoch_time(f), dtype='int64')

src/pyiodaconv/def_jedi_utils.py

@@ -22,6 +22,7 @@
 epoch = datetime(1970, 1, 1, tzinfo=timezone.utc)
 ioda_float_type = 'float32'
 ioda_int_type = 'int32'
+double_missing_value = iconv.get_default_fill_val(np.float64)
 float_missing_value = iconv.get_default_fill_val(np.float32)
 int_missing_value = iconv.get_default_fill_val(np.int32)
 long_missing_value = iconv.get_default_fill_val(np.int64)
@@ -73,7 +74,7 @@ def set_obspace_attributes(VarAttrs):
     return VarAttrs
 
 
-def compute_scan_angle(instr_scan_ang, sensor_altitude, sensor_zenith, qc_flag=[None]):
+def compute_scan_angle(sensor_altitude, sensor_zenith, qc_flag=None):
 
     # should come from standard table
     earth_mean_radius_km = 6378.1370  # WGS84
@@ -86,18 +87,23 @@ def compute_scan_angle(instr_scan_ang, sensor_altitude, sensor_zenith, qc_flag=[
     d2r = np.pi/180.
     r2d = 180./np.pi
 
-    # do we need a missing here
-    ratio = np.empty_like(sensor_altitude)
+    # Initialize output with missing values
+    scanang = np.full_like(sensor_altitude, float_missing_value, dtype=float)
 
-    # compute scan angle
-    if not qc_flag[0]:
-        qc_flag = np.zeros_like(sensor_altitude)
-    good = qc_flag[:] == 0
-    if sum(good) > 0:
-        ratio[good] = earth_mean_radius_km/(earth_mean_radius_km + sensor_altitude[good]/1000.)
-
-    # γ = arcsin(R / (R + h) * sin(theta)),h: sat alt; theta: sat zenith angle
-    scanang = np.arcsin(ratio*np.sin(abs(sensor_zenith)*d2r))*r2d
+    # Handle qc_flag
+    if qc_flag is None:
+        good = np.ones_like(sensor_altitude, dtype=bool)
+    else:
+        qc_flag = np.asarray(qc_flag)
+        good = (qc_flag == 0)
+
+    # compute scan angle only for good observations
+    if np.any(good):
+        # γ = arcsin(R / (R + h) * sin(theta)),h: sat alt; theta: sat zenith angle
+        ratio = earth_mean_radius_km/(earth_mean_radius_km + sensor_altitude[good]/1000.)
+        sin_angle = ratio * np.sin(np.abs(sensor_zenith[good]) * d2r)
+        sin_angle = np.clip(sin_angle, -1.0, 1.0)   # avoid arcsin domain error
+        scanang[good] = np.arcsin(sin_angle)*r2d
 
     return scanang