diff --git a/src/compo/viirs_l1bnc2ioda.py b/src/compo/viirs_l1bnc2ioda.py index b6dac3ad7..e20356ce9 100755 --- a/src/compo/viirs_l1bnc2ioda.py +++ b/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 diff --git a/src/hdf5/amsr2_2ioda.py b/src/hdf5/amsr2_2ioda.py index 8935258f8..7e6f735d7 100755 --- a/src/hdf5/amsr2_2ioda.py +++ b/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) diff --git a/src/hdf5/cowvr_hdf5_2ioda.py b/src/hdf5/cowvr_hdf5_2ioda.py index c18208f65..d4571a796 100755 --- a/src/hdf5/cowvr_hdf5_2ioda.py +++ b/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) diff --git a/src/hdf5/gmi_2ioda.py b/src/hdf5/gmi_2ioda.py index 9892486d9..406bce6a1 100755 --- a/src/hdf5/gmi_2ioda.py +++ b/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) diff --git a/src/hdf5/ssmis_upp_2ioda.py b/src/hdf5/ssmis_upp_2ioda.py index 545a49efc..69031aa80 100755 --- a/src/hdf5/ssmis_upp_2ioda.py +++ b/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)]]) diff --git a/src/hdf5/tropics_2ioda.py b/src/hdf5/tropics_2ioda.py index 1eb0229d6..9a1358f0c 100755 --- a/src/hdf5/tropics_2ioda.py +++ b/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') diff --git a/src/pyiodaconv/def_jedi_utils.py b/src/pyiodaconv/def_jedi_utils.py index d8a9cfca5..8b7fc55d1 100644 --- a/src/pyiodaconv/def_jedi_utils.py +++ b/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