From 2985fe0b92290722b682423ac48f019b69792b82 Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Tue, 3 Dec 2019 18:52:17 +0100 Subject: [PATCH 01/12] Add LAC support --- pygac/gac_klm.py | 527 ++------------------- pygac/gac_pod.py | 408 +---------------- pygac/gac_reader.py | 807 +------------------------------- pygac/klm_reader.py | 732 +++++++++++++++++++++++++++++ pygac/lac_klm.py | 227 +++++++++ pygac/lac_pod.py | 93 ++++ pygac/lac_reader.py | 40 ++ pygac/pod_reader.py | 491 ++++++++++++++++++++ pygac/pygac_geotiepoints.py | 16 +- pygac/reader.py | 891 ++++++++++++++++++++++++++++++++++++ pygac/tests/test_pod.py | 11 +- pygac/tests/test_reader.py | 6 +- 12 files changed, 2572 insertions(+), 1677 deletions(-) create mode 100644 pygac/klm_reader.py create mode 100644 pygac/lac_klm.py create mode 100644 pygac/lac_pod.py create mode 100644 pygac/lac_reader.py create mode 100644 pygac/pod_reader.py create mode 100644 pygac/reader.py diff --git a/pygac/gac_klm.py b/pygac/gac_klm.py index d6702abb..47b49e9b 100644 --- a/pygac/gac_klm.py +++ b/pygac/gac_klm.py @@ -23,6 +23,7 @@ # along with this program. If not, see . """Read a gac file. + Reads L1b GAC data from KLM series of satellites (NOAA-15 and later) and does most of the computations. Format specification can be found here: http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c8/sec83142-1.htm @@ -31,325 +32,14 @@ from __future__ import print_function -import numpy as np -from pygac.correct_tsm_issue import TSM_AFFECTED_INTERVALS_KLM, get_tsm_idx -from pygac.gac_reader import GACReader, inherit_doc -import datetime -from pygac import gac_io import logging -LOG = logging.getLogger(__name__) - -# GAC header object +import numpy as np -header = np.dtype([("data_set_creation_site_id", "S3"), - ("ascii_blank_=_x20", "S1"), - ("noaa_level_1b_format_version_number:", ">u2"), - ("noaa_level_1b_format_version_year", ">u2"), - ("noaa_level_1b_format_version_day_of_year", ">u2"), - ("reserved_for_logical_record_length", ">u2"), - ("reserved_for_block_size", ">u2"), - ("count_of_header_records", ">u2"), - ("zero_fill0", ">i2", (3, )), - ("data_set_name", "S42"), - ("processing_block_identification", "S8"), - ("noaa_spacecraft_identification_code", ">u2"), - ("instrument_id", ">u2"), - ("data_type_code", ">u2"), - ("tip_source_code", ">u2"), - ("start_of_data_set_day_count_starting_from_0_at_00h,_1_jan_1950", - ">u4"), - ("start_of_data_set_year", ">u2"), - ("start_of_data_set_day_of_year", ">u2"), - ("start_of_data_set_utc_time_of_day", ">u4"), - ("end_of_data_set_day_count_starting_from_0_at_00h,_1_jan_1950", - ">u4"), - ("end_of_data_set_year", ">u2"), - ("end_of_data_set_day_of_year", ">u2"), - ("end_of_data_set_utc_time_of_day", ">u4"), - ("year_of_last_cpids_update", ">u2"), - ("day_of_year_of_last_cpids_update", ">u2"), - ("zero_fill1", ">i2", (4, )), - # DATA SET QUALITY INDICATORS - ("instrument_status", ">u4"), - ("zero_fill2", ">i2"), - ("record_number_of_status_change", ">u2"), - ("second_instrument_status", ">u4"), - ("count_of_data_records", ">u2"), - ("count_of_calibrated,_earth_located_scan_lines", ">u2"), - ("count_of_missing_scan_lines", ">u2"), - ("count_of_data_gaps", ">u2"), - ("filler0", ">u2"), - ("count_of_pacs_detected_tip_parity_errors", ">u2"), - ("filler1", ">u2"), - ("time_sequence_error", ">u2"), - ("time_sequence_error_code", ">u2"), - ("socc_clock_update_indicator", ">u2"), - ("earth_location_error_indicator", ">u2"), - ("earth_location_error_code", ">u2"), - ("pacs_status_bit_field", ">u2"), - ("data_source", ">u2"), - ("zero_fill3", ">i4"), - ("reserved_for_the_ingester", "S8"), - ("reserved_for_decommutation", "S8"), - ("zero_fill4", ">i2", (5, )), - # CALIBRATION - ("ramp_auto_calibration_indicators_bit_field", ">u2"), - ("year_of_most_recent_solar_channel_calibration", ">u2"), - ("day_of_year_of_most_recent_solar_channel_calibration", - ">u2"), - ("primary_calibration_algorithm_id", ">u2"), - ("primary_calibration_algorithm_selected_options", ">u2"), - ("secondary_calibration_algorithm_id", ">u2"), - ("secondary_calibration_algorithm_selected_options", ">u2"), - ("ir_target_temperature_1_conversion_coefficient_1", ">i2"), - ("ir_target_temperature_1_conversion_coefficient_2", ">i2"), - ("ir_target_temperature_1_conversion_coefficient_3", ">i2"), - ("ir_target_temperature_1_conversion_coefficient_4", ">i2"), - ("ir_target_temperature_1_conversion_coefficient_5", ">i2"), - ("ir_target_temperature_1_conversion_coefficient_6", ">i2"), - ("ir_target_temperature_2_conversion_coefficient_1", ">i2"), - ("ir_target_temperature_2_conversion_coefficient_2", ">i2"), - ("ir_target_temperature_2_conversion_coefficient_3", ">i2"), - ("ir_target_temperature_2_conversion_coefficient_4", ">i2"), - ("ir_target_temperature_2_conversion_coefficient_5", ">i2"), - ("ir_target_temperature_2_conversion_coefficient_6", ">i2"), - ("ir_target_temperature_3_conversion_coefficient_1", ">i2"), - ("ir_target_temperature_3_conversion_coefficient_2", ">i2"), - ("ir_target_temperature_3_conversion_coefficient_3", ">i2"), - ("ir_target_temperature_3_conversion_coefficient_4", ">i2"), - ("ir_target_temperature_3_conversion_coefficient_5", ">i2"), - ("ir_target_temperature_3_conversion_coefficient_6", ">i2"), - ("ir_target_temperature_4_conversion_coefficient_1", ">i2"), - ("ir_target_temperature_4_conversion_coefficient_2", ">i2"), - ("ir_target_temperature_4_conversion_coefficient_3", ">i2"), - ("ir_target_temperature_4_conversion_coefficient_4", ">i2"), - ("ir_target_temperature_4_conversion_coefficient_5", ">i2"), - ("ir_target_temperature_4_conversion_coefficient_6", ">i2"), - ("zero_fill5", ">i4", (2, )), - # RADIANCE CONVERSION - ("ch_1_solar_filtered_irradiance_in_wavelength", ">i4"), - ("ch_1_equivalent_filter_width_in_wavelength", ">i4"), - ("ch_2_solar_filtered_irradiance_in_wavelength", ">i4"), - ("ch_2_equivalent_filter_width_in_wavelength", ">i4"), - ("ch_3a_solar_filtered_irradiance_in_wavelength", ">i4"), - ("ch_3a_equivalent_filter_width_in_wavelength", ">i4"), - ("ch_3b_central_wavenumber_(see", ">i4"), - ("ch_3b_constant_1", ">i4"), - ("ch_3b_constant_2", ">i4"), - ("ch_4_central_wavenumber", ">i4"), - ("ch_4_constant_1", ">i4"), - ("ch_4_constant_2", ">i4"), - ("ch_5_central_wavenumber", ">i4"), - ("ch_5_constant_1", ">i4"), - ("ch_5_constant_2", ">i4"), - ("zero_fill6", ">i4", (3, )), - # NAVIGATION - ("reference_ellipsoid_model_id", "S8"), - ("nadir_earth_location_tolerance", ">u2"), - ("earth_location_bit_field", ">u2"), - ("zero_fill7", ">i2"), - ("constant_roll_attitude_error", ">i2"), - ("constant_pitch_attitude_error", ">i2"), - ("constant_yaw_attitude_error", ">i2"), - ("epoch_year_for_orbit_vector", ">u2"), - ("day_of_epoch_year_for_orbit_vector", ">u2"), - ("epoch_utc_time_of_day_for_orbit_vector", ">u4"), - ("semi-major_axis", ">i4"), - ("eccentricity", ">i4"), - ("inclination", ">i4"), - ("argument_of_perigee", ">i4"), - ("right_ascension_of_the_ascending_node", ">i4"), - ("mean_anomaly", ">i4"), - ("position_vector_x_component", ">i4"), - ("position_vector_y_component", ">i4"), - ("position_vector_z_component", ">i4"), - ("velocity_vector_x-dot_component", ">i4"), - ("velocity_vector_y-dot_component", ">i4"), - ("velocity_vector_z-dot_component", ">i4"), - ("earth_sun_distance_ratio", ">u4"), - ("zero_fill8", ">i4", (4, )), - # ANALOG TELEMETRY CONVERSION - ("patch_temperature_conversion_coefficient_1", ">i2"), - ("patch_temperature_conversion_coefficient_2", ">i2"), - ("patch_temperature_conversion_coefficient_3", ">i2"), - ("patch_temperature_conversion_coefficient_4", ">i2"), - ("patch_temperature_conversion_coefficient_5", ">i2"), - ("reserved0", ">i2"), - ("patch_temperature_extended_conversion_coefficient_1", - ">i2"), - ("patch_temperature_extended_conversion_coefficient_2", - ">i2"), - ("patch_temperature_extended_conversion_coefficient_3", - ">i2"), - ("patch_temperature_extended_conversion_coefficient_4", - ">i2"), - ("patch_temperature_extended_conversion_coefficient_5", - ">i2"), - ("reserved1", ">i2"), - ("patch_power_conversion_coefficient_1", ">i2"), - ("patch_power_conversion_coefficient_2", ">i2"), - ("patch_power_conversion_coefficient_3", ">i2"), - ("patch_power_conversion_coefficient_4", ">i2"), - ("patch_power_conversion_coefficient_5", ">i2"), - ("reserved2", ">i2"), - ("radiator_temperature_conversion_coefficient_1", ">i2"), - ("radiator_temperature_conversion_coefficient_2", ">i2"), - ("radiator_temperature_conversion_coefficient_3", ">i2"), - ("radiator_temperature_conversion_coefficient_4", ">i2"), - ("radiator_temperature_conversion_coefficient_5", ">i2"), - ("reserved3", ">i2"), - ("blackbody_temperature_1_conversion_coefficient_1", ">i2"), - ("blackbody_temperature_1_conversion_coefficient_2", ">i2"), - ("blackbody_temperature_1_conversion_coefficient_3", ">i2"), - ("blackbody_temperature_1_conversion_coefficient_4", ">i2"), - ("blackbody_temperature_1_conversion_coefficient_5", ">i2"), - ("reserved4", ">i2"), - ("blackbody_temperature_2_conversion_coefficient_1", ">i2"), - ("blackbody_temperature_2_conversion_coefficient_2", ">i2"), - ("blackbody_temperature_2_conversion_coefficient_3", ">i2"), - ("blackbody_temperature_2_conversion_coefficient_4", ">i2"), - ("blackbody_temperature_2_conversion_coefficient_5", ">i2"), - ("reserved5", ">i2"), - ("blackbody_temperature_3_conversion_coefficient_1", ">i2"), - ("blackbody_temperature_3_conversion_coefficient_2", ">i2"), - ("blackbody_temperature_3_conversion_coefficient_3", ">i2"), - ("blackbody_temperature_3_conversion_coefficient_4", ">i2"), - ("blackbody_temperature_3_conversion_coefficient_5", ">i2"), - ("reserved6", ">i2"), - ("blackbody_temperature_4_conversion_coefficient_1", ">i2"), - ("blackbody_temperature_4_conversion_coefficient_2", ">i2"), - ("blackbody_temperature_4_conversion_coefficient_3", ">i2"), - ("blackbody_temperature_4_conversion_coefficient_4", ">i2"), - ("blackbody_temperature_4_conversion_coefficient_5", ">i2"), - ("reserved7", ">i2"), - ("electronics_current_conversion_coefficient_1", ">i2"), - ("electronics_current_conversion_coefficient_2", ">i2"), - ("electronics_current_conversion_coefficient_3", ">i2"), - ("electronics_current_conversion_coefficient_4", ">i2"), - ("electronics_current_conversion_coefficient_5", ">i2"), - ("reserved8", ">i2"), - ("motor_current_conversion_coefficient_1", ">i2"), - ("motor_current_conversion_coefficient_2", ">i2"), - ("motor_current_conversion_coefficient_3", ">i2"), - ("motor_current_conversion_coefficient_4", ">i2"), - ("motor_current_conversion_coefficient_5", ">i2"), - ("reserved9", ">i2"), - ("earth_shield_position_conversion_coefficient_1", ">i2"), - ("earth_shield_position_conversion_coefficient_2", ">i2"), - ("earth_shield_position_conversion_coefficient_3", ">i2"), - ("earth_shield_position_conversion_coefficient_4", ">i2"), - ("earth_shield_position_conversion_coefficient_5", ">i2"), - ("reserved10", ">i2"), - ("electronics_temperature_conversion_coefficient_1", ">i2"), - ("electronics_temperature_conversion_coefficient_2", ">i2"), - ("electronics_temperature_conversion_coefficient_3", ">i2"), - ("electronics_temperature_conversion_coefficient_4", ">i2"), - ("electronics_temperature_conversion_coefficient_5", ">i2"), - ("reserved11", ">i2"), - ("cooler_housing_temperature_conversion_coefficient_1", - ">i2"), - ("cooler_housing_temperature_conversion_coefficient_2", - ">i2"), - ("cooler_housing_temperature_conversion_coefficient_3", - ">i2"), - ("cooler_housing_temperature_conversion_coefficient_4", - ">i2"), - ("cooler_housing_temperature_conversion_coefficient_5", - ">i2"), - ("reserved12", ">i2"), - ("baseplate_temperature_conversion_coefficient_1", ">i2"), - ("baseplate_temperature_conversion_coefficient_2", ">i2"), - ("baseplate_temperature_conversion_coefficient_3", ">i2"), - ("baseplate_temperature_conversion_coefficient_4", ">i2"), - ("baseplate_temperature_conversion_coefficient_5", ">i2"), - ("reserved13", ">i2"), - ("motor_housing_temperature_conversion_coefficient_1", - ">i2"), - ("motor_housing_temperature_conversion_coefficient_2", - ">i2"), - ("motor_housing_temperature_conversion_coefficient_3", - ">i2"), - ("motor_housing_temperature_conversion_coefficient_4", - ">i2"), - ("motor_housing_temperature_conversion_coefficient_5", - ">i2"), - ("reserved14", ">i2"), - ("a/d_converter_temperature_conversion_coefficient_1", - ">i2"), - ("a/d_converter_temperature_conversion_coefficient_2", - ">i2"), - ("a/d_converter_temperature_conversion_coefficient_3", - ">i2"), - ("a/d_converter_temperature_conversion_coefficient_4", - ">i2"), - ("a/d_converter_temperature_conversion_coefficient_5", - ">i2"), - ("reserved15", ">i2"), - ("detector_#4_bias_voltage_conversion_coefficient_1", - ">i2"), - ("detector_#4_bias_voltage_conversion_coefficient_2", - ">i2"), - ("detector_#4_bias_voltage_conversion_coefficient_3", - ">i2"), - ("detector_#4_bias_voltage_conversion_coefficient_4", - ">i2"), - ("detector_#4_bias_voltage_conversion_coefficient_5", - ">i2"), - ("reserved16", ">i2"), - ("detector_#5_bias_voltage_conversion_coefficient_1", - ">i2"), - ("detector_#5_bias_voltage_conversion_coefficient_2", - ">i2"), - ("detector_#5_bias_voltage_conversion_coefficient_3", - ">i2"), - ("detector_#5_bias_voltage_conversion_coefficient_4", - ">i2"), - ("detector_#5_bias_voltage_conversion_coefficient_5", - ">i2"), - ("reserved17", ">i2"), - ("channel_3b_blackbody_view_conversion_coefficient_1", - ">i2"), - ("channel_3b_blackbody_view_conversion_coefficient_2", - ">i2"), - ("channel_3b_blackbody_view_conversion_coefficient_3", - ">i2"), - ("channel_3b_blackbody_view_conversion_coefficient_4", - ">i2"), - ("channel_3b_blackbody_view_conversion_coefficient_5", - ">i2"), - ("reserved18", ">i2"), - ("channel_4_blackbody_view_conversion_coefficient_1", - ">i2"), - ("channel_4_blackbody_view_conversion_coefficient_2", - ">i2"), - ("channel_4_blackbody_view_conversion_coefficient_3", - ">i2"), - ("channel_4_blackbody_view_conversion_coefficient_4", - ">i2"), - ("channel_4_blackbody_view_conversion_coefficient_5", - ">i2"), - ("reserved19", ">i2"), - ("channel_5_blackbody_view_conversion_coefficient_1", - ">i2"), - ("channel_5_blackbody_view_conversion_coefficient_2", - ">i2"), - ("channel_5_blackbody_view_conversion_coefficient_3", - ">i2"), - ("channel_5_blackbody_view_conversion_coefficient_4", - ">i2"), - ("channel_5_blackbody_view_conversion_coefficient_5", - ">i2"), - ("reserved20", ">i2"), - ("reference_voltage_conversion_coefficient_1", ">i2"), - ("reference_voltage_conversion_coefficient_2", ">i2"), - ("reference_voltage_conversion_coefficient_3", ">i2"), - ("reference_voltage_conversion_coefficient_4", ">i2"), - ("reference_voltage_conversion_coefficient_5", ">i2"), - ("reserved21", ">i2")]) -# FILLER -# ("zero_fill9", ">i4", (4608 - 688, ))]) +from pygac.gac_reader import GACReader +from pygac.klm_reader import KLMReader +LOG = logging.getLogger(__name__) # video data object @@ -362,7 +52,10 @@ ("zero_fill0", ">i2", (5, )), # QUALITY INDICATORS ("quality_indicator_bit_field", ">u4"), - ("scan_line_quality_flags", ">u4"), + ("scan_line_quality_flags", [("reserved", ">u1"), + ("time_problem_code", ">u1"), + ("calibration_problem_code", ">u1"), + ("earth_location_problem_code", ">u1")]), ("calibration_quality_flags", ">u2", (3, )), ("count_of_bit_errors_in_frame_sync", ">u2"), ("zero_fill1", ">i4", (2, )), @@ -430,8 +123,10 @@ ("ir_test_cal_ch_5_coefficient_1", ">i4"), ("ir_test_cal_ch_5_coefficient_2", ">i4"), ("ir_test_cal_ch_5_coefficient_3", ">i4"), - ("zero_fill2", ">i4", (3, )), # NAVIGATION + ("computed_yaw_steering", ">i2", (3,)), # only in version 5 + ("total_applied_attitude_correction", + ">i2", (3,)), # only in version 5 ("navigation_status_bit_field", ">u4"), ("time_associated_with_tip_euler_angles", ">u4"), ("tip_euler_angles", ">i2", (3, )), @@ -444,7 +139,10 @@ ("frame_sync", ">u2", (6, )), ("id", ">u2", (2, )), ("time_code", ">u2", (4, )), - ("telemetry", ">u2", (10, )), + ('telemetry', [("ramp_calibration", '>u2', (5, )), + ("PRT", '>u2', (3, )), + ("ch3_patch_temp", '>u2'), + ("spare", '>u2'), ]), ("back_scan", ">u2", (30, )), ("space_data", ">u2", (50, )), ("sync_delta", ">u2"), @@ -458,7 +156,28 @@ ("zero_fill7", ">i4", (3, )), # ANALOG HOUSEKEEPING DATA (TIP) ("invalid_word_bit_flags2", ">u4"), - ("word_1_patch_temperature", ">u1", (22, )), + ("patch_temperature_range", ">u1"), + ("patch_temperature_extended", ">u1"), + ("patch_power", ">u1"), + ("radiator_temperature", ">u1"), + ("blackbody_temperature_1", ">u1"), + ("blackbody_temperature_2", ">u1"), + ("blackbody_temperature_3", ">u1"), + ("blackbody_temperature_4", ">u1"), + ("electronics_current", ">u1"), + ("motor_current", ">u1"), + ("earth_shield_position", ">u1"), + ("electronics_temperature", ">u1"), + ("cooler_housing_temperature", ">u1"), + ("baseplate_temperature", ">u1"), + ("motor_housing_temperature", ">u1"), + ("a/d_converter_temperature", ">u1"), + ("detector_#4_bias_voltage", ">u1"), + ("detector_#5_bias_voltage", ">u1"), + ("blackbody_temperature_channel3b", ">u1"), + ("blackbody_temperature_channel4", ">u1"), + ("blackbody_temperature_channel5", ">u1"), + ("reference_voltage", ">u1"), ("zero_fill8", ">i2", (3, )), # CLOUDS FROM AVHRR (CLAVR) ("reserved0", ">u4"), @@ -468,167 +187,11 @@ ("zero_fill9", ">i4", (112, ))]) -@inherit_doc -class GACKLMReader(GACReader): - - spacecraft_names = {4: 'noaa15', - 2: 'noaa16', - 6: 'noaa17', - 7: 'noaa18', - 8: 'noaa19', - 12: 'metopa', - 11: 'metopb', - } - spacecrafts_orbital = {4: 'noaa 15', - 2: 'noaa 16', - 6: 'noaa 17', - 7: 'noaa 18', - 8: 'noaa 19', - 12: 'metop 02', - 11: 'metop 01', - } - - tsm_affected_intervals = TSM_AFFECTED_INTERVALS_KLM - - def read(self, filename): - super(GACKLMReader, self).read(filename=filename) - - with open(filename) as fd_: - self.head = np.fromfile(fd_, dtype=header, count=1)[0] - fd_.seek(4608, 0) - self.scans = np.fromfile( - fd_, dtype=scanline, count=self.head["count_of_data_records"]) - self.correct_scan_line_numbers() - self.spacecraft_id = self.head["noaa_spacecraft_identification_code"] - self.spacecraft_name = self.spacecraft_names[self.spacecraft_id] - - def get_telemetry(self): - prt_counts = np.mean(self.scans["telemetry"][:, 5:8], axis=1) - - # getting ICT counts - - ict_counts = np.zeros((len(self.scans), 3)) - ict_counts[:, 0] = np.mean(self.scans["back_scan"][:, 0::3], axis=1) - ict_counts[:, 1] = np.mean(self.scans["back_scan"][:, 1::3], axis=1) - ict_counts[:, 2] = np.mean(self.scans["back_scan"][:, 2::3], axis=1) - - # getting space counts - - space_counts = np.zeros((len(self.scans), 3)) - space_counts[:, 0] = np.mean(self.scans["space_data"][:, 2::5], axis=1) - space_counts[:, 1] = np.mean(self.scans["space_data"][:, 3::5], axis=1) - space_counts[:, 2] = np.mean(self.scans["space_data"][:, 4::5], axis=1) - - return prt_counts, ict_counts, space_counts - - def _get_lonlat(self): - lats = self.scans["earth_location"][:, 0::2] / 1e4 - lons = self.scans["earth_location"][:, 1::2] / 1e4 - return lons, lats - - def get_header_timestamp(self): - year = self.head['start_of_data_set_year'] - jday = self.head['start_of_data_set_day_of_year'] - msec = self.head['start_of_data_set_utc_time_of_day'] - try: - return self.to_datetime(self.to_datetime64(year=year, jday=jday, - msec=msec)) - except ValueError as err: - raise ValueError('Corrupt header timestamp: {0}'.format(err)) - - def _get_times(self): - year = self.scans["scan_line_year"] - jday = self.scans["scan_line_day_of_year"] - msec = self.scans["scan_line_utc_time_of_day"] - return year, jday, msec - - def get_ch3_switch(self): - """Channel 3 identification. - - 0: Channel 3b (Brightness temperature - 1: Channel 3a (Reflectance) - 2: Transition (No data) - """ - return self.scans["scan_line_bit_field"][:] & 3 - - def _get_corrupt_mask(self): - """Get mask for corrupt scanlines.""" - mask = ((self.scans["quality_indicator_bit_field"] >> 31) | - ((self.scans["quality_indicator_bit_field"] << 3) >> 31) | - ((self.scans["quality_indicator_bit_field"] << 4) >> 31)) - return mask.astype(bool) - - def get_qual_flags(self): - """Read quality flags.""" - number_of_scans = self.scans["telemetry"].shape[0] - qual_flags = np.zeros((int(number_of_scans), 7)) - qual_flags[:, 0] = self.scans["scan_line_number"] - qual_flags[:, 1] = (self.scans["quality_indicator_bit_field"] >> 31) - qual_flags[:, 2] = ( - (self.scans["quality_indicator_bit_field"] << 3) >> 31) - qual_flags[:, 3] = ( - (self.scans["quality_indicator_bit_field"] << 4) >> 31) - qual_flags[:, 4] = ( - (self.scans["quality_indicator_bit_field"] << 24) >> 30) - qual_flags[:, 5] = ( - (self.scans["quality_indicator_bit_field"] << 26) >> 30) - qual_flags[:, 6] = ( - (self.scans["quality_indicator_bit_field"] << 28) >> 30) - - return qual_flags - - def postproc(self, channels): - """Apply KLM specific postprocessing. - - Masking out 3a/3b/transition. - """ - switch = self.get_ch3_switch() - channels[:, :, 2][switch == 0] = np.nan - channels[:, :, 3][switch == 1] = np.nan - channels[:, :, 2][switch == 2] = np.nan - channels[:, :, 3][switch == 2] = np.nan - - def _adjust_clock_drift(self): - """Clock drift correction is only applied to POD satellites.""" - pass - - def get_tsm_pixels(self, channels): - """Determine pixels affected by the scan motor issue. - - Uses channels 1, 2, 4 and 5. Neither 3a, nor 3b. - """ - return get_tsm_idx(channels[:, :, 0], channels[:, :, 1], - channels[:, :, 4], channels[:, :, 5]) - - -def main(filename, start_line, end_line): - tic = datetime.datetime.now() - reader = GACKLMReader() - reader.read(filename) - reader.get_lonlat() - channels = reader.get_calibrated_channels() - sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles() - qual_flags = reader.get_qual_flags() - if np.all(reader.mask): - print("ERROR: All data is masked out. Stop processing") - raise ValueError("All data is masked out.") - - gac_io.save_gac(reader.spacecraft_name, - reader.utcs, - reader.lats, reader.lons, - channels[:, :, 0], channels[:, :, 1], - channels[:, :, 2], - channels[:, :, 3], - channels[:, :, 4], - channels[:, :, 5], - sun_zen, sat_zen, sun_azi, sat_azi, rel_azi, - qual_flags, start_line, end_line, - reader.filename, - reader.get_midnight_scanline(), - reader.get_miss_lines()) - LOG.info("pygac took: %s", str(datetime.datetime.now() - tic)) - +class GACKLMReader(GACReader, KLMReader): + """The GAC KLM reader class.""" -if __name__ == "__main__": - import sys - main(sys.argv[1], sys.argv[2], sys.argv[3]) + def __init__(self, *args, **kwargs): + """Init the GAC KLM reader.""" + GACReader.__init__(self, *args, **kwargs) + self.scanline_type = scanline + self.offset = 4608 diff --git a/pygac/gac_pod.py b/pygac/gac_pod.py index 552724c5..3c344781 100644 --- a/pygac/gac_pod.py +++ b/pygac/gac_pod.py @@ -1,13 +1,12 @@ -#!/usr/bin/env python - -# -*- coding: utf-8 -*- - -# Copyright (c) 2014, 2015 Abhay Devasthale +#!/usr/bin/python +# Copyright (c) 2014-2016. +# # Author(s): # Abhay Devasthale # Adam Dybbroe +# Sajid Pareeth # Martin Raspaud # This work was done in the framework of ESA-CCI-Clouds phase I @@ -25,127 +24,24 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . - - -"""Read a gac file. +"""Read a GAC POD file. Format specification can be found here: http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/c2/sec2-0.htm http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/c3/sec3-1.htm - """ -from __future__ import print_function - -import numpy as np -from pygac.gac_reader import GACReader, inherit_doc -from pygac.correct_tsm_issue import TSM_AFFECTED_INTERVALS_POD, get_tsm_idx -import datetime -from pygac import gac_io - import logging -LOG = logging.getLogger(__name__) - -# common header -header0 = np.dtype([("noaa_spacecraft_identification_code", ">u1"), - ("data_type_code", ">u1"), - ("start_time", ">u2", (3, )), - ("number_of_scans", ">u2"), - ("end_time", ">u2", (3, ))]) -# until 8 september 1992 -header1 = np.dtype([("noaa_spacecraft_identification_code", ">u1"), - ("data_type_code", ">u1"), - ("start_time", ">u2", (3, )), - ("number_of_scans", ">u2"), - ("end_time", ">u2", (3, )), - ("processing_block_id", "S7"), - ("ramp_auto_calibration", ">u1"), - ("number_of_data_gaps", ">u2"), - ("dacs_quality", ">u1", (6, )), - ("calibration_parameter_id", ">i2"), - ("dacs_status", ">u1"), - ("spare1", ">i1", (5, )), - ("data_set_name", "S44"), - ("spare2", ">u2", (1568, ))]) - -header2 = np.dtype([("noaa_spacecraft_identification_code", ">u1"), - ("data_type_code", ">u1"), - ("start_time", ">u2", (3, )), - ("number_of_scans", ">u2"), - ("end_time", ">u2", (3, )), - ("processing_block_id", "S7"), - ("ramp_auto_calibration", ">u1"), - ("number_of_data_gaps", ">u2"), - ("dacs_quality", ">u1", (6, )), - ("calibration_parameter_id", ">i2"), - ("dacs_status", ">u1"), - ("spare1", ">i1", (5, )), - ("data_set_name", "S42"), - ("blankfill", "S2"), - ("year_of_epoch", ">u2"), - ("julian_day_of_epoch", ">u2"), - ("millisecond_utc_epoch_time_of_day", ">u4"), - # Keplerian orbital elements - ("semi_major_axis", ">f8"), - ("eccentricity", ">f8"), - ("inclination", ">f8"), - ("argument_of_perigee", ">f8"), - ("right_ascension", ">f8"), - ("mean_anomaly", ">f8"), - # cartesian inertial true date of elements - ("x_component_of_position_vector", ">f8"), - ("y_component_of_position_vector", ">f8"), - ("z_component_of_position_vector", ">f8"), - ("x_dot_component_of_position_vector", ">f8"), - ("y_dot_component_of_position_vector", ">f8"), - ("z_dot_component_of_position_vector", ">f8"), - ("spare2", ">u2", (1516, ))]) +import numpy as np -header3 = np.dtype([("noaa_spacecraft_identification_code", ">u1"), - ("data_type_code", ">u1"), - ("start_time", ">u2", (3, )), - ("number_of_scans", ">u2"), - ("end_time", ">u2", (3, )), - ("processing_block_id", "S7"), - ("ramp_auto_calibration", ">u1"), - ("number_of_data_gaps", ">u2"), - ("dacs_quality", ">u1", (6, )), - ("calibration_parameter_id", ">i2"), - ("dacs_status", ">u1"), - ("reserved_for_mounting_and_fixed_attitude_correction_indicator", - ">i1"), - ("nadir_earth_location_tolerance", ">i1"), - ("spare1", ">i1"), - ("start_of_data_set_year", ">u2"), - ("data_set_name", "S44"), - ("year_of_epoch", ">u2"), - ("julian_day_of_epoch", ">u2"), - ("millisecond_utc_epoch_time_of_day", ">u4"), - # Keplerian orbital elements - ("semi_major_axis", ">i4"), - ("eccentricity", ">i4"), - ("inclination", ">i4"), - ("argument_of_perigee", ">i4"), - ("right_ascension", ">i4"), - ("mean_anomaly", ">i4"), - # cartesian inertial true date of elements - ("x_component_of_position_vector", ">i4"), - ("y_component_of_position_vector", ">i4"), - ("z_component_of_position_vector", ">i4"), - ("x_dot_component_of_position_vector", ">i4"), - ("y_dot_component_of_position_vector", ">i4"), - ("z_dot_component_of_position_vector", ">i4"), - # future use - ("yaw_fixed_error_correction", ">i2"), - ("roll_fixed_error_correction", ">i2"), - ("pitch_fixed_error_correction", ">i2"), - ("spare2", ">u2", (1537, ))]) +from pygac.pod_reader import PODReader +from pygac.gac_reader import GACReader +LOG = logging.getLogger(__name__) scanline = np.dtype([("scan_line_number", ">i2"), ("time_code", ">u2", (3, )), - # ("time_code", ">u1", (6, )), ("quality_indicators", ">u4"), ("calibration_coefficients", ">i4", (10, )), ("number_of_meaningful_zenith_angles_and_earth_location_appended", @@ -159,282 +55,12 @@ ("spare3", "u2", (11, ))]) -@inherit_doc -class GACPODReader(GACReader): - - spacecrafts_orbital = {25: 'tiros n', - 2: 'noaa 6', - 4: 'noaa 7', - 6: 'noaa 8', - 7: 'noaa 9', - 8: 'noaa 10', - 1: 'noaa 11', - 5: 'noaa 12', - 3: 'noaa 14', - } - spacecraft_names = {25: 'tirosn', - 2: 'noaa6', - 4: 'noaa7', - 6: 'noaa8', - 7: 'noaa9', - 8: 'noaa10', - 1: 'noaa11', - 5: 'noaa12', - 3: 'noaa14', - } - - tsm_affected_intervals = TSM_AFFECTED_INTERVALS_POD - - def correct_scan_line_numbers(self): - # Perform common corrections first. - super(GACPODReader, self).correct_scan_line_numbers() - - # cleaning up the data - min_scanline_number = np.amin(np.absolute(self.scans["scan_line_number"][:])) - if self.scans["scan_line_number"][0] == self.scans["scan_line_number"][-1] + 1: - while self.scans["scan_line_number"][0] != min_scanline_number: - self.scans = np.roll(self.scans, -1) - else: - while self.scans["scan_line_number"][0] != min_scanline_number: - self.scans = self.scans[1:] - - self.scans = self.scans[self.scans["scan_line_number"] != 0] - - def read(self, filename): - super(GACPODReader, self).read(filename=filename) - - # choose the right header depending on the date - with open(filename) as fd_: - head = np.fromfile(fd_, dtype=header0, count=1)[0] - year, jday, _ = self.decode_timestamps(head["start_time"]) - - start_date = (datetime.date(year, 1, 1) + datetime.timedelta(days=int(jday) - 1)) - - if start_date < datetime.date(1992, 9, 8): - header = header1 - elif start_date <= datetime.date(1994, 11, 15): - header = header2 - else: - header = header3 - - with open(filename) as fd_: - self.head = np.fromfile(fd_, dtype=header, count=1)[0] - self.scans = np.fromfile(fd_, - dtype=scanline, - count=self.head["number_of_scans"]) - self.correct_scan_line_numbers() - self.spacecraft_id = self.head["noaa_spacecraft_identification_code"] - if self.spacecraft_id == 1 and start_date < datetime.date(1982, 1, 1): - self.spacecraft_id = 25 - self.spacecraft_name = self.spacecraft_names[self.spacecraft_id] - LOG.info( - "Reading %s data", self.spacecrafts_orbital[self.spacecraft_id]) - - return self.head, self.scans - - def get_header_timestamp(self): - year, jday, msec = self.decode_timestamps(self.head["start_time"]) - try: - return self.to_datetime(self.to_datetime64(year=year, jday=jday, - msec=msec)) - except ValueError as err: - raise ValueError('Corrupt header timestamp: {0}'.format(err)) - - @staticmethod - def decode_timestamps(encoded): - """ - Decode timestamps. - @return: year, day of year, milliseconds since 00:00 - """ - ndims = len(encoded.shape) - if ndims == 1: - # Single header timestamp - enc0 = encoded[0] - enc1 = encoded[1] - enc2 = encoded[2] - elif ndims == 2: - # Scanline timestamps - enc0 = encoded[:, 0] - enc1 = encoded[:, 1] - enc2 = encoded[:, 2] - else: - raise ValueError('Invalid timestamp dimension') - - year = enc0 >> 9 - year = np.where(year > 75, year + 1900, year + 2000) - jday = (enc0 & 0x1FF) - msec = ((np.uint32(enc1 & 2047) << 16) | (np.uint32(enc2))) - - return year, jday, msec - - def _get_times(self): - return self.decode_timestamps(self.scans["time_code"]) - - def _adjust_clock_drift(self): - """Adjust the geolocation to compensate for the clock error. - - TODO: bad things might happen when scanlines are skipped. - """ - tic = datetime.datetime.now() - self.get_times() - from pygac.clock_offsets_converter import get_offsets - try: - offset_times, clock_error = get_offsets(self.spacecraft_name) - except KeyError: - LOG.info("No clock drift info available for %s", - self.spacecraft_name) - else: - offset_times = np.array(offset_times, dtype='datetime64[ms]') - offsets = np.interp(self.utcs.astype(np.uint64), - offset_times.astype(np.uint64), - clock_error) - LOG.info("Adjusting for clock drift of %s to %s", - str(min(offsets)), - str(max(offsets))) - self.times = (self.utcs + - offsets.astype('timedelta64[s]')).astype(datetime.datetime) - offsets *= -2 - - int_offsets = np.floor(offsets).astype(np.int) - - # filling out missing geolocations with computation from pyorbital. - line_indices = (self.scans["scan_line_number"] - + int_offsets) - - missed = sorted((set(line_indices) | - set(line_indices + 1)) - - set(self.scans["scan_line_number"])) - - min_idx = min(line_indices) - max_idx = max(max(line_indices), - max(self.scans["scan_line_number"] - min_idx)) + 1 - idx_len = max_idx - min_idx + 2 - - complete_lons = np.full((idx_len, self.lats.shape[1]), np.nan, - dtype=np.float) - complete_lats = np.full((idx_len, self.lats.shape[1]), np.nan, - dtype=np.float) - - complete_lons[self.scans["scan_line_number"] - min_idx] = self.lons - complete_lats[self.scans["scan_line_number"] - min_idx] = self.lats - - missed_utcs = ((np.array(missed) - 1) * np.timedelta64(500, "ms") - + self.utcs[0]) - - mlons, mlats = self.compute_lonlat(width=self.lats.shape[1], - utcs=missed_utcs, - clock_drift_adjust=True) - - complete_lons[missed - min_idx] = mlons - complete_lats[missed - min_idx] = mlats - - from pygac.slerp import slerp - off = offsets - np.floor(offsets) - res = slerp(complete_lons[line_indices - min_idx, :], - complete_lats[line_indices - min_idx, :], - complete_lons[line_indices - min_idx + 1, :], - complete_lats[line_indices - min_idx + 1, :], - off[:, np.newaxis, np.newaxis]) - - self.lons = res[:, :, 0] - self.lats = res[:, :, 1] - self.utcs += offsets.astype('timedelta64[s]') - - toc = datetime.datetime.now() - LOG.debug("clock drift adjustment took %s", str(toc - tic)) - - def _get_lonlat(self): - lats = self.scans["earth_location"][:, 0::2] / 128.0 - lons = self.scans["earth_location"][:, 1::2] / 128.0 - return lons, lats - - def get_telemetry(self): - number_of_scans = self.scans["telemetry"].shape[0] - decode_tele = np.zeros((int(number_of_scans), 105)) - decode_tele[:, ::3] = (self.scans["telemetry"] >> 20) & 1023 - decode_tele[:, 1::3] = (self.scans["telemetry"] >> 10) & 1023 - decode_tele[:, 2::3] = self.scans["telemetry"] & 1023 - - prt_counts = np.mean(decode_tele[:, 17:20], axis=1) - - # getting ICT counts - - ict_counts = np.zeros((int(number_of_scans), 3)) - ict_counts[:, 0] = np.mean(decode_tele[:, 22:50:3], axis=1) - ict_counts[:, 1] = np.mean(decode_tele[:, 23:51:3], axis=1) - ict_counts[:, 2] = np.mean(decode_tele[:, 24:52:3], axis=1) - - # getting space counts - - space_counts = np.zeros((int(number_of_scans), 3)) - space_counts[:, 0] = np.mean(decode_tele[:, 54:100:5], axis=1) - space_counts[:, 1] = np.mean(decode_tele[:, 55:101:5], axis=1) - space_counts[:, 2] = np.mean(decode_tele[:, 56:102:5], axis=1) - - return prt_counts, ict_counts, space_counts - - def _get_corrupt_mask(self): - """Get mask for corrupt scanlines.""" - mask = ((self.scans["quality_indicators"] >> 31) | - ((self.scans["quality_indicators"] << 4) >> 31) | - ((self.scans["quality_indicators"] << 5) >> 31)) - return mask.astype(bool) - - def get_qual_flags(self): - """Read quality flags.""" - number_of_scans = self.scans["telemetry"].shape[0] - qual_flags = np.zeros((int(number_of_scans), 7)) - qual_flags[:, 0] = self.scans["scan_line_number"] - qual_flags[:, 1] = (self.scans["quality_indicators"] >> 31) - qual_flags[:, 2] = ((self.scans["quality_indicators"] << 4) >> 31) - qual_flags[:, 3] = ((self.scans["quality_indicators"] << 5) >> 31) - qual_flags[:, 4] = ((self.scans["quality_indicators"] << 13) >> 31) - qual_flags[:, 5] = ((self.scans["quality_indicators"] << 14) >> 31) - qual_flags[:, 6] = ((self.scans["quality_indicators"] << 15) >> 31) - - return qual_flags - - def postproc(self, channels): - """No POD specific postprocessing to be done.""" - pass - - def get_tsm_pixels(self, channels): - """Determine pixels affected by the scan motor issue. - - Uses channels 1, 2, 4 and 5. Neither 3a, nor 3b. - """ - return get_tsm_idx(channels[:, :, 0], channels[:, :, 1], - channels[:, :, 3], channels[:, :, 4]) - - -def main(filename, start_line, end_line): - tic = datetime.datetime.now() - reader = GACPODReader() - reader.read(filename) - reader.get_lonlat() - channels = reader.get_calibrated_channels() - sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles() - qual_flags = reader.get_qual_flags() - if np.all(reader.mask): - print("ERROR: All data is masked out. Stop processing") - raise ValueError("All data is masked out.") - - gac_io.save_gac(reader.spacecraft_name, - reader.utcs, - reader.lats, reader.lons, - channels[:, :, 0], channels[:, :, 1], - np.full_like(channels[:, :, 0], np.nan), - channels[:, :, 2], - channels[:, :, 3], - channels[:, :, 4], - sun_zen, sat_zen, sun_azi, sat_azi, rel_azi, - qual_flags, start_line, end_line, - reader.filename, - reader.get_midnight_scanline(), - reader.get_miss_lines()) - LOG.info("pygac took: %s", str(datetime.datetime.now() - tic)) - +class GACPODReader(GACReader, PODReader): + """The GAC POD reader class.""" -if __name__ == "__main__": - import sys - main(sys.argv[1], sys.argv[2], sys.argv[3]) + def __init__(self, *args, **kwargs): + """Init the GAC POD reader.""" + GACReader.__init__(self, *args, **kwargs) + self.scanline_type = scanline + self.offset = 3220 + self.scan_points = np.arange(3.5, 2048, 5) diff --git a/pygac/gac_reader.py b/pygac/gac_reader.py index e8b75cf3..4e72b2e3 100644 --- a/pygac/gac_reader.py +++ b/pygac/gac_reader.py @@ -20,804 +20,23 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -"""Generic reader for GAC data. Can't be used as is, has to be subclassed to add -specific read functions. -""" -from abc import ABCMeta, abstractmethod, abstractproperty -import logging -import numpy as np -import os -import six -import types +"""Generic reader for GAC data. -from pygac import (CONFIG_FILE, centered_modulus, - get_absolute_azimuth_angle_diff) -try: - import ConfigParser -except ImportError: - import configparser as ConfigParser +Can't be used as is, has to be subclassed to add specific read functions. +""" -from pyorbital.orbital import Orbital -from pyorbital import astronomy -import datetime -from pygac.calibration import calibrate_solar, calibrate_thermal +from pygac.reader import Reader import pygac.pygac_geotiepoints as gtp -LOG = logging.getLogger(__name__) - - -class GACReader(six.with_metaclass(ABCMeta)): - - scan_freq = 2.0/1000.0 - """Scanning frequency (scanlines per millisecond)""" - - def __init__(self, interpolate_coords=True, adjust_clock_drift=True, - tle_dir=None, tle_name=None, tle_thresh=7): - """ - Args: - interpolate_coords: Interpolate coordinates from every eighth pixel - to all pixels. - adjust_clock_drift: Adjust the geolocation to compensate for the - clock error (POD satellites only). - tle_dir: Directory holding TLE files - tle_name: Filename pattern of TLE files. - tle_thresh: Maximum number of days between observation and nearest - TLE - """ - self.interpolate_coords = interpolate_coords - self.adjust_clock_drift = adjust_clock_drift - self.tle_dir = tle_dir - self.tle_name = tle_name - self.tle_thresh = tle_thresh - self.head = None - self.scans = None - self.spacecraft_name = None - self.spacecraft_id = None - self.utcs = None - self.lats = None - self.lons = None - self.times = None - self.tle_lines = None - self.filename = None - self._mask = None - - @abstractmethod - def read(self, filename): - """Read GAC data. - - Args: - filename (str): Specifies the GAC file to be read. - """ - self.filename = os.path.basename(filename) - LOG.info('Reading %s', self.filename) - - @abstractmethod - def get_header_timestamp(self): - """Read start timestamp from the header. - - Returns: - datetime.datetime: Start timestamp - """ - raise NotImplementedError - - def get_counts(self): - packed_data = self.scans["sensor_data"] - gac_counts = np.zeros((len(self.scans), 409 * 5)) - gac_counts[:, 0::3] = (packed_data >> 20) & 1023 - gac_counts[:, 1::3] = (packed_data >> 10) & 1023 - gac_counts[:, 2::3] = (packed_data & 1023)[:, :-1] - gac_counts = gac_counts.reshape((-1, 409, 5)) - try: - switch = self.get_ch3_switch() - except AttributeError: - return gac_counts - else: - channels = np.zeros((len(self.scans), 409, 6), - dtype=gac_counts.dtype) - channels[:, :, :2] = gac_counts[:, :, :2] - channels[:, :, -2:] = gac_counts[:, :, -2:] - channels[:, :, 2][switch == 1] = gac_counts[:, :, 2][switch == 1] - channels[:, :, 3][switch == 0] = gac_counts[:, :, 2][switch == 0] - return channels - - @abstractmethod - def _get_times(self): - """Specifies how to read scanline timestamps from GAC data. - - Returns: - int: year - int: day of year - int: milliseconds since 00:00 - """ - raise NotImplementedError - - def get_times(self): - """Read scanline timestamps and try to correct invalid values. - - Note: - Also sets self.utcs and self.times! - - Returns: - UTC timestamps - """ - if self.utcs is None: - # Read timestamps - year, jday, msec = self._get_times() - - # Correct invalid values - year, jday, msec = self.correct_times_median(year=year, jday=jday, - msec=msec) - self.utcs = self.to_datetime64(year=year, jday=jday, msec=msec) - self.correct_times_thresh() - - # Convert timestamps to datetime objects - self.times = self.to_datetime(self.utcs) - - return self.utcs - - @staticmethod - def to_datetime64(year, jday, msec): - """Convert timestamps to numpy.datetime64 - - Args: - year: Year - jday: Day of the year (1-based) - msec: Milliseconds since 00:00 - - Returns: - numpy.datetime64: Converted timestamps - """ - return (((year - 1970).astype('datetime64[Y]') - + (jday - 1).astype('timedelta64[D]')).astype('datetime64[ms]') - + msec.astype('timedelta64[ms]')) - - @staticmethod - def to_datetime(datetime64): - """Convert numpy.datetime64 to datetime.datetime - - Args: - datetime64 (numpy.datetime64): Numpy timestamp to be converted. - - Returns: - datetime.datetime: Converted timestamp - """ - return datetime64.astype(datetime.datetime) - - def lineno2msec(self, scan_line_number): - """Compute ideal scanline timestamp based on the scanline number. - - Assumes a constant scanning frequency. - - Args: - scan_line_number: Specifies the scanline number (1-based) - - Returns: - Corresponding timestamps in milliseconds since 1970-01-01 00:00, - i.e. the first scanline has timestamp 0. - """ - return (scan_line_number - 1) / self.scan_freq - - def compute_lonlat(self, width, utcs=None, clock_drift_adjust=True): - """Compute lat/lon coordinates. - - Args: - width: Number of coordinates per scanlines - utcs: Scanline timestamps - clock_drift_adjust: If True, adjust clock drift. - """ - tle1, tle2 = self.get_tle_lines() - - scan_points = np.arange(3.5, 2048, 5) - - if utcs is None: - utcs = self.get_times() - - # adjusting clock for drift - tic = datetime.datetime.now() - if clock_drift_adjust: - from pygac.clock_offsets_converter import get_offsets - try: - offset_times, clock_error = get_offsets(self.spacecraft_name) - except KeyError: - LOG.info("No clock drift info available for %s", - self.spacecraft_name) - else: - offset_times = np.array(offset_times, dtype='datetime64[ms]') - offsets = np.interp(utcs.astype(np.uint64), - offset_times.astype(np.uint64), - clock_error) - utcs -= (offsets * 1000).astype('timedelta64[ms]') - - t = utcs[0].astype(datetime.datetime) - - if "constant_yaw_attitude_error" in self.head.dtype.fields: - rpy = np.deg2rad([self.head["constant_roll_attitude_error"] / 1e3, - self.head["constant_pitch_attitude_error"] / 1e3, - self.head["constant_yaw_attitude_error"] / 1e3]) - else: - rpy = [0, 0, 0] - - LOG.info("Using rpy: %s", str(rpy)) - - from pyorbital.geoloc_instrument_definitions import avhrr_gac - from pyorbital.geoloc import compute_pixels, get_lonlatalt - sgeom = avhrr_gac(utcs.astype(datetime.datetime), scan_points, 55.385) - s_times = sgeom.times(t) - - pixels_pos = compute_pixels((tle1, tle2), sgeom, s_times, rpy) - pos_time = get_lonlatalt(pixels_pos, s_times) - - toc = datetime.datetime.now() - - LOG.warning("Computation of geolocation: %s", str(toc - tic)) - - lons, lats = pos_time[:2] - - return lons.reshape(-1, width), lats.reshape(-1, width) - - def get_calibrated_channels(self): - channels = self.get_counts() - self.get_times() - year = self.times[0].year - delta = self.times[0].date() - datetime.date(year, 1, 1) - jday = delta.days + 1 - - # Earth-Sun distance correction factor - corr = 1.0 - 0.0334 * np.cos(2.0 * np.pi * (jday - 2) / 365.25) - - # how many reflective channels are there ? - tot_ref = channels.shape[2] - 3 - - channels[:, :, 0:tot_ref] = calibrate_solar(channels[:, :, 0:tot_ref], - np.arange(tot_ref), - year, jday, - self.spacecraft_name, - corr) - prt, ict, space = self.get_telemetry() - for chan in [3, 4, 5]: - channels[:, :, chan - 6] = calibrate_thermal( - channels[:, :, chan - 6], - prt, - ict[:, chan - 3], - space[:, chan - 3], - self.scans["scan_line_number"], - chan, - self.spacecraft_name) - - # Mask out corrupt values - channels[self.mask] = np.nan - - # Apply KLM/POD specific postprocessing - self.postproc(channels) - - # Mask pixels affected by scan motor issue - if self.is_tsm_affected(): - LOG.info('Correcting for temporary scan motor issue') - self.mask_tsm_pixels(channels) - - return channels - - def get_lonlat(self): - """Compute lat/lon coordinates. - - TODO: Switch to faster interpolator? - """ - if self.lons is None and self.lats is None: - self.lons, self.lats = self._get_lonlat() - - # Interpolate from every eighth pixel to all pixels. - if self.interpolate_coords: - self.lons, self.lats = gtp.Gac_Lat_Lon_Interpolator( - self.lons, self.lats) - - # Adjust clock drift - if self.adjust_clock_drift: - self._adjust_clock_drift() - - # Mask out corrupt scanlines - self.lons[self.mask] = np.nan - self.lats[self.mask] = np.nan - - # Mask values outside the valid range - self.lats[np.fabs(self.lats) > 90.0] = np.nan - self.lons[np.fabs(self.lons) > 180.0] = np.nan - - return self.lons, self.lats - - @abstractmethod - def _get_lonlat(self): - """KLM/POD specific readout of lat/lon coordinates.""" - raise NotImplementedError - - @property - def mask(self): - """Mask for corrupt scanlines.""" - if self._mask is None: - self._mask = self._get_corrupt_mask() - return self._mask - - @abstractmethod - def _get_corrupt_mask(self): - """KLM/POD specific readout of corrupt scanline mask.""" - raise NotImplementedError - - @abstractmethod - def postproc(self, channels): - """Apply KLM/POD specific postprocessing.""" - raise NotImplementedError - - @abstractmethod - def _adjust_clock_drift(self): - """Adjust clock drift.""" - raise NotImplementedError - - @staticmethod - def tle2datetime64(times): - """Convert TLE timestamps to numpy.datetime64 - - Args: - times (float): TLE timestamps as %y%j.1234, e.g. 18001.25 - """ - # Convert %y%j.12345 to %Y%j.12345 (valid for 1950-2049) - times = np.where(times > 50000, times + 1900000, times + 2000000) - - # Convert float to datetime64 - doys = (times % 1000).astype('int') - 1 - years = (times // 1000).astype('int') - msecs = np.rint(24 * 3600 * 1000 * (times % 1)) - times64 = (years - 1970).astype('datetime64[Y]').astype('datetime64[ms]') - times64 += doys.astype('timedelta64[D]') - times64 += msecs.astype('timedelta64[ms]') - - return times64 - - def get_tle_file(self): - """Find TLE file for the current satellite.""" - tle_dir, tle_name = self.tle_dir, self.tle_name - - # If user didn't specify TLE dir/name, try config file - if tle_dir is None or tle_name is None: - conf = ConfigParser.ConfigParser() - try: - conf.read(CONFIG_FILE) - except ConfigParser.NoSectionError: - LOG.exception('Failed reading configuration file: %s', - str(CONFIG_FILE)) - raise - - options = {} - for option, value in conf.items('tle', raw=True): - options[option] = value - - tle_dir = options['tledir'] - tle_name = options['tlename'] - - values = {"satname": self.spacecraft_name, } - tle_filename = os.path.join(tle_dir, tle_name % values) - LOG.info('TLE filename = ' + str(tle_filename)) - - return tle_filename - - def read_tle_file(self, tle_filename): - """Read TLE file.""" - with open(tle_filename, 'r') as fp_: - return fp_.readlines() - - def get_tle_lines(self): - """Find closest two line elements (TLEs) for the current orbit - - Raises: - IndexError, if the closest TLE is more than :meth:`pygac.GACReader.tle_thresh` days apart - """ - if self.tle_lines is not None: - return self.tle_lines - - self.get_times() - tle_data = self.read_tle_file(self.get_tle_file()) - sdate = np.datetime64(self.times[0], '[ms]') - dates = self.tle2datetime64( - np.array([float(line[18:32]) for line in tle_data[::2]])) - - # Find index "iindex" such that dates[iindex-1] < sdate <= dates[iindex] - # Notes: - # 1. If sdate < dates[0] then iindex = 0 - # 2. If sdate > dates[-1] then iindex = len(dates), beyond the right boundary! - iindex = np.searchsorted(dates, sdate) - - if iindex in (0, len(dates)): - if iindex == len(dates): - # Reset index if beyond the right boundary (see note 2. above) - iindex -= 1 - elif abs(sdate - dates[iindex - 1]) < abs(sdate - dates[iindex]): - # Choose the closest of the two surrounding dates - iindex -= 1 - - # Make sure the TLE we found is within the threshold - delta_days = abs(sdate - dates[iindex]) / np.timedelta64(1, 'D') - if delta_days > self.tle_thresh: - raise IndexError( - "Can't find tle data for %s within +/- %d days around %s" % - (self.spacecraft_name, self.tle_thresh, sdate)) - - if delta_days > 3: - LOG.warning("Found TLE data for %s that is %f days appart", - sdate, delta_days) - else: - LOG.debug("Found TLE data for %s that is %f days appart", - sdate, delta_days) - - # Select TLE data - tle1 = tle_data[iindex * 2] - tle2 = tle_data[iindex * 2 + 1] - self.tle_lines = tle1, tle2 - return tle1, tle2 - - def get_angles(self): - """Get azimuth and zenith angles. - - Azimuth angle definition is the same as in pyorbital, but with - different units (degrees not radians for sun azimuth angles) - and different ranges. - - Returns: - sat_azi: satellite azimuth angle - degree clockwise from north in range ]-180, 180], - sat_zentih: satellite zenith angles in degrees in range [0,90], - sun_azi: sun azimuth angle - degree clockwise from north in range ]-180, 180], - sun_zentih: sun zenith angles in degrees in range [0,90], - rel_azi: absolute azimuth angle difference in degrees between sun - and sensor in range [0, 180] - - """ - self.get_times() - self.get_lonlat() - tle1, tle2 = self.get_tle_lines() - orb = Orbital(self.spacecrafts_orbital[self.spacecraft_id], - line1=tle1, line2=tle2) - - sat_azi, sat_elev = orb.get_observer_look(self.times[:, np.newaxis], - self.lons, self.lats, 0) - - sat_zenith = 90 - sat_elev - - sun_zenith = astronomy.sun_zenith_angle(self.times[:, np.newaxis], - self.lons, self.lats) - - alt, sun_azi = astronomy.get_alt_az(self.times[:, np.newaxis], - self.lons, self.lats) - del alt - sun_azi = np.rad2deg(sun_azi) - rel_azi = get_absolute_azimuth_angle_diff(sun_azi, sat_azi) - - # Scale angles range to half open interval ]-180, 180] - sat_azi = centered_modulus(sat_azi, 360.0) - sun_azi = centered_modulus(sun_azi, 360.0) - - # Mask corrupt scanlines - for arr in (sat_azi, sat_zenith, sun_azi, sun_zenith, rel_azi): - arr[self.mask] = np.nan - - return sat_azi, sat_zenith, sun_azi, sun_zenith, rel_azi - - def correct_times_median(self, year, jday, msec): - """Replace invalid timestamps with statistical estimates (using median). - - Assumes that the majority of timestamps is ok. - - Args: - year: Year - jday: Day of the year - msec: Milliseconds since 00:00 - - Returns: - Corrected year - Corrected day of the year - Corrected milliseconds - """ - # Estimate ideal timestamps based on the scanline number. Still without - # offset, e.g. the first scanline has timestamp 1970-01-01 00:00 - msec_lineno = self.lineno2msec(self.scans["scan_line_number"]) - - jday = np.where(np.logical_or(jday < 1, jday > 366), np.median(jday), jday) - if_wrong_jday = np.ediff1d(jday, to_begin=0) - jday = np.where(if_wrong_jday < 0, max(jday), jday) - - if_wrong_msec = np.where(msec < 1) - if_wrong_msec = if_wrong_msec[0] - if len(if_wrong_msec) > 0: - if if_wrong_msec[0] != 0: - msec = msec[0] + msec_lineno - else: - msec0 = np.median(msec - msec_lineno) - msec = msec0 + msec_lineno - - if_wrong_msec = np.ediff1d(msec, to_begin=0) - msec = np.where(np.logical_and(np.logical_or(if_wrong_msec < -1000, if_wrong_msec > 1000), if_wrong_jday != 1), - msec[0] + msec_lineno, msec) - - # checking if year value is out of valid range - if_wrong_year = np.where( - np.logical_or(year < 1978, year > datetime.datetime.now().year)) - if_wrong_year = if_wrong_year[0] - if len(if_wrong_year) > 0: - # if the first scanline has valid time stamp - if if_wrong_year[0] != 0: - year = year[0] - jday = jday[0] - msec = msec[0] + msec_lineno - # Otherwise use median time stamp - else: - year = np.median(year) - jday = np.median(jday) - msec0 = np.median(msec - msec_lineno) - msec = msec0 + msec_lineno - - return year, jday, msec - - def correct_scan_line_numbers(self): - """Remove scanlines with corrupted scanline numbers - - This includes: - - Scanline numbers outside the valide range - - Scanline numbers deviating more than a certain threshold from the - ideal case (1,2,3,...N) - - Example files having corrupt scanline numbers: - - NSS.GHRR.NJ.D96144.S2000.E2148.B0720102.GC - - NSS.GHRR.NJ.D96064.S0043.E0236.B0606162.WI - - NSS.GHRR.NJ.D99286.S1818.E2001.B2466869.WI - - Returns: - Intermediate and final results (for plotting purpose) - """ - along_track = np.arange(1, len(self.scans["scan_line_number"])+1) - results = {'along_track': along_track, - 'n_orig': self.scans['scan_line_number'].copy()} - - # Remove scanlines whose scanline number is outside the valid range - within_range = np.logical_and(self.scans["scan_line_number"] < 15000, - self.scans["scan_line_number"] >= 0) - self.scans = self.scans[within_range] - - # Remove scanlines deviating more than a certain threshold from the - # ideal case (1,2,3,...N). - ideal = np.arange(1, len(self.scans["scan_line_number"])+1) - - # ... Estimate possible offset (in case some scanlines are missing in - # the beginning of the scan) - offsets = self.scans["scan_line_number"] - ideal - med_offset = np.median(offsets) - - # ... Compute difference to ideal case (1,2,3,...N) + offset - diffs = np.abs(self.scans["scan_line_number"] - (ideal + med_offset)) - - # ... Remove those scanlines whose difference is larger than a certain - # threshold. For the threshold computation we only regard nonzero - # differences. - nz_diffs = diffs[diffs > 0] - if len(nz_diffs) < 50: - # Not enough differences for reliable statistics. Use fixed - # threshold. - thresh = 500 - else: - mean_nz_diffs = np.mean(nz_diffs) - std_nz_diffs = np.std(nz_diffs) - med_nz_diffs = np.median(nz_diffs) - mad_nz_diffs = np.median(np.abs(nz_diffs - med_nz_diffs)) - if mean_nz_diffs / float(med_nz_diffs) < 3: - # Relatively small variation, keep (almost) everything - thresh = mean_nz_diffs + 3*std_nz_diffs - else: - # Large variation, filter more agressively. Use median and - # median absolute deviation (MAD) as they are less sensitive to - # outliers. However, allow differences < 500 scanlines as they - # occur quite often. - thresh = max(500, med_nz_diffs + 3*mad_nz_diffs) - self.scans = self.scans[diffs <= thresh] - - LOG.debug('Removed %s scanline(s) with corrupt scanline numbers', - str(len(along_track) - len(self.scans))) - - results.update({'n_corr': self.scans['scan_line_number'], - 'within_range': within_range, - 'diffs': diffs, - 'thresh': thresh, - 'nz_diffs': nz_diffs}) - return results - - def correct_times_thresh(self, max_diff_from_t0_head=6*60*1000, - min_frac_near_t0_head=0.01, - max_diff_from_ideal_t=10*1000): - """Correct corrupted timestamps using a threshold approach. - - The threshold approach is based on the scanline number and the header - timestamp. It also works if the majority of scanlines has corrupted - timestamps. - - The header timestamp is used as a guideline to estimate the offset - between timestamps computed from the scanline number and the actual - scanline timestamps in the data. If header timestamp and scanline - timestamps do not match, no correction is applied. - - Once the offset has been estimated, one can calculate the ideal - timestamps based on the scanline number. Timestamps deviating more than - a certain threshold from the ideal timestamps are replaced by - the ideal timestamps. - - Example files having corrupt timestamps: - - NSS.GHRR.NA.D81193.S2329.E0116.B1061214.WI - - NSS.GHRR.NL.D01035.S2342.E0135.B0192627.WI - - Args: - max_diff_from_t0_head (int): Threshold for offset estimation: A - scanline timestamp matches the header timestamp t0_head if it is - within the interval - - [t0_head - max_diff_from_t0_head, - t0_head + max_diff_from_t0_head] - - around the header timestamp. - min_frac_near_t0_head (float): Specifies the minimum fraction of - scanline timestamps matching the header timestamp required for - applying the correction. - max_diff_from_ideal_t (float): Threshold for timestamp correction: - If a scanline timestamp deviates more than max_diff_from_ideal_t - from the ideal timestamp, it is regarded as corrupt and will be - replaced with the ideal timestamp. - Returns: - Intermediate and final results (for plotting purpose) - """ - results = {} - dt64_msec = ">M8[ms]" - - # Check whether scanline number increases monotonically - n = self.scans["scan_line_number"] - results.update({'t': self.utcs.copy(), 'n': n}) - if np.any(np.diff(n) < 0): - LOG.error("Cannot perform timestamp correction. Scanline number " - "does not increase monotonically.") - results['fail_reason'] = "Scanline number jumps backwards" - return results - - # Convert time to milliseconds since 1970-01-01 - t = self.utcs.astype("i8") - try: - t0_head = np.array([self.get_header_timestamp().isoformat()], - dtype="datetime64[ms]").astype("i8")[0] - except ValueError as err: - LOG.error("Cannot perform timestamp correction: %s", err) - return - - # Compute ideal timestamps based on the scanline number. Still - # without offset, i.e. scanline 0 has timestamp 1970-01-01 00:00 - tn = self.lineno2msec(self.scans["scan_line_number"]) - - # Try to determine the timestamp t0 of the first scanline. Since none - # of the actual timestamps is trustworthy, use the header timestamp - # as a guideline. However, the header timestamp may also be corrupted, - # so we only apply corrections if there is a minimum fraction of - # scanlines whose timestamps match the header timestamp. - # - # 1) Compute offsets between actual timestamps and idealized timestamps - offsets = t - tn - - # 2) If the offsets of a certain minimum fraction of scanlines are - # within a certain interval around the header timestamp, estimate - # t0 by calculating the median offset among these timestamps. If not, - # we do not have reliable information and cannot proceed. - near_t0_head = np.where( - np.fabs(offsets - t0_head) <= max_diff_from_t0_head)[0] - results.update({'offsets': offsets, - 't0_head': t0_head, - 'max_diff_from_t0_head': max_diff_from_t0_head}) - if near_t0_head.size / float(n.size) >= min_frac_near_t0_head: - t0 = np.median(offsets[near_t0_head]) - else: - LOG.error("Timestamp mismatch. Cannot perform correction.") - results['fail_reason'] = "Timestamp mismatch" - return results - - # Add estimated offset to the ideal timestamps - tn += t0 - - # Replace timestamps deviating more than a certain threshold from the - # ideal timestamp with the ideal timestamp. - corrupt_lines = np.where(np.fabs(t - tn) > max_diff_from_ideal_t) - self.utcs[corrupt_lines] = tn[corrupt_lines].astype(dt64_msec) - LOG.debug("Corrected %s timestamp(s)", str(len(corrupt_lines[0]))) - - results.update({'tn': tn, 'tcorr': self.utcs, 't0': t0}) - return results - - @abstractproperty - def tsm_affected_intervals(self): - """Specifies time intervals being affected by the scan motor problem. - - Returns: - dict: Affected time intervals. A dictionary containing a list of - (start, end) tuples for each affected platform. Both start and - end must be datetime.datetime objects. - """ - raise NotImplementedError - - def is_tsm_affected(self): - """Determine whether this orbit is affected by the scan motor problem. - - Returns: - bool: True if the orbit is affected, False otherwise. - """ - self.get_times() - ts = self.times[0] - te = self.times[-1] - try: - for interval in self.tsm_affected_intervals[self.spacecraft_id]: - if ts >= interval[0] and te <= interval[1]: - # Found a matching interval - return True - - # No matching interval, orbit is not affected - return False - except KeyError: - # Platform is not affected at all - return False - - def get_midnight_scanline(self): - """Find the scanline where the UTC date increases by one day. - - Returns: - int: The midnight scanline if it exists and is unique. - None, else. - """ - self.get_times() - d0 = np.datetime64(datetime.date(1970, 1, 1), 'D') - days = (self.utcs.astype('datetime64[D]') - d0).astype(int) - incr = np.where(np.diff(days) == 1)[0] - if len(incr) != 1: - if len(incr) > 1: - LOG.warning('Unable to determine midnight scanline: ' - 'UTC date increases more than once. ') - return None - else: - return incr[0] - - def get_miss_lines(self): - """Find missing scanlines, i.e. scanlines which were dropped for some - reason or were never recorded. - - Returns: - Indices of missing scanlines - """ - # Compare scanline number against the ideal case (1, 2, 3, ...) and - # find the missing line numbers. - ideal = set(range(1, self.scans['scan_line_number'][-1] + 1)) - missing = sorted(ideal.difference(set(self.scans['scan_line_number']))) - return np.array(missing, dtype=int) - - def mask_tsm_pixels(self, channels): - """Mask pixels affected by the scan motor issue.""" - idx = self.get_tsm_pixels(channels) - channels[idx] = np.nan - - @abstractmethod - def get_tsm_pixels(self, channels): - """Determine pixels affected by the scan motor issue. - - Channel selection is POD/KLM specific. - """ - raise NotImplementedError +class GACReader(Reader): + """Reader for GAC data.""" -def inherit_doc(cls): - """Make a class method inherit its docstring from the parent class. + # Scanning frequency (scanlines per millisecond) + scan_freq = 2.0 / 1000.0 - Copied from http://stackoverflow.com/a/8101598/5703449 . - """ - for name, func in vars(cls).items(): - if isinstance(func, types.FunctionType) and not func.__doc__: - for parent in cls.__bases__: - parfunc = getattr(parent, name, None) - if parfunc and getattr(parfunc, '__doc__', None): - func.__doc__ = parfunc.__doc__ - break - return cls + def __init__(self, *args, **kwargs): + """Init the GAC reader.""" + super(GACReader, self).__init__(*args, **kwargs) + self.scan_width = 409 + self.lonlat_interpolator = gtp.gac_lat_lon_interpolator diff --git a/pygac/klm_reader.py b/pygac/klm_reader.py new file mode 100644 index 00000000..247d1429 --- /dev/null +++ b/pygac/klm_reader.py @@ -0,0 +1,732 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +# Copyright (c) 2014-2019 + +# Author(s): + +# Abhay Devasthale +# Sajid Pareeth +# Martin Raspaud +# Adam Dybbroe + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +"""Read KLM data. + +Reads L1b GAC/LAC data from KLM series of satellites (NOAA-15 and later) and does most of the computations. +Format specification can be found here: +http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c8/sec83142-1.htm +""" + +import logging + +import numpy as np + +from pygac.correct_tsm_issue import TSM_AFFECTED_INTERVALS_KLM, get_tsm_idx +from pygac.reader import Reader + +LOG = logging.getLogger(__name__) + +# GAC header object + +header = np.dtype([("data_set_creation_site_id", "S3"), + ("ascii_blank_=_x20", "S1"), + ("noaa_level_1b_format_version_number", ">u2"), + ("noaa_level_1b_format_version_year", ">u2"), + ("noaa_level_1b_format_version_day_of_year", ">u2"), + ("reserved_for_logical_record_length", ">u2"), + ("reserved_for_block_size", ">u2"), + ("count_of_header_records", ">u2"), + ("zero_fill0", ">i2", (3, )), + ("data_set_name", "S42"), + ("processing_block_identification", "S8"), + ("noaa_spacecraft_identification_code", ">u2"), + ("instrument_id", ">u2"), + ("data_type_code", ">u2"), + ("tip_source_code", ">u2"), + ("start_of_data_set_day_count_starting_from_0_at_00h,_1_jan_1950", + ">u4"), + ("start_of_data_set_year", ">u2"), + ("start_of_data_set_day_of_year", ">u2"), + ("start_of_data_set_utc_time_of_day", ">u4"), + ("end_of_data_set_day_count_starting_from_0_at_00h,_1_jan_1950", + ">u4"), + ("end_of_data_set_year", ">u2"), + ("end_of_data_set_day_of_year", ">u2"), + ("end_of_data_set_utc_time_of_day", ">u4"), + ("year_of_last_cpids_update", ">u2"), + ("day_of_year_of_last_cpids_update", ">u2"), + ("zero_fill1", ">i2", (4, )), + # DATA SET QUALITY INDICATORS + ("instrument_status", ">u4"), + ("zero_fill2", ">i2"), + ("record_number_of_status_change", ">u2"), + ("second_instrument_status", ">u4"), + ("count_of_data_records", ">u2"), + ("count_of_calibrated,_earth_located_scan_lines", ">u2"), + ("count_of_missing_scan_lines", ">u2"), + ("count_of_data_gaps", ">u2"), + # only in version 5 + ("count_of_data_frames_without_frame_sync_word", ">u2"), + ("count_of_pacs_detected_tip_parity_errors", ">u2"), + # only in version 5 + ("sum_of_all_auxiliary_sync_errors_detected", ">u2"), + ("time_sequence_error", ">u2"), + ("time_sequence_error_code", ">u2"), + ("socc_clock_update_indicator", ">u2"), + ("earth_location_error_indicator", ">u2"), + ("earth_location_error_code", ">u2"), + ("pacs_status_bit_field", ">u2"), + ("data_source", ">u2"), + ("zero_fill3", ">i4"), + ("reserved_for_the_ingester", "S8"), + ("reserved_for_decommutation", "S8"), + ("zero_fill4", ">i2", (5, )), + # CALIBRATION + ("ramp_auto_calibration_indicators_bit_field", ">u2"), + ("year_of_most_recent_solar_channel_calibration", ">u2"), + ("day_of_year_of_most_recent_solar_channel_calibration", + ">u2"), + ("primary_calibration_algorithm_id", ">u2"), + ("primary_calibration_algorithm_selected_options", ">u2"), + ("secondary_calibration_algorithm_id", ">u2"), + ("secondary_calibration_algorithm_selected_options", ">u2"), + ("ir_target_temperature_1_conversion_coefficient_1", ">i2"), + ("ir_target_temperature_1_conversion_coefficient_2", ">i2"), + ("ir_target_temperature_1_conversion_coefficient_3", ">i2"), + ("ir_target_temperature_1_conversion_coefficient_4", ">i2"), + ("ir_target_temperature_1_conversion_coefficient_5", ">i2"), + ("ir_target_temperature_1_conversion_coefficient_6", ">i2"), + ("ir_target_temperature_2_conversion_coefficient_1", ">i2"), + ("ir_target_temperature_2_conversion_coefficient_2", ">i2"), + ("ir_target_temperature_2_conversion_coefficient_3", ">i2"), + ("ir_target_temperature_2_conversion_coefficient_4", ">i2"), + ("ir_target_temperature_2_conversion_coefficient_5", ">i2"), + ("ir_target_temperature_2_conversion_coefficient_6", ">i2"), + ("ir_target_temperature_3_conversion_coefficient_1", ">i2"), + ("ir_target_temperature_3_conversion_coefficient_2", ">i2"), + ("ir_target_temperature_3_conversion_coefficient_3", ">i2"), + ("ir_target_temperature_3_conversion_coefficient_4", ">i2"), + ("ir_target_temperature_3_conversion_coefficient_5", ">i2"), + ("ir_target_temperature_3_conversion_coefficient_6", ">i2"), + ("ir_target_temperature_4_conversion_coefficient_1", ">i2"), + ("ir_target_temperature_4_conversion_coefficient_2", ">i2"), + ("ir_target_temperature_4_conversion_coefficient_3", ">i2"), + ("ir_target_temperature_4_conversion_coefficient_4", ">i2"), + ("ir_target_temperature_4_conversion_coefficient_5", ">i2"), + ("ir_target_temperature_4_conversion_coefficient_6", ">i2"), + ("zero_fill5", ">i4", (2, )), + # RADIANCE CONVERSION + ("ch_1_solar_filtered_irradiance_in_wavelength", ">i4"), + ("ch_1_equivalent_filter_width_in_wavelength", ">i4"), + ("ch_2_solar_filtered_irradiance_in_wavelength", ">i4"), + ("ch_2_equivalent_filter_width_in_wavelength", ">i4"), + ("ch_3a_solar_filtered_irradiance_in_wavelength", ">i4"), + ("ch_3a_equivalent_filter_width_in_wavelength", ">i4"), + ("ch_3b_central_wavenumber", ">i4"), + ("ch_3b_constant_1", ">i4"), + ("ch_3b_constant_2", ">i4"), + ("ch_4_central_wavenumber", ">i4"), + ("ch_4_constant_1", ">i4"), + ("ch_4_constant_2", ">i4"), + ("ch_5_central_wavenumber", ">i4"), + ("ch_5_constant_1", ">i4"), + ("ch_5_constant_2", ">i4"), + ("zero_fill6", ">i4", (3, )), + # NAVIGATION + ("reference_ellipsoid_model_id", "S8"), + ("nadir_earth_location_tolerance", ">u2"), + ("earth_location_bit_field", ">u2"), + ("zero_fill7", ">i2"), + ("constant_roll_attitude_error", ">i2"), + ("constant_pitch_attitude_error", ">i2"), + ("constant_yaw_attitude_error", ">i2"), + ("epoch_year_for_orbit_vector", ">u2"), + ("day_of_epoch_year_for_orbit_vector", ">u2"), + ("epoch_utc_time_of_day_for_orbit_vector", ">u4"), + ("semi-major_axis", ">i4"), + ("eccentricity", ">i4"), + ("inclination", ">i4"), + ("argument_of_perigee", ">i4"), + ("right_ascension_of_the_ascending_node", ">i4"), + ("mean_anomaly", ">i4"), + ("position_vector_x_component", ">i4"), + ("position_vector_y_component", ">i4"), + ("position_vector_z_component", ">i4"), + ("velocity_vector_x-dot_component", ">i4"), + ("velocity_vector_y-dot_component", ">i4"), + ("velocity_vector_z-dot_component", ">i4"), + ("earth_sun_distance_ratio", ">u4"), + ("zero_fill8", ">i4", (4, ))]) + +# ANALOG TELEMETRY CONVERSION +analog_telemetry_v2 = np.dtype([("patch_temperature_conversion_coefficient_1", ">i2"), + ("patch_temperature_conversion_coefficient_2", ">i2"), + ("patch_temperature_conversion_coefficient_3", ">i2"), + ("patch_temperature_conversion_coefficient_4", ">i2"), + ("patch_temperature_conversion_coefficient_5", ">i2"), + ("reserved0", ">i2"), + ("patch_temperature_extended_conversion_coefficient_1", + ">i2"), + ("patch_temperature_extended_conversion_coefficient_2", + ">i2"), + ("patch_temperature_extended_conversion_coefficient_3", + ">i2"), + ("patch_temperature_extended_conversion_coefficient_4", + ">i2"), + ("patch_temperature_extended_conversion_coefficient_5", + ">i2"), + ("reserved1", ">i2"), + ("patch_power_conversion_coefficient_1", ">i2"), + ("patch_power_conversion_coefficient_2", ">i2"), + ("patch_power_conversion_coefficient_3", ">i2"), + ("patch_power_conversion_coefficient_4", ">i2"), + ("patch_power_conversion_coefficient_5", ">i2"), + ("reserved2", ">i2"), + ("radiator_temperature_conversion_coefficient_1", ">i2"), + ("radiator_temperature_conversion_coefficient_2", ">i2"), + ("radiator_temperature_conversion_coefficient_3", ">i2"), + ("radiator_temperature_conversion_coefficient_4", ">i2"), + ("radiato r_temperature_conversion_coefficient_5", ">i2"), + ("reserved3", ">i2"), + ("blackbody_temperature_1_conversion_coefficient_1", ">i2"), + ("blackbody_temperature_1_conversion_coefficient_2", ">i2"), + ("blackbody_temperature_1_conversion_coefficient_3", ">i2"), + ("blackbody_temperature_1_conversion_coefficient_4", ">i2"), + ("blackbody_temperature_1_conversion_coefficient_5", ">i2"), + ("reserved4", ">i2"), + ("blackbody_temperature_2_conversion_coefficient_1", ">i2"), + ("blackbody_temperature_2_conversion_coefficient_2", ">i2"), + ("blackbody_temperature_2_conversion_coefficient_3", ">i2"), + ("blackbody_temperature_2_conversion_coefficient_4", ">i2"), + ("blackbody_temperature_2_conversion_coefficient_5", ">i2"), + ("reserved5", ">i2"), + ("blackbody_temperature_3_conversion_coefficient_1", ">i2"), + ("blackbody_temperature_3_conversion_coefficient_2", ">i2"), + ("blackbody_temperature_3_conversion_coefficient_3", ">i2"), + ("blackbody_temperature_3_conversion_coefficient_4", ">i2"), + ("blackbody_temperature_3_conversion_coefficient_5", ">i2"), + ("reserved6", ">i2"), + ("blackbody_temperature_4_conversion_coefficient_1", ">i2"), + ("blackbody_temperature_4_conversion_coefficient_2", ">i2"), + ("blackbody_temperature_4_conversion_coefficient_3", ">i2"), + ("blackbody_temperature_4_conversion_coefficient_4", ">i2"), + ("blackbody_temperature_4_conversion_coefficient_5", ">i2"), + ("reserved7", ">i2"), + ("electronics_current_conversion_coefficient_1", ">i2"), + ("electronics_current_conversion_coefficient_2", ">i2"), + ("electronics_current_conversion_coefficient_3", ">i2"), + ("electronics_current_conversion_coefficient_4", ">i2"), + ("electronics_current_conversion_coefficient_5", ">i2"), + ("reserved8", ">i2"), + ("motor_current_conversion_coefficient_1", ">i2"), + ("motor_current_conversion_coefficient_2", ">i2"), + ("motor_current_conversion_coefficient_3", ">i2"), + ("motor_current_conversion_coefficient_4", ">i2"), + ("motor_current_conversion_coefficient_5", ">i2"), + ("reserved9", ">i2"), + ("earth_shield_position_conversion_coefficient_1", ">i2"), + ("earth_shield_position_conversion_coefficient_2", ">i2"), + ("earth_shield_position_conversion_coefficient_3", ">i2"), + ("earth_shield_position_conversion_coefficient_4", ">i2"), + ("earth_shield_position_conversion_coefficient_5", ">i2"), + ("reserved10", ">i2"), + ("electronics_temperature_conversion_coefficient_1", ">i2"), + ("electronics_temperature_conversion_coefficient_2", ">i2"), + ("electronics_temperature_conversion_coefficient_3", ">i2"), + ("electronics_temperature_conversion_coefficient_4", ">i2"), + ("electronics_temperature_conversion_coefficient_5", ">i2"), + ("reserved11", ">i2"), + ("cooler_housing_temperature_conversion_coefficient_1", + ">i2"), + ("cooler_housing_temperature_conversion_coefficient_2", + ">i2"), + ("cooler_housing_temperature_conversion_coefficient_3", + ">i2"), + ("cooler_housing_temperature_conversion_coefficient_4", + ">i2"), + ("cooler_housing_temperature_conversion_coefficient_5", + ">i2"), + ("reserved12", ">i2"), + ("baseplate_temperature_conversion_coefficient_1", ">i2"), + ("baseplate_temperature_conversion_coefficient_2", ">i2"), + ("baseplate_temperature_conversion_coefficient_3", ">i2"), + ("baseplate_temperature_conversion_coefficient_4", ">i2"), + ("baseplate_temperature_conversion_coefficient_5", ">i2"), + ("reserved13", ">i2"), + ("motor_housing_temperature_conversion_coefficient_1", + ">i2"), + ("motor_housing_temperature_conversion_coefficient_2", + ">i2"), + ("motor_housing_temperature_conversion_coefficient_3", + ">i2"), + ("motor_housing_temperature_conversion_coefficient_4", + ">i2"), + ("motor_housing_temperature_conversion_coefficient_5", + ">i2"), + ("reserved14", ">i2"), + ("a/d_converter_temperature_conversion_coefficient_1", + ">i2"), + ("a/d_converter_temperature_conversion_coefficient_2", + ">i2"), + ("a/d_converter_temperature_conversion_coefficient_3", + ">i2"), + ("a/d_converter_temperature_conversion_coefficient_4", + ">i2"), + ("a/d_converter_temperature_conversion_coefficient_5", + ">i2"), + ("reserved15", ">i2"), + ("detector_#4_bias_voltage_conversion_coefficient_1", + ">i2"), + ("detector_#4_bias_voltage_conversion_coefficient_2", + ">i2"), + ("detector_#4_bias_voltage_conversion_coefficient_3", + ">i2"), + ("detector_#4_bias_voltage_conversion_coefficient_4", + ">i2"), + ("detector_#4_bias_voltage_conversion_coefficient_5", + ">i2"), + ("reserved16", ">i2"), + ("detector_#5_bias_voltage_conversion_coefficient_1", + ">i2"), + ("detector_#5_bias_voltage_conversion_coefficient_2", + ">i2"), + ("detector_#5_bias_voltage_conversion_coefficient_3", + ">i2"), + ("detector_#5_bias_voltage_conversion_coefficient_4", + ">i2"), + ("detector_#5_bias_voltage_conversion_coefficient_5", + ">i2"), + ("reserved17", ">i2"), + ("channel_3b_blackbody_view_conversion_coefficient_1", + ">i2"), + ("channel_3b_blackbody_view_conversion_coefficient_2", + ">i2"), + ("channel_3b_blackbody_view_conversion_coefficient_3", + ">i2"), + ("channel_3b_blackbody_view_conversion_coefficient_4", + ">i2"), + ("channel_3b_blackbody_view_conversion_coefficient_5", + ">i2"), + ("reserved18", ">i2"), + ("channel_4_blackbody_view_conversion_coefficient_1", + ">i2"), + ("channel_4_blackbody_view_conversion_coefficient_2", + ">i2"), + ("channel_4_blackbody_view_conversion_coefficient_3", + ">i2"), + ("channel_4_blackbody_view_conversion_coefficient_4", + ">i2"), + ("channel_4_blackbody_view_conversion_coefficient_5", + ">i2"), + ("reserved19", ">i2"), + ("channel_5_blackbody_view_conversion_coefficient_1", + ">i2"), + ("channel_5_blackbody_view_conversion_coefficient_2", + ">i2"), + ("channel_5_blackbody_view_conversion_coefficient_3", + ">i2"), + ("channel_5_blackbody_view_conversion_coefficient_4", + ">i2"), + ("channel_5_blackbody_view_conversion_coefficient_5", + ">i2"), + ("reserved20", ">i2"), + ("reference_voltage_conversion_coefficient_1", ">i2"), + ("reference_voltage_conversion_coefficient_2", ">i2"), + ("reference_voltage_conversion_coefficient_3", ">i2"), + ("reference_voltage_conversion_coefficient_4", ">i2"), + ("reference_voltage_conversion_coefficient_5", ">i2"), + ("reserved21", ">i2")]) +# FILLER +# ("zero_fill9", ">i4", (4608 - 688, ))]) + + +analog_telemetry_v5 = np.dtype([("patch_temperature_conversion_coefficient_1", ">i4"), + ("patch_temperature_conversion_coefficient_2", ">i4"), + ("patch_temperature_conversion_coefficient_3", ">i4"), + ("patch_temperature_conversion_coefficient_4", ">i4"), + ("patch_temperature_conversion_coefficient_5", ">i4"), + ("patch_temperature_conversion_coefficient_6", ">i4"), + ("patch_temperature_extended_conversion_coefficient_1", + ">i4"), + ("patch_temperature_extended_conversion_coefficient_2", + ">i4"), + ("patch_temperature_extended_conversion_coefficient_3", + ">i4"), + ("patch_temperature_extended_conversion_coefficient_4", + ">i4"), + ("patch_temperature_extended_conversion_coefficient_5", + ">i4"), + ("patch_temperature_extended_conversion_coefficient_6", ">i4"), + ("patch_power_conversion_coefficient_1", ">i4"), + ("patch_power_conversion_coefficient_2", ">i4"), + ("patch_power_conversion_coefficient_3", ">i4"), + ("patch_power_conversion_coefficient_4", ">i4"), + ("patch_power_conversion_coefficient_5", ">i4"), + ("patch_power_conversion_coefficient_6", ">i4"), + ("radiator_temperature_conversion_coefficient_1", ">i4"), + ("radiator_temperature_conversion_coefficient_2", ">i4"), + ("radiator_temperature_conversion_coefficient_3", ">i4"), + ("radiator_temperature_conversion_coefficient_4", ">i4"), + ("radiator_temperature_conversion_coefficient_5", ">i4"), + ("radiator_temperature_conversion_coefficient_6", ">i4"), + ("blackbody_temperature_1_conversion_coefficient_1", ">i4"), + ("blackbody_temperature_1_conversion_coefficient_2", ">i4"), + ("blackbody_temperature_1_conversion_coefficient_3", ">i4"), + ("blackbody_temperature_1_conversion_coefficient_4", ">i4"), + ("blackbody_temperature_1_conversion_coefficient_5", ">i4"), + ("blackbody_temperature_1_conversion_coefficient_6", ">i4"), + ("blackbody_temperature_2_conversion_coefficient_1", ">i4"), + ("blackbody_temperature_2_conversion_coefficient_2", ">i4"), + ("blackbody_temperature_2_conversion_coefficient_3", ">i4"), + ("blackbody_temperature_2_conversion_coefficient_4", ">i4"), + ("blackbody_temperature_2_conversion_coefficient_5", ">i4"), + ("blackbody_temperature_2_conversion_coefficient_6", ">i4"), + ("blackbody_temperature_3_conversion_coefficient_1", ">i4"), + ("blackbody_temperature_3_conversion_coefficient_2", ">i4"), + ("blackbody_temperature_3_conversion_coefficient_3", ">i4"), + ("blackbody_temperature_3_conversion_coefficient_4", ">i4"), + ("blackbody_temperature_3_conversion_coefficient_5", ">i4"), + ("blackbody_temperature_3_conversion_coefficient_6", ">i4"), + ("blackbody_temperature_4_conversion_coefficient_1", ">i4"), + ("blackbody_temperature_4_conversion_coefficient_2", ">i4"), + ("blackbody_temperature_4_conversion_coefficient_3", ">i4"), + ("blackbody_temperature_4_conversion_coefficient_4", ">i4"), + ("blackbody_temperature_4_conversion_coefficient_5", ">i4"), + ("blackbody_temperature_4_conversion_coefficient_6", ">i4"), + ("electronics_current_conversion_coefficient_1", ">i4"), + ("electronics_current_conversion_coefficient_2", ">i4"), + ("electronics_current_conversion_coefficient_3", ">i4"), + ("electronics_current_conversion_coefficient_4", ">i4"), + ("electronics_current_conversion_coefficient_5", ">i4"), + ("electronics_current_conversion_coefficient_6", ">i4"), + ("motor_current_conversion_coefficient_1", ">i4"), + ("motor_current_conversion_coefficient_2", ">i4"), + ("motor_current_conversion_coefficient_3", ">i4"), + ("motor_current_conversion_coefficient_4", ">i4"), + ("motor_current_conversion_coefficient_5", ">i4"), + ("motor_current_conversion_coefficient_6", ">i4"), + ("earth_shield_position_conversion_coefficient_1", ">i4"), + ("earth_shield_position_conversion_coefficient_2", ">i4"), + ("earth_shield_position_conversion_coefficient_3", ">i4"), + ("earth_shield_position_conversion_coefficient_4", ">i4"), + ("earth_shield_position_conversion_coefficient_5", ">i4"), + ("earth_shield_position_conversion_coefficient_6", ">i4"), + ("electronics_temperature_conversion_coefficient_1", ">i4"), + ("electronics_temperature_conversion_coefficient_2", ">i4"), + ("electronics_temperature_conversion_coefficient_3", ">i4"), + ("electronics_temperature_conversion_coefficient_4", ">i4"), + ("electronics_temperature_conversion_coefficient_5", ">i4"), + ("electronics_temperature_conversion_coefficient_6", ">i4"), + ("cooler_housing_temperature_conversion_coefficient_1", ">i4"), + ("cooler_housing_temperature_conversion_coefficient_2", ">i4"), + ("cooler_housing_temperature_conversion_coefficient_3", ">i4"), + ("cooler_housing_temperature_conversion_coefficient_4", ">i4"), + ("cooler_housing_temperature_conversion_coefficient_5", ">i4"), + ("cooler_housing_temperature_conversion_coefficient_6", ">i4"), + ("baseplate_temperature_conversion_coefficient_1", ">i4"), + ("baseplate_temperature_conversion_coefficient_2", ">i4"), + ("baseplate_temperature_conversion_coefficient_3", ">i4"), + ("baseplate_temperature_conversion_coefficient_4", ">i4"), + ("baseplate_temperature_conversion_coefficient_5", ">i4"), + ("baseplate_temperature_conversion_coefficient_6", ">i4"), + ("motor_housing_temperature_conversion_coefficient_1", ">i4"), + ("motor_housing_temperature_conversion_coefficient_2", ">i4"), + ("motor_housing_temperature_conversion_coefficient_3", ">i4"), + ("motor_housing_temperature_conversion_coefficient_4", ">i4"), + ("motor_housing_temperature_conversion_coefficient_5", ">i4"), + ("motor_housing_temperature_conversion_coefficient_6", ">i4"), + ("a/d_converter_temperature_conversion_coefficient_1", ">i4"), + ("a/d_converter_temperature_conversion_coefficient_2", ">i4"), + ("a/d_converter_temperature_conversion_coefficient_3", ">i4"), + ("a/d_converter_temperature_conversion_coefficient_4", ">i4"), + ("a/d_converter_temperature_conversion_coefficient_5", ">i4"), + ("a/d_converter_temperature_conversion_coefficient_6", ">i4"), + ("detector_#4_bias_voltage_conversion_coefficient_1", ">i4"), + ("detector_#4_bias_voltage_conversion_coefficient_2", ">i4"), + ("detector_#4_bias_voltage_conversion_coefficient_3", ">i4"), + ("detector_#4_bias_voltage_conversion_coefficient_4", ">i4"), + ("detector_#4_bias_voltage_conversion_coefficient_5", ">i4"), + ("detector_#4_bias_voltage_conversion_coefficient_6", ">i4"), + ("detector_#5_bias_voltage_conversion_coefficient_1", ">i4"), + ("detector_#5_bias_voltage_conversion_coefficient_2", ">i4"), + ("detector_#5_bias_voltage_conversion_coefficient_3", ">i4"), + ("detector_#5_bias_voltage_conversion_coefficient_4", ">i4"), + ("detector_#5_bias_voltage_conversion_coefficient_5", ">i4"), + ("detector_#5_bias_voltage_conversion_coefficient_6", ">i4"), + ("blackbody_temperature_channel3b_conversion_coefficient_1", ">i4"), + ("blackbody_temperature_channel3b_conversion_coefficient_2", ">i4"), + ("blackbody_temperature_channel3b_conversion_coefficient_3", ">i4"), + ("blackbody_temperature_channel3b_conversion_coefficient_4", ">i4"), + ("blackbody_temperature_channel3b_conversion_coefficient_5", ">i4"), + ("blackbody_temperature_channel3b_conversion_coefficient_6", ">i4"), + ("blackbody_temperature_channel4_conversion_coefficient_1", + ">i4"), + ("blackbody_temperature_channel4_conversion_coefficient_2", + ">i4"), + ("blackbody_temperature_channel4_conversion_coefficient_3", + ">i4"), + ("blackbody_temperature_channel4_conversion_coefficient_4", + ">i4"), + ("blackbody_temperature_channel4_conversion_coefficient_5", + ">i4"), + ("blackbody_temperature_channel4_conversion_coefficient_6", + ">i4"), + ("blackbody_temperature_channel5_conversion_coefficient_1", + ">i4"), + ("blackbody_temperature_channel5_conversion_coefficient_2", + ">i4"), + ("blackbody_temperature_channel5_conversion_coefficient_3", + ">i4"), + ("blackbody_temperature_channel5_conversion_coefficient_4", + ">i4"), + ("blackbody_temperature_channel5_conversion_coefficient_5", + ">i4"), + ("blackbody_temperature_channel5_conversion_coefficient_6", + ">i4"), + ("reference_voltage_conversion_coefficient_1", ">i4"), + ("reference_voltage_conversion_coefficient_2", ">i4"), + ("reference_voltage_conversion_coefficient_3", ">i4"), + ("reference_voltage_conversion_coefficient_4", ">i4"), + ("reference_voltage_conversion_coefficient_5", ">i4"), + ("reference_voltage_conversion_coefficient_6", ">i4"), + # METOP MANEUVERS IDENTIFICATION + ("start_of_maneuver_year", ">u2"), + ("start_of_maneuver_day", ">u2"), + ("start_of_maneuver_utc_time_of_day", ">u4"), + ("end_of_manuever_year", ">u2"), + ("end_of_maneuver_day", ">u2"), + ("end_of_maneuver_utc_time_of_day", ">u4"), + ("change_in_spacecraft_velocity", ">i4"), + ("change in spacecraft_mass", ">u4"), + # CLOUDS FROM AVHRR(CLAVR) + ("clavr_status_bit_field", ">u2"), + ("zero_fill9", ">i2")]) + + +ars_header = np.dtype([('COST_number', 'S6'), + ('SAA_number', 'S8'), + ('order_creation_year', 'S4'), + ('order_creation_day_of_year', 'S3'), + ('processing_site_code', 'S1'), + ('processing_software', 'S8'), + # data selection criteria + ('data_set_name', 'S42'), + ('ascii_blank_', 'S2'), + ('select_flag', 'S1'), + ('beginning_latitude', 'S3'), + ('ending_latitude', 'S3'), + ('beginning_longitude', 'S4'), + ('ending_longitude', 'S4'), + ('start_hour', 'S2'), + ('start_minute', 'S2'), + ('number_of_minutes', 'S3'), + ('appended_data_flag', 'S1'), + ('channel_select_flag', 'S1', (20, )), + # dataset summary + ('ascii_blank__', 'S29'), + ('ascend_descend_flag', 'S1'), + ('first_latitude', 'S3'), + ('last_latitude', 'S3'), + ('first_longitude', 'S4'), + ('last_longitude', 'S4'), + ('data_format', 'S20'), + ('size_of_record', 'S6'), + ('number_of_records', 'S6'), + # filler + ('ascii_blank', 'S319') + ]) + + +class KLMReader(Reader): + """Reader for KLM data.""" + + spacecraft_names = {4: 'noaa15', + 2: 'noaa16', + 6: 'noaa17', + 7: 'noaa18', + 8: 'noaa19', + 12: 'metopa', + 11: 'metopb', + } + spacecrafts_orbital = {4: 'noaa 15', + 2: 'noaa 16', + 6: 'noaa 17', + 7: 'noaa 18', + 8: 'noaa 19', + 12: 'metop 02', + 11: 'metop 01', + } + + tsm_affected_intervals = TSM_AFFECTED_INTERVALS_KLM + + def read(self, filename): + """Read the data. + + Args: + filename: string + The filename to read from. + + Returns: + header: numpy record array + The header metadata + scans: numpy record array + The scanlines + + """ + super(KLMReader, self).read(filename=filename) + with open(filename) as fd_: + self.ars_head = np.fromfile(fd_, dtype=ars_header, count=1)[0] + if not self.ars_head['data_set_name'].startswith(self.creation_site): + fd_.seek(0) + self.ars_head = None + ars_offset = 0 + else: + ars_offset = ars_header.itemsize + self.head = np.fromfile(fd_, dtype=header, count=1)[0] + self.header_version = self.head[ + "noaa_level_1b_format_version_number"] + if self.header_version >= 5: + self.analog_telemetry = np.fromfile( + fd_, dtype=analog_telemetry_v5, count=1)[0] + else: + self.analog_telemetry = np.fromfile( + fd_, dtype=analog_telemetry_v2, count=1)[0] + # LAC: 1, GAC: 2, ... + self.data_type = self.head['data_type_code'] + fd_.seek(self.offset + ars_offset, 0) + self.scans = np.fromfile(fd_, dtype=self.scanline_type, + count=self.head["count_of_data_records"]) + + self.correct_scan_line_numbers() + self.spacecraft_id = self.head["noaa_spacecraft_identification_code"] + self.spacecraft_name = self.spacecraft_names[self.spacecraft_id] + + def get_telemetry(self): + """Get the telemetry. + + Returns: + prt_counts: np.array + ict_counts: np.array + space_counts: np.array + + """ + prt_counts = np.mean(self.scans["telemetry"]['PRT'], axis=1) + + # getting ICT counts + + ict_counts = np.zeros((len(self.scans), 3)) + ict_counts[:, 0] = np.mean(self.scans["back_scan"][:, 0::3], axis=1) + ict_counts[:, 1] = np.mean(self.scans["back_scan"][:, 1::3], axis=1) + ict_counts[:, 2] = np.mean(self.scans["back_scan"][:, 2::3], axis=1) + + # getting space counts + + space_counts = np.zeros((len(self.scans), 3)) + space_counts[:, 0] = np.mean(self.scans["space_data"][:, 2::5], axis=1) + space_counts[:, 1] = np.mean(self.scans["space_data"][:, 3::5], axis=1) + space_counts[:, 2] = np.mean(self.scans["space_data"][:, 4::5], axis=1) + + return prt_counts, ict_counts, space_counts + + def _get_lonlat(self): + """Get the longitudes and latitudes.""" + lats = self.scans["earth_location"][:, 0::2] / 1e4 + lons = self.scans["earth_location"][:, 1::2] / 1e4 + return lons, lats + + def get_header_timestamp(self): + """Get the timestamp from the header. + + Returns: + A datetime object containing the timestamp from the header. + + Raises: + A ValueError if the timestamp is corrupt. + + """ + year = self.head['start_of_data_set_year'] + jday = self.head['start_of_data_set_day_of_year'] + msec = self.head['start_of_data_set_utc_time_of_day'] + try: + return self.to_datetime(self.to_datetime64(year=year, jday=jday, + msec=msec)) + except ValueError as err: + raise ValueError('Corrupt header timestamp: {0}'.format(err)) + + def _get_times(self): + """Get the times of the scanlines.""" + year = self.scans["scan_line_year"] + jday = self.scans["scan_line_day_of_year"] + msec = self.scans["scan_line_utc_time_of_day"] + return year, jday, msec + + def get_ch3_switch(self): + """Channel 3 identification. + + 0: Channel 3b (Brightness temperature + 1: Channel 3a (Reflectance) + 2: Transition (No data) + """ + return self.scans["scan_line_bit_field"][:] & 3 + + def _get_corrupt_mask(self): + """Get mask for corrupt scanlines.""" + mask = ((self.scans["quality_indicator_bit_field"] >> 31) | + ((self.scans["quality_indicator_bit_field"] << 3) >> 31) | + ((self.scans["quality_indicator_bit_field"] << 4) >> 31)) + return mask.astype(bool) + + def get_qual_flags(self): + """Read quality flags.""" + number_of_scans = self.scans["telemetry"].shape[0] + qual_flags = np.zeros((int(number_of_scans), 7)) + qual_flags[:, 0] = self.scans["scan_line_number"] + qual_flags[:, 1] = (self.scans["quality_indicator_bit_field"] >> 31) + qual_flags[:, 2] = ( + (self.scans["quality_indicator_bit_field"] << 3) >> 31) + qual_flags[:, 3] = ( + (self.scans["quality_indicator_bit_field"] << 4) >> 31) + qual_flags[:, 4] = ( + (self.scans["quality_indicator_bit_field"] << 24) >> 30) + qual_flags[:, 5] = ( + (self.scans["quality_indicator_bit_field"] << 26) >> 30) + qual_flags[:, 6] = ( + (self.scans["quality_indicator_bit_field"] << 28) >> 30) + + return qual_flags + + def postproc(self, channels): + """Apply KLM specific postprocessing. + + Masking out 3a/3b/transition. + """ + switch = self.get_ch3_switch() + channels[:, :, 2][switch == 0] = np.nan + channels[:, :, 3][switch == 1] = np.nan + channels[:, :, 2][switch == 2] = np.nan + channels[:, :, 3][switch == 2] = np.nan + + def _adjust_clock_drift(self): + """Clock drift correction is only applied to POD satellites.""" + pass + + def get_tsm_pixels(self, channels): + """Determine pixels affected by the scan motor issue. + + Uses channels 1, 2, 4 and 5. Neither 3a, nor 3b. + """ + return get_tsm_idx(channels[:, :, 0], channels[:, :, 1], + channels[:, :, 4], channels[:, :, 5]) diff --git a/pygac/lac_klm.py b/pygac/lac_klm.py new file mode 100644 index 00000000..56727a0e --- /dev/null +++ b/pygac/lac_klm.py @@ -0,0 +1,227 @@ +#!/usr/bin/python +# Copyright (c) 2014-2019 +# + +# Author(s): + +# Abhay Devasthale +# Sajid Pareeth +# Martin Raspaud +# Adam Dybbroe + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +"""Reader for LAC KLM data.""" + +import datetime +import logging + +import numpy as np + +from pygac.klm_reader import KLMReader +from pygac.lac_reader import LACReader + +LOG = logging.getLogger(__name__) + + +# video data object + +scanline = np.dtype([("scan_line_number", ">u2"), + ("scan_line_year", ">u2"), + ("scan_line_day_of_year", ">u2"), + ("satellite_clock_drift_delta", ">i2"), + ("scan_line_utc_time_of_day", ">u4"), + ("scan_line_bit_field", ">u2"), + ("zero_fill0", ">i2", (5, )), + # QUALITY INDICATORS + ("quality_indicator_bit_field", ">u4"), + ("scan_line_quality_flags", [("reserved", ">u1"), + ("time_problem_code", ">u1"), + ("calibration_problem_code", ">u1"), + ("earth_location_problem_code", ">u1")]), + ("calibration_quality_flags", ">u2", (3, )), + ("count_of_bit_errors_in_frame_sync", ">u2"), + ("zero_fill1", ">i4", (2, )), + # CALIBRATION COEFFICIENTS + ("visible_operational_cal_ch_1_slope_1", ">i4"), + ("visible_operational_cal_ch_1_intercept_1", ">i4"), + ("visible_operational_cal_ch_1_slope_2", ">i4"), + ("visible_operational_cal_ch_1_intercept_2", ">i4"), + ("visible_operational_cal_ch_1_intersection", ">i4"), + ("visible_test_cal_ch_1_slope_1", ">i4"), + ("visible_test_cal_ch_1_intercept_1", ">i4"), + ("visible_test_cal_ch_1_slope_2", ">i4"), + ("visible_test_cal_ch_1_intercept_2", ">i4"), + ("visible_test_cal_ch_1_intersection", ">i4"), + ("visible_prelaunch_cal_ch_1_slope_1", ">i4"), + ("visible_prelaunch_cal_ch_1_intercept_1", ">i4"), + ("visible_prelaunch_cal_ch_1_slope_2", ">i4"), + ("visible_prelaunch_cal_ch_1_intercept_2", ">i4"), + ("visible_prelaunch_cal_ch_1_intersection", ">i4"), + ("visible_operational_cal_ch_2_slope_1", ">i4"), + ("visible_operational_cal_ch_2_intercept_1", ">i4"), + ("visible_operational_cal_ch_2_slope_2", ">i4"), + ("visible_operational_cal_ch_2_intercept_2", ">i4"), + ("visible_operational_cal_ch_2_intersection", ">i4"), + ("visible_test_cal_ch_2_slope_1", ">i4"), + ("visible_test_cal_ch_2_intercept_1", ">i4"), + ("visible_test_cal_ch_2_slope_2", ">i4"), + ("visible_test_cal_ch_2_intercept_2", ">i4"), + ("visible_test_cal_ch_2_intersection", ">i4"), + ("visible_prelaunch_cal_ch_2_slope_1", ">i4"), + ("visible_prelaunch_cal_ch_2_intercept_1", ">i4"), + ("visible_prelaunch_cal_ch_2_slope_2", ">i4"), + ("visible_prelaunch_cal_ch_2_intercept_2", ">i4"), + ("visible_prelaunch_cal_ch_2_intersection", ">i4"), + ("visible_operational_cal_ch_3a_slope_1", ">i4"), + ("visible_operational_cal_ch_3a_intercept_1", ">i4"), + ("visible_operational_cal_ch_3a_slope_2", ">i4"), + ("visible_operational_cal_ch_3a_intercept_2", ">i4"), + ("visible_operational_cal_ch_3a_intersection", ">i4"), + ("visible_test_cal_ch_3a_slope_1", ">i4"), + ("visible_test_cal_ch_3a_intercept_1", ">i4"), + ("visible_test_cal_ch_3a_slope_2", ">i4"), + ("visible_test_cal_ch_3a_intercept_2", ">i4"), + ("visible_test_cal_ch_3a_intersection", ">i4"), + ("visible_prelaunch_cal_ch_3a_slope_1", ">i4"), + ("visible_prelaunch_cal_ch_3a_intercept_1", ">i4"), + ("visible_prelaunch_cal_ch_3a_slope_2", ">i4"), + ("visible_prelaunch_cal_ch_3a_intercept_2", ">i4"), + ("visible_prelaunch_cal_ch_3a_intersection", ">i4"), + ("ir_operational_cal_ch_3b_coefficient_1", ">i4"), + ("ir_operational_cal_ch_3b_coefficient_2", ">i4"), + ("ir_operational_cal_ch_3b_coefficient_3", ">i4"), + ("ir_test_cal_ch_3b_coefficient_1", ">i4"), + ("ir_test_cal_ch_3b_coefficient_2", ">i4"), + ("ir_test_cal_ch_3b_coefficient_3", ">i4"), + ("ir_operational_cal_ch_4_coefficient_1", ">i4"), + ("ir_operational_cal_ch_4_coefficient_2", ">i4"), + ("ir_operational_cal_ch_4_coefficient_3", ">i4"), + ("ir_test_cal_ch_4_coefficient_1", ">i4"), + ("ir_test_cal_ch_4_coefficient_2", ">i4"), + ("ir_test_cal_ch_4_coefficient_3", ">i4"), + ("ir_operational_cal_ch_5_coefficient_1", ">i4"), + ("ir_operational_cal_ch_5_coefficient_2", ">i4"), + ("ir_operational_cal_ch_5_coefficient_3", ">i4"), + ("ir_test_cal_ch_5_coefficient_1", ">i4"), + ("ir_test_cal_ch_5_coefficient_2", ">i4"), + ("ir_test_cal_ch_5_coefficient_3", ">i4"), + # NAVIGATION + ("computed_yaw_steering", ">i2", (3,)), # only in version 5 + ("total_applied_attitude_correction", + ">i2", (3,)), # only in version 5 + ("navigation_status_bit_field", ">u4"), + ("time_associated_with_tip_euler_angles", ">u4"), + ("tip_euler_angles", ">i2", (3, )), + ("spacecraft_altitude_above_reference_ellipsoid", ">u2"), + ("angular_relationships", ">i2", (153, )), + ("zero_fill2", ">i2", (3, )), + ("earth_location", ">i4", (102, )), + ("zero_fill3", ">i4", (2, )), + # HRPT MINOR FRAME TELEMETRY + ("frame_sync", ">u2", (6, )), + ("id", ">u2", (2, )), + ("time_code", ">u2", (4, )), + ('telemetry', [("ramp_calibration", '>u2', (5, )), + ("PRT", '>u2', (3, )), + ("ch3_patch_temp", '>u2'), + ("spare", '>u2'), ]), + ("back_scan", ">u2", (30, )), + ("space_data", ">u2", (50, )), + ("sync_delta", ">u2"), + ("zero_fill4", ">i2"), + # EARTH OBSERVATIONS + ("sensor_data", ">u4", (3414,)), + ("zero_fill5", ">i4", (2,)), + # DIGITAL B HOUSEKEEPING TELEMETRY + ("digital_b_telemetry_update_flags", ">u2"), + ("avhrr_digital_b_data", ">u2"), + ("zero_fill6", ">i4", (3,)), + # ANALOG HOUSEKEEPING DATA (TIP) + ("analog_telemetry_update_flags", ">u4"), + ("patch_temperature_range", ">u1"), + ("patch_temperature_extended", ">u1"), + ("patch_power", ">u1"), + ("radiator_temperature", ">u1"), + ("blackbody_temperature_1", ">u1"), + ("blackbody_temperature_2", ">u1"), + ("blackbody_temperature_3", ">u1"), + ("blackbody_temperature_4", ">u1"), + ("electronics_current", ">u1"), + ("motor_current", ">u1"), + ("earth_shield_position", ">u1"), + ("electronics_temperature", ">u1"), + ("cooler_housing_temperature", ">u1"), + ("baseplate_temperature", ">u1"), + ("motor_housing_temperature", ">u1"), + ("a/d_converter_temperature", ">u1"), + ("detector_#4_bias_voltage", ">u1"), + ("detector_#5_bias_voltage", ">u1"), + ("blackbody_temperature_channel3b", ">u1"), + ("blackbody_temperature_channel4", ">u1"), + ("blackbody_temperature_channel5", ">u1"), + ("reference_voltage", ">u1"), + ("zero_fill7", ">i2", (3,)), + # CLOUDS FROM AVHRR (CLAVR) + ("reserved_clavr_status_bit_field", ">u4"), + ("reserved_clavr", ">u4"), + ("reserved_clavr_ccm", ">u2", (256,)), + # FILLER + ("zero_fill8", ">i4", (94,))]) + + +class LACKLMReader(LACReader, KLMReader): + """The LAC KLM reader.""" + + def __init__(self, *args, **kwargs): + """Init the LAC KLM reader.""" + LACReader.__init__(self, *args, **kwargs) + self.scanline_type = scanline + self.offset = 15872 + # packed: 15872 + # self.offset = 22528 + # unpacked: 22528 + + +def main(filename, start_line, end_line): + """Generate a l1c file.""" + from pygac import gac_io + tic = datetime.datetime.now() + reader = LACKLMReader() + reader.read(filename) + reader.get_lonlat() + channels = reader.get_calibrated_channels() + sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles() + + mask, qual_flags = reader.get_corrupt_mask() + if (np.all(mask)): + print("ERROR: All data is masked out. Stop processing") + raise ValueError("All data is masked out.") + + gac_io.save_gac(reader.spacecraft_name, + reader.utcs, + reader.lats, reader.lons, + channels[:, :, 0], channels[:, :, 1], + channels[:, :, 2], + channels[:, :, 3], + channels[:, :, 4], + channels[:, :, 5], + sun_zen, sat_zen, sun_azi, sat_azi, rel_azi, + mask, qual_flags, start_line, end_line, reader.get_ch3_switch()) + LOG.info("pygac took: %s", str(datetime.datetime.now() - tic)) + + +if __name__ == "__main__": + import sys + main(sys.argv[1], sys.argv[2], sys.argv[3]) diff --git a/pygac/lac_pod.py b/pygac/lac_pod.py new file mode 100644 index 00000000..88f6ba32 --- /dev/null +++ b/pygac/lac_pod.py @@ -0,0 +1,93 @@ +#!/usr/bin/python +# Copyright (c) 2014-2019. +# + +# Author(s): + +# Abhay Devasthale +# Adam Dybbroe +# Sajid Pareeth +# Martin Raspaud + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +"""The LAC POD reader routines.""" + +import datetime +import logging + +import numpy as np + +from pygac import gac_io +from pygac.pod_reader import PODReader +from pygac.lac_reader import LACReader + +LOG = logging.getLogger(__name__) + +scanline = np.dtype([("scan_line_number", ">i2"), + ("time_code", ">u2", (3, )), + ("quality_indicators", ">u4"), + ("calibration_coefficients", ">i4", (10, )), + ("number_of_meaningful_zenith_angles_and_earth_location_appended", + ">u1"), + ("solar_zenith_angles", "i1", (51, )), + ("earth_location", ">i2", (102, )), + ("telemetry", ">u4", (35, )), + ("sensor_data", ">u4", (3414, )), + ("add_on_zenith", ">u2", (10, )), + ("clock_drift_delta", ">u2"), # only in newest version + ("spare3", "u2", (337, ))]) + + +class LACPODReader(LACReader, PODReader): + """The LAC POD reader.""" + + def __init__(self, *args, **kwargs): + """Init the LAC POD reader.""" + LACReader.__init__(self, *args, **kwargs) + self.scanline_type = scanline + self.offset = 14800 + self.scan_points = np.arange(2048) + + +def main(filename, start_line, end_line): + """Generate a l1c file.""" + tic = datetime.datetime.now() + reader = LACPODReader() + reader.read(filename) + reader.get_lonlat() + reader.adjust_clock_drift() + channels = reader.get_calibrated_channels() + sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles() + + mask, qual_flags = reader.get_corrupt_mask() + if (np.all(mask)): + print("ERROR: All data is masked out. Stop processing") + raise ValueError("All data is masked out.") + gac_io.save_gac(reader.spacecraft_name, + reader.utcs, + reader.lats, reader.lons, + channels[:, :, 0], channels[:, :, 1], + np.ones_like(channels[:, :, 0]) * -1, + channels[:, :, 2], + channels[:, :, 3], + channels[:, :, 4], + sun_zen, sat_zen, sun_azi, sat_azi, rel_azi, + mask, qual_flags, start_line, end_line) + LOG.info("pygac took: %s", str(datetime.datetime.now() - tic)) + + +if __name__ == "__main__": + import sys + main(sys.argv[1], sys.argv[2], sys.argv[3]) diff --git a/pygac/lac_reader.py b/pygac/lac_reader.py new file mode 100644 index 00000000..0907d8ba --- /dev/null +++ b/pygac/lac_reader.py @@ -0,0 +1,40 @@ +#!/usr/bin/python +# Copyright (c) 2014-2019. +# + +# Author(s): + +# Abhay Devasthale +# Adam Dybbroe +# Sajid Pareeth +# Martin Raspaud + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +"""The LAC reader.""" + +from pygac.reader import Reader +import pygac.pygac_geotiepoints as gtp + + +class LACReader(Reader): + """Reader for LAC data.""" + + # Scanning frequency (scanlines per millisecond) + scan_freq = 6.0 / 1000.0 + + def __init__(self, *args, **kwargs): + """Init the LAC reader.""" + super(LACReader, self).__init__(*args, **kwargs) + self.scan_width = 2048 + self.lonlat_interpolator = gtp.lac_lat_lon_interpolator diff --git a/pygac/pod_reader.py b/pygac/pod_reader.py new file mode 100644 index 00000000..ed597f10 --- /dev/null +++ b/pygac/pod_reader.py @@ -0,0 +1,491 @@ +#!/usr/bin/python +# Copyright (c) 2014-2019 Pygac developers +# + +# Author(s): + +# Abhay Devasthale +# Adam Dybbroe +# Sajid Pareeth +# Martin Raspaud + +# This work was done in the framework of ESA-CCI-Clouds phase I + + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + + +"""POD file reading. + +Format specification can be found here: +http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/c2/sec2-0.htm +http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/c3/sec3-1.htm +""" + +import datetime +import logging + +import numpy as np + +from pygac.correct_tsm_issue import TSM_AFFECTED_INTERVALS_POD, get_tsm_idx +from pygac.reader import Reader + +LOG = logging.getLogger(__name__) + +# common header +header0 = np.dtype([("noaa_spacecraft_identification_code", ">u1"), + ("data_type_code", ">u1"), + ("start_time", ">u2", (3, )), + ("number_of_scans", ">u2"), + ("end_time", ">u2", (3, ))]) + + +# For L1B data before September 8, 1992 +# http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/k/app-k.htm +header1 = np.dtype([("noaa_spacecraft_identification_code", ">u1"), + ("data_type_code", ">u1"), + ("start_time", ">u2", (3, )), + ("number_of_scans", ">u2"), + ("end_time", ">u2", (3, )), + ("processing_block_id", "S7"), + ("ramp_auto_calibration", ">u1"), + ("number_of_data_gaps", ">u2"), + ("dacs_quality", ">u1", (6, )), + ("calibration_parameter_id", ">i2"), + ("dacs_status", ">u1"), + ("spare1", ">i1", (5, )), + ("data_set_name", "S44")]) + +# For L1B data between October 21, 1992 to November 15, 1994 +# http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/l/app-l.htm +header2 = np.dtype([("noaa_spacecraft_identification_code", ">u1"), + ("data_type_code", ">u1"), + ("start_time", ">u2", (3, )), + ("number_of_scans", ">u2"), + ("end_time", ">u2", (3, )), + ("processing_block_id", "S7"), + ("ramp_auto_calibration", ">u1"), + ("number_of_data_gaps", ">u2"), + ("dacs_quality", ">u1", (6, )), + ("calibration_parameter_id", ">i2"), + ("dacs_status", ">u1"), + ("spare1", ">i1", (5, )), + ("data_set_name", "S42"), + ("blankfill", "S2"), + ("julian_year_of_epoch", ">u2"), + ("julian_day_of_epoch", ">u2"), + ("millisecond_utc_epoch_time_of_day", ">u4"), + # Keplerian orbital elements + ("semi_major_axis", ">f8"), + ("eccentricity", ">f8"), + ("inclination", ">f8"), + ("argument_of_perigee", ">f8"), + ("right_ascension", ">f8"), + ("mean_anomaly", ">f8"), + # cartesian inertial true date of elements + ("x_component_of_position_vector", ">f8"), + ("y_component_of_position_vector", ">f8"), + ("z_component_of_position_vector", ">f8"), + ("x_dot_component_of_position_vector", ">f8"), + ("y_dot_component_of_position_vector", ">f8"), + ("z_dot_component_of_position_vector", ">f8")]) + +# For L1B data post November 15, 1994 +# http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/c2/sec2-0.htm +header3 = np.dtype([("noaa_spacecraft_identification_code", ">u1"), + ("data_type_code", ">u1"), + ("start_time", ">u2", (3, )), + ("number_of_scans", ">u2"), + ("end_time", ">u2", (3, )), + ("processing_block_id", "S7"), + ("ramp_auto_calibration", ">u1"), + ("number_of_data_gaps", ">u2"), + ("dacs_quality", ">u1", (6, )), + ("calibration_parameter_id", ">i2"), + ("dacs_status", ">u1"), + ("reserved_for_mounting_and_fixed_attitude_correction_indicator", + ">i1"), + ("nadir_earth_location_tolerance", ">i1"), + ("spare1", ">i1"), + ("start_of_data_set_year", ">u2"), + ("data_set_name", "S44"), + ("year_of_epoch", ">u2"), + ("julian_day_of_epoch", ">u2"), + ("millisecond_utc_epoch_time_of_day", ">u4"), + # Keplerian orbital elements + ("semi_major_axis", ">i4"), + ("eccentricity", ">i4"), + ("inclination", ">i4"), + ("argument_of_perigee", ">i4"), + ("right_ascension", ">i4"), + ("mean_anomaly", ">i4"), + # cartesian inertial true date of elements + ("x_component_of_position_vector", ">i4"), + ("y_component_of_position_vector", ">i4"), + ("z_component_of_position_vector", ">i4"), + ("x_dot_component_of_position_vector", ">i4"), + ("y_dot_component_of_position_vector", ">i4"), + ("z_dot_component_of_position_vector", ">i4"), + # future use + ("yaw_fixed_error_correction", ">i2"), + ("roll_fixed_error_correction", ">i2"), + ("pitch_fixed_error_correction", ">i2")]) + +# archive header +tbm_header = np.dtype([('fill', 'S30'), + ('data_set_name', 'S44'), + ('select_flag', 'S1'), + ('beginning_latitude', 'S3'), + ('ending_latitude', 'S3'), + ('beginning_longitude', 'S4'), + ('ending_longitude', 'S4'), + ('start_hour', 'S2'), + ('start_minute', 'S2'), + ('number_of_minutes', 'S3'), + ('appended_data_flag', 'S1'), + ('channel_select_flag', 'S1', (20, )), + ('sensor_data_word_size', 'S2'), + ('fill2', 'S3')]) + + +class PODReader(Reader): + """The POD reader.""" + + spacecrafts_orbital = {25: 'tiros n', + 2: 'noaa 6', + 4: 'noaa 7', + 6: 'noaa 8', + 7: 'noaa 9', + 8: 'noaa 10', + 1: 'noaa 11', + 5: 'noaa 12', + 3: 'noaa 14', + } + spacecraft_names = {25: 'tirosn', + 2: 'noaa6', + 4: 'noaa7', + 6: 'noaa8', + 7: 'noaa9', + 8: 'noaa10', + 1: 'noaa11', + 5: 'noaa12', + 3: 'noaa14', + } + + tsm_affected_intervals = TSM_AFFECTED_INTERVALS_POD + + def correct_scan_line_numbers(self): + """Correct the scan line numbers.""" + # Perform common corrections first. + super(PODReader, self).correct_scan_line_numbers() + + # cleaning up the data + min_scanline_number = np.amin(np.absolute(self.scans["scan_line_number"][:])) + if self.scans["scan_line_number"][0] == self.scans["scan_line_number"][-1] + 1: + while self.scans["scan_line_number"][0] != min_scanline_number: + self.scans = np.roll(self.scans, -1) + else: + while self.scans["scan_line_number"][0] != min_scanline_number: + self.scans = self.scans[1:] + + self.scans = self.scans[self.scans["scan_line_number"] != 0] + + def read(self, filename): + """Read the data. + + Args: + filename: string + The filename to read from. + + Returns: + header: numpy record array + The header metadata + scans: numpy record array + The scanlines + + """ + super(PODReader, self).read(filename=filename) + # choose the right header depending on the date + with open(filename) as fd_: + # read archive header + self.tbm_head = np.fromfile(fd_, dtype=tbm_header, count=1)[0] + if ((not self.tbm_head['data_set_name'].startswith(self.creation_site + b'.')) and + (self.tbm_head['data_set_name'] != b'\x00' * 42 + b' ')): + fd_.seek(0) + self.tbm_head = None + tbm_offset = 0 + else: + tbm_offset = tbm_header.itemsize + + head = np.fromfile(fd_, dtype=header0, count=1)[0] + year, jday, _ = self.decode_timestamps(head["start_time"]) + + start_date = (datetime.date(year, 1, 1) + datetime.timedelta(days=int(jday) - 1)) + + if start_date < datetime.date(1992, 9, 8): + header = header1 + elif start_date <= datetime.date(1994, 11, 15): + header = header2 + else: + header = header3 + + fd_.seek(tbm_offset, 0) + self.head = np.fromfile(fd_, dtype=header, count=1)[0] + fd_.seek(self.offset + tbm_offset, 0) + self.scans = np.fromfile(fd_, + dtype=self.scanline_type, + count=self.head["number_of_scans"]) + + self.correct_scan_line_numbers() + self.spacecraft_id = self.head["noaa_spacecraft_identification_code"] + if self.spacecraft_id == 1 and start_date < datetime.date(1982, 1, 1): + self.spacecraft_id = 25 + self.spacecraft_name = self.spacecraft_names[self.spacecraft_id] + LOG.info( + "Reading %s data", self.spacecrafts_orbital[self.spacecraft_id]) + + return self.head, self.scans + + def get_header_timestamp(self): + """Get the timestamp from the header. + + Returns: + A datetime object containing the timestamp from the header. + + Raises: + A ValueError if the timestamp is corrupt. + + """ + year, jday, msec = self.decode_timestamps(self.head["start_time"]) + try: + return self.to_datetime(self.to_datetime64(year=year, jday=jday, + msec=msec)) + except ValueError as err: + raise ValueError('Corrupt header timestamp: {0}'.format(err)) + + @staticmethod + def decode_timestamps(encoded): + """Decode timestamps. + + Returns: + year + day of year + milliseconds since 00:00 + + """ + ndims = len(encoded.shape) + if ndims == 1: + # Single header timestamp + enc0 = encoded[0] + enc1 = encoded[1] + enc2 = encoded[2] + elif ndims == 2: + # Scanline timestamps + enc0 = encoded[:, 0] + enc1 = encoded[:, 1] + enc2 = encoded[:, 2] + else: + raise ValueError('Invalid timestamp dimension') + + year = enc0 >> 9 + year = np.where(year > 75, year + 1900, year + 2000) + jday = (enc0 & 0x1FF) + msec = ((np.uint32(enc1 & 2047) << 16) | (np.uint32(enc2))) + + return year, jday, msec + + def _get_times(self): + return self.decode_timestamps(self.scans["time_code"]) + + def _adjust_clock_drift(self): + """Adjust the geolocation to compensate for the clock error. + + TODO: bad things might happen when scanlines are skipped. + """ + tic = datetime.datetime.now() + self.get_times() + from pygac.clock_offsets_converter import get_offsets + try: + offset_times, clock_error = get_offsets(self.spacecraft_name) + except KeyError: + LOG.info("No clock drift info available for %s", + self.spacecraft_name) + else: + offset_times = np.array(offset_times, dtype='datetime64[ms]') + offsets = np.interp(self.utcs.astype(np.uint64), + offset_times.astype(np.uint64), + clock_error) + LOG.info("Adjusting for clock drift of %s to %s", + str(min(offsets)), + str(max(offsets))) + self.times = (self.utcs + + offsets.astype('timedelta64[s]')).astype(datetime.datetime) + offsets *= -2 + + int_offsets = np.floor(offsets).astype(np.int) + + # filling out missing geolocations with computation from pyorbital. + line_indices = (self.scans["scan_line_number"] + + int_offsets) + + missed = sorted((set(line_indices) | + set(line_indices + 1)) + - set(self.scans["scan_line_number"])) + + min_idx = min(line_indices) + max_idx = max(max(line_indices), + max(self.scans["scan_line_number"] - min_idx)) + 1 + idx_len = max_idx - min_idx + 2 + + complete_lons = np.full((idx_len, self.lats.shape[1]), np.nan, + dtype=np.float) + complete_lats = np.full((idx_len, self.lats.shape[1]), np.nan, + dtype=np.float) + + complete_lons[self.scans["scan_line_number"] - min_idx] = self.lons + complete_lats[self.scans["scan_line_number"] - min_idx] = self.lats + + missed_utcs = ((np.array(missed) - 1) * np.timedelta64(500, "ms") + + self.utcs[0]) + try: + mlons, mlats = self.compute_lonlat(width=self.lats.shape[1], + utcs=missed_utcs, + clock_drift_adjust=True) + except IndexError as err: + LOG.warning('Cannot perform clock drift correction: %s', str(err)) + return + + complete_lons[missed - min_idx] = mlons + complete_lats[missed - min_idx] = mlats + + from pygac.slerp import slerp + off = offsets - np.floor(offsets) + res = slerp(complete_lons[line_indices - min_idx, :], + complete_lats[line_indices - min_idx, :], + complete_lons[line_indices - min_idx + 1, :], + complete_lats[line_indices - min_idx + 1, :], + off[:, np.newaxis, np.newaxis]) + + self.lons = res[:, :, 0] + self.lats = res[:, :, 1] + self.utcs += offsets.astype('timedelta64[s]') + + toc = datetime.datetime.now() + LOG.debug("clock drift adjustment took %s", str(toc - tic)) + + def _get_lonlat(self): + lats = self.scans["earth_location"][:, 0::2] / 128.0 + lons = self.scans["earth_location"][:, 1::2] / 128.0 + return lons, lats + + def get_telemetry(self): + """Get the telemetry. + + Returns: + prt_counts: np.array + ict_counts: np.array + space_counts: np.array + + """ + number_of_scans = self.scans["telemetry"].shape[0] + decode_tele = np.zeros((int(number_of_scans), 105)) + decode_tele[:, ::3] = (self.scans["telemetry"] >> 20) & 1023 + decode_tele[:, 1::3] = (self.scans["telemetry"] >> 10) & 1023 + decode_tele[:, 2::3] = self.scans["telemetry"] & 1023 + + prt_counts = np.mean(decode_tele[:, 17:20], axis=1) + + # getting ICT counts + + ict_counts = np.zeros((int(number_of_scans), 3)) + ict_counts[:, 0] = np.mean(decode_tele[:, 22:50:3], axis=1) + ict_counts[:, 1] = np.mean(decode_tele[:, 23:51:3], axis=1) + ict_counts[:, 2] = np.mean(decode_tele[:, 24:52:3], axis=1) + + # getting space counts + + space_counts = np.zeros((int(number_of_scans), 3)) + space_counts[:, 0] = np.mean(decode_tele[:, 54:100:5], axis=1) + space_counts[:, 1] = np.mean(decode_tele[:, 55:101:5], axis=1) + space_counts[:, 2] = np.mean(decode_tele[:, 56:102:5], axis=1) + + return prt_counts, ict_counts, space_counts + + def _get_corrupt_mask(self): + """Get mask for corrupt scanlines.""" + mask = ((self.scans["quality_indicators"] >> 31) | + ((self.scans["quality_indicators"] << 4) >> 31) | + ((self.scans["quality_indicators"] << 5) >> 31)) + return mask.astype(bool) + + def get_qual_flags(self): + """Read quality flags.""" + number_of_scans = self.scans["telemetry"].shape[0] + qual_flags = np.zeros((int(number_of_scans), 7)) + qual_flags[:, 0] = self.scans["scan_line_number"] + qual_flags[:, 1] = (self.scans["quality_indicators"] >> 31) + qual_flags[:, 2] = ((self.scans["quality_indicators"] << 4) >> 31) + qual_flags[:, 3] = ((self.scans["quality_indicators"] << 5) >> 31) + qual_flags[:, 4] = ((self.scans["quality_indicators"] << 13) >> 31) + qual_flags[:, 5] = ((self.scans["quality_indicators"] << 14) >> 31) + qual_flags[:, 6] = ((self.scans["quality_indicators"] << 15) >> 31) + + return qual_flags + + def postproc(self, channels): + """No POD specific postprocessing to be done.""" + pass + + def get_tsm_pixels(self, channels): + """Determine pixels affected by the scan motor issue. + + Uses channels 1, 2, 4 and 5. Neither 3a, nor 3b. + """ + return get_tsm_idx(channels[:, :, 0], channels[:, :, 1], + channels[:, :, 3], channels[:, :, 4]) + + +def main(filename, start_line, end_line): + """Generate a l1c file.""" + from pygac import gac_io + tic = datetime.datetime.now() + reader = PODReader() + reader.read(filename) + reader.get_lonlat() + channels = reader.get_calibrated_channels() + sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles() + qual_flags = reader.get_qual_flags() + if np.all(reader.mask): + print("ERROR: All data is masked out. Stop processing") + raise ValueError("All data is masked out.") + + gac_io.save_gac(reader.spacecraft_name, + reader.utcs, + reader.lats, reader.lons, + channels[:, :, 0], channels[:, :, 1], + np.full_like(channels[:, :, 0], np.nan), + channels[:, :, 2], + channels[:, :, 3], + channels[:, :, 4], + sun_zen, sat_zen, sun_azi, sat_azi, rel_azi, + qual_flags, start_line, end_line, + reader.filename, + reader.get_midnight_scanline(), + reader.get_miss_lines()) + LOG.info("pygac took: %s", str(datetime.datetime.now() - tic)) + + +if __name__ == "__main__": + import sys + main(sys.argv[1], sys.argv[2], sys.argv[3]) diff --git a/pygac/pygac_geotiepoints.py b/pygac/pygac_geotiepoints.py index 974e1256..96e5ffb5 100644 --- a/pygac/pygac_geotiepoints.py +++ b/pygac/pygac_geotiepoints.py @@ -29,7 +29,7 @@ import numpy as np -def Gac_Lat_Lon_Interpolator(lons_subset, lats_subset): +def gac_lat_lon_interpolator(lons_subset, lats_subset): """Interpolate lat-lon values in the AVHRR GAC data. Each GAC row has 409 pixels. @@ -37,10 +37,20 @@ def Gac_Lat_Lon_Interpolator(lons_subset, lats_subset): ranging from 5 to 405. Interpolate to full resolution. """ - # cols_subset = np.arange(0, 404, 8) - # cols_full = np.arange(405) cols_subset = np.arange(4, 405, 8) cols_full = np.arange(409) + return lat_lon_interpolator(lons_subset, lats_subset, cols_subset, cols_full) + + +def lac_lat_lon_interpolator(lons_subset, lats_subset): + """Interpolate lat-lon values in the AVHRR LAC data.""" + cols_subset = np.arange(24, 2048, 40) + cols_full = np.arange(2048) + return lat_lon_interpolator(lons_subset, lats_subset, cols_subset, cols_full) + + +def lat_lon_interpolator(lons_subset, lats_subset, cols_subset, cols_full): + """Interpolate lat-lon values in the AVHRR data.""" lines = lats_subset.shape[0] rows_subset = np.arange(lines) rows_full = np.arange(lines) diff --git a/pygac/reader.py b/pygac/reader.py new file mode 100644 index 00000000..4083fe29 --- /dev/null +++ b/pygac/reader.py @@ -0,0 +1,891 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +# Copyright (c) 2014, 2019 Pygac Developers + +# Author(s): + +# Martin Raspaud + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +"""Generic reader for GAC and LAC data. + +Can't be used as is, has to be subclassed to add specific read functions. +""" +from abc import ABCMeta, abstractmethod, abstractproperty +import logging +import numpy as np +import os +import six +import types + +from pygac import (CONFIG_FILE, centered_modulus, + get_absolute_azimuth_angle_diff) +try: + import ConfigParser +except ImportError: + import configparser as ConfigParser + +from pyorbital.orbital import Orbital +from pyorbital import astronomy +import datetime +from pygac.calibration import calibrate_solar, calibrate_thermal + +LOG = logging.getLogger(__name__) + +# rpy values from +# here:http://yyy.rsmas.miami.edu/groups/rrsl/pathfinder/Processing/proc_app_a.html +rpy_coeffs = { + 'noaa7': {'roll': 0.000, + 'pitch': 0.000, + 'yaw': 0.000, + }, + 'noaa9': {'roll': 0.000, + 'pitch': 0.0025, + 'yaw': 0.000, + }, + 'noaa10': {'roll': 0.000, + 'pitch': 0.000, + 'yaw': 0.000, + }, + 'noaa11': {'roll': -0.0019, + 'pitch': -0.0037, + 'yaw': 0.000, + }, + 'noaa12': {'roll': 0.000, + 'pitch': 0.000, + 'yaw': 0.000, + }, + 'noaa14': {'roll': 0.000, + 'pitch': 0.000, + 'yaw': 0.000, + }} + + +class Reader(six.with_metaclass(ABCMeta)): + """Reader for Gac and Lac, POD and KLM data.""" + + def __init__(self, interpolate_coords=True, adjust_clock_drift=True, + tle_dir=None, tle_name=None, tle_thresh=7, creation_site=None): + """Init the reader. + + Args: + interpolate_coords: Interpolate coordinates from tiepoint grid + to all pixels. + adjust_clock_drift: Adjust the geolocation to compensate for the + clock error (POD satellites only). + tle_dir: Directory holding TLE files + tle_name: Filename pattern of TLE files. + tle_thresh: Maximum number of days between observation and nearest + TLE + creation_site: The three-letter identifier of the creation site (eg 'NSS') + + """ + self.interpolate_coords = interpolate_coords + self.adjust_clock_drift = adjust_clock_drift + self.tle_dir = tle_dir + self.tle_name = tle_name + self.tle_thresh = tle_thresh + self.creation_site = (creation_site or 'NSS').encode('utf-8') + self.head = None + self.scans = None + self.spacecraft_name = None + self.spacecraft_id = None + self.utcs = None + self.lats = None + self.lons = None + self.times = None + self.tle_lines = None + self.filename = None + self._mask = None + + @abstractmethod + def read(self, filename): + """Read the GAC/LAC data. + + Args: + filename (str): Specifies the GAC/LAC file to be read. + """ + self.filename = os.path.basename(filename) + LOG.info('Reading %s', self.filename) + + @abstractmethod + def get_header_timestamp(self): + """Read start timestamp from the header. + + Returns: + datetime.datetime: Start timestamp + + """ + raise NotImplementedError + + def get_counts(self): + """Get the counts. + + Returns: + np.array: The counts, with channel 3a and 3b split if necessary. + + """ + packed_data = self.scans["sensor_data"] + counts = np.zeros((len(self.scans), self.scan_width * 5)) + counts_nb = (self.scan_width * 5) // 3 + remainder = (self.scan_width * 5) % 3 + if remainder == 0: + nb1 = nb2 = nb3 = counts_nb + elif remainder == 1: + nb1 = counts_nb + 1 + nb2 = nb3 = counts_nb + elif remainder == 2: + nb1 = nb2 = counts_nb + 1 + nb3 = counts_nb + counts[:, 0::3] = ((packed_data >> 20) & 1023)[:, :nb1] + counts[:, 1::3] = ((packed_data >> 10) & 1023)[:, :nb2] + counts[:, 2::3] = (packed_data & 1023)[:, :nb3] + counts = counts.reshape((-1, self.scan_width, 5)) + try: + switch = self.get_ch3_switch() + except AttributeError: + return counts + else: + channels = np.zeros((len(self.scans), self.scan_width, 6), + dtype=counts.dtype) + channels[:, :, :2] = counts[:, :, :2] + channels[:, :, -2:] = counts[:, :, -2:] + channels[:, :, 2][switch == 1] = counts[:, :, 2][switch == 1] + channels[:, :, 3][switch == 0] = counts[:, :, 2][switch == 0] + return channels + + @abstractmethod + def _get_times(self): + """Specify how to read scanline timestamps from GAC data. + + Returns: + int: year + int: day of year + int: milliseconds since 00:00 + + """ + raise NotImplementedError + + def get_times(self): + """Read scanline timestamps and try to correct invalid values. + + Note: + Also sets self.utcs and self.times! + + Returns: + UTC timestamps + + """ + if self.utcs is None: + # Read timestamps + year, jday, msec = self._get_times() + + # Correct invalid values + year, jday, msec = self.correct_times_median(year=year, jday=jday, + msec=msec) + self.utcs = self.to_datetime64(year=year, jday=jday, msec=msec) + self.correct_times_thresh() + + # Convert timestamps to datetime objects + self.times = self.to_datetime(self.utcs) + + return self.utcs + + @staticmethod + def to_datetime64(year, jday, msec): + """Convert timestamps to numpy.datetime64. + + Args: + year: Year + jday: Day of the year (1-based) + msec: Milliseconds since 00:00 + + Returns: + numpy.datetime64: Converted timestamps + + """ + return (((year - 1970).astype('datetime64[Y]') + + (jday - 1).astype('timedelta64[D]')).astype('datetime64[ms]') + + msec.astype('timedelta64[ms]')) + + @staticmethod + def to_datetime(datetime64): + """Convert numpy.datetime64 to datetime.datetime. + + Args: + datetime64 (numpy.datetime64): Numpy timestamp to be converted. + + Returns: + datetime.datetime: Converted timestamp + + """ + return datetime64.astype(datetime.datetime) + + def lineno2msec(self, scan_line_number): + """Compute ideal scanline timestamp based on the scanline number. + + Assumes a constant scanning frequency. + + Args: + scan_line_number: Specifies the scanline number (1-based) + + Returns: + Corresponding timestamps in milliseconds since 1970-01-01 00:00, + i.e. the first scanline has timestamp 0. + + """ + return (scan_line_number - 1) / self.scan_freq + + def compute_lonlat(self, width, utcs=None, clock_drift_adjust=True): + """Compute lat/lon coordinates. + + Args: + width: Number of coordinates per scanlines + utcs: Scanline timestamps + clock_drift_adjust: If True, adjust clock drift. + + """ + if utcs is None: + utcs = self.get_times() + + # adjusting clock for drift + tic = datetime.datetime.now() + if clock_drift_adjust: + from pygac.clock_offsets_converter import get_offsets + try: + offset_times, clock_error = get_offsets(self.spacecraft_name) + except KeyError: + LOG.info("No clock drift info available for %s", + self.spacecraft_name) + else: + offset_times = np.array(offset_times, dtype='datetime64[ms]') + offsets = np.interp(utcs.astype(np.uint64), + offset_times.astype(np.uint64), + clock_error) + utcs = utcs - (offsets * 1000).astype('timedelta64[ms]') + + t = utcs[0].astype(datetime.datetime) + + if "constant_yaw_attitude_error" in self.head.dtype.fields: + rpy = np.deg2rad([self.head["constant_roll_attitude_error"] / 1e3, + self.head["constant_pitch_attitude_error"] / 1e3, + self.head["constant_yaw_attitude_error"] / 1e3]) + else: + try: + rpy_spacecraft = rpy_coeffs[self.spacecraft_name] + rpy = [rpy_spacecraft['roll'], + rpy_spacecraft['pitch'], + rpy_spacecraft['yaw']] + LOG.debug("Using static attitude correction") + except KeyError: + LOG.debug("Not applying attitude correction") + rpy = [0, 0, 0] + + LOG.info("Using rpy: %s", str(rpy)) + + from pyorbital.geoloc_instrument_definitions import avhrr_gac + from pyorbital.geoloc import compute_pixels, get_lonlatalt + # TODO: Are we sure all satellites have this scan width in degrees ? + sgeom = avhrr_gac(utcs.astype(datetime.datetime), self.scan_points, 55.385) + s_times = sgeom.times(t) + tle1, tle2 = self.get_tle_lines() + + pixels_pos = compute_pixels((tle1, tle2), sgeom, s_times, rpy) + pos_time = get_lonlatalt(pixels_pos, s_times) + + toc = datetime.datetime.now() + + LOG.warning("Computation of geolocation: %s", str(toc - tic)) + + lons, lats = pos_time[:2] + + return lons.reshape(-1, width), lats.reshape(-1, width) + + def get_calibrated_channels(self): + """Calibrate and return the channels.""" + channels = self.get_counts() + self.get_times() + year = self.times[0].year + delta = self.times[0].date() - datetime.date(year, 1, 1) + jday = delta.days + 1 + + # Earth-Sun distance correction factor + corr = 1.0 - 0.0334 * np.cos(2.0 * np.pi * (jday - 2) / 365.25) + + # how many reflective channels are there ? + tot_ref = channels.shape[2] - 3 + + channels[:, :, 0:tot_ref] = calibrate_solar(channels[:, :, 0:tot_ref], + np.arange(tot_ref), + year, jday, + self.spacecraft_name, + corr) + prt, ict, space = self.get_telemetry() + for chan in [3, 4, 5]: + channels[:, :, chan - 6] = calibrate_thermal( + channels[:, :, chan - 6], + prt, + ict[:, chan - 3], + space[:, chan - 3], + self.scans["scan_line_number"], + chan, + self.spacecraft_name) + + # Mask out corrupt values + channels[self.mask] = np.nan + + # Apply KLM/POD specific postprocessing + self.postproc(channels) + + # Mask pixels affected by scan motor issue + if self.is_tsm_affected(): + LOG.info('Correcting for temporary scan motor issue') + self.mask_tsm_pixels(channels) + + return channels + + def get_lonlat(self): + """Compute lat/lon coordinates. + + TODO: Switch to faster interpolator? + """ + if self.lons is None and self.lats is None: + self.lons, self.lats = self._get_lonlat() + + # Interpolate from every eighth pixel to all pixels. + if self.interpolate_coords: + self.lons, self.lats = self.lonlat_interpolator( + self.lons, self.lats) + + # Adjust clock drift + if self.adjust_clock_drift: + self._adjust_clock_drift() + + # Mask out corrupt scanlines + self.lons[self.mask] = np.nan + self.lats[self.mask] = np.nan + + # Mask values outside the valid range + self.lats[np.fabs(self.lats) > 90.0] = np.nan + self.lons[np.fabs(self.lons) > 180.0] = np.nan + + return self.lons, self.lats + + @abstractmethod + def _get_lonlat(self): + """KLM/POD specific readout of lat/lon coordinates.""" + raise NotImplementedError + + @property + def mask(self): + """Mask for corrupt scanlines.""" + if self._mask is None: + self._mask = self._get_corrupt_mask() + return self._mask + + @abstractmethod + def _get_corrupt_mask(self): + """KLM/POD specific readout of corrupt scanline mask.""" + raise NotImplementedError + + @abstractmethod + def postproc(self, channels): + """Apply KLM/POD specific postprocessing.""" + raise NotImplementedError + + @abstractmethod + def _adjust_clock_drift(self): + """Adjust clock drift.""" + raise NotImplementedError + + @staticmethod + def tle2datetime64(times): + """Convert TLE timestamps to numpy.datetime64. + + Args: + times (float): TLE timestamps as %y%j.1234, e.g. 18001.25 + """ + # Convert %y%j.12345 to %Y%j.12345 (valid for 1950-2049) + times = np.where(times > 50000, times + 1900000, times + 2000000) + + # Convert float to datetime64 + doys = (times % 1000).astype('int') - 1 + years = (times // 1000).astype('int') + msecs = np.rint(24 * 3600 * 1000 * (times % 1)) + times64 = (years - 1970).astype('datetime64[Y]').astype('datetime64[ms]') + times64 += doys.astype('timedelta64[D]') + times64 += msecs.astype('timedelta64[ms]') + + return times64 + + def get_tle_file(self): + """Find TLE file for the current satellite.""" + tle_dir, tle_name = self.tle_dir, self.tle_name + + # If user didn't specify TLE dir/name, try config file + if tle_dir is None or tle_name is None: + conf = ConfigParser.ConfigParser() + try: + conf.read(CONFIG_FILE) + except ConfigParser.NoSectionError: + LOG.exception('Failed reading configuration file: %s', + str(CONFIG_FILE)) + raise + + options = {} + for option, value in conf.items('tle', raw=True): + options[option] = value + + tle_dir = options['tledir'] + tle_name = options['tlename'] + + values = {"satname": self.spacecraft_name, } + tle_filename = os.path.join(tle_dir, tle_name % values) + LOG.info('TLE filename = ' + str(tle_filename)) + + return tle_filename + + def read_tle_file(self, tle_filename): + """Read TLE file.""" + with open(tle_filename, 'r') as fp_: + return fp_.readlines() + + def get_tle_lines(self): + """Find closest two line elements (TLEs) for the current orbit. + + Raises: + IndexError, if the closest TLE is more than :meth:`pygac.GACReader.tle_thresh` days apart + + """ + if self.tle_lines is not None: + return self.tle_lines + self.get_times() + tle_data = self.read_tle_file(self.get_tle_file()) + sdate = np.datetime64(self.times[0], '[ms]') + dates = self.tle2datetime64( + np.array([float(line[18:32]) for line in tle_data[::2]])) + + # Find index "iindex" such that dates[iindex-1] < sdate <= dates[iindex] + # Notes: + # 1. If sdate < dates[0] then iindex = 0 + # 2. If sdate > dates[-1] then iindex = len(dates), beyond the right boundary! + iindex = np.searchsorted(dates, sdate) + + if iindex in (0, len(dates)): + if iindex == len(dates): + # Reset index if beyond the right boundary (see note 2. above) + iindex -= 1 + elif abs(sdate - dates[iindex - 1]) < abs(sdate - dates[iindex]): + # Choose the closest of the two surrounding dates + iindex -= 1 + + # Make sure the TLE we found is within the threshold + delta_days = abs(sdate - dates[iindex]) / np.timedelta64(1, 'D') + if delta_days > self.tle_thresh: + raise IndexError( + "Can't find tle data for %s within +/- %d days around %s" % + (self.spacecraft_name, self.tle_thresh, sdate)) + + if delta_days > 3: + LOG.warning("Found TLE data for %s that is %f days appart", + sdate, delta_days) + else: + LOG.debug("Found TLE data for %s that is %f days appart", + sdate, delta_days) + + # Select TLE data + tle1 = tle_data[iindex * 2] + tle2 = tle_data[iindex * 2 + 1] + self.tle_lines = tle1, tle2 + return tle1, tle2 + + def get_angles(self): + """Get azimuth and zenith angles. + + Azimuth angle definition is the same as in pyorbital, but with + different units (degrees not radians for sun azimuth angles) + and different ranges. + + Returns: + sat_azi: satellite azimuth angle + degree clockwise from north in range ]-180, 180], + sat_zentih: satellite zenith angles in degrees in range [0,90], + sun_azi: sun azimuth angle + degree clockwise from north in range ]-180, 180], + sun_zentih: sun zenith angles in degrees in range [0,90], + rel_azi: absolute azimuth angle difference in degrees between sun + and sensor in range [0, 180] + + """ + self.get_times() + self.get_lonlat() + tle1, tle2 = self.get_tle_lines() + orb = Orbital(self.spacecrafts_orbital[self.spacecraft_id], + line1=tle1, line2=tle2) + + sat_azi, sat_elev = orb.get_observer_look(self.times[:, np.newaxis], + self.lons, self.lats, 0) + + sat_zenith = 90 - sat_elev + + sun_zenith = astronomy.sun_zenith_angle(self.times[:, np.newaxis], + self.lons, self.lats) + + alt, sun_azi = astronomy.get_alt_az(self.times[:, np.newaxis], + self.lons, self.lats) + del alt + sun_azi = np.rad2deg(sun_azi) + rel_azi = get_absolute_azimuth_angle_diff(sun_azi, sat_azi) + + # Scale angles range to half open interval ]-180, 180] + sat_azi = centered_modulus(sat_azi, 360.0) + sun_azi = centered_modulus(sun_azi, 360.0) + + # Mask corrupt scanlines + for arr in (sat_azi, sat_zenith, sun_azi, sun_zenith, rel_azi): + arr[self.mask] = np.nan + + return sat_azi, sat_zenith, sun_azi, sun_zenith, rel_azi + + def correct_times_median(self, year, jday, msec): + """Replace invalid timestamps with statistical estimates (using median). + + Assumes that the majority of timestamps is ok. + + Args: + year: Year + jday: Day of the year + msec: Milliseconds since 00:00 + + Returns: + Corrected year + Corrected day of the year + Corrected milliseconds + + """ + # Estimate ideal timestamps based on the scanline number. Still without + # offset, e.g. the first scanline has timestamp 1970-01-01 00:00 + msec_lineno = self.lineno2msec(self.scans["scan_line_number"]) + + jday = np.where(np.logical_or(jday < 1, jday > 366), np.median(jday), jday) + if_wrong_jday = np.ediff1d(jday, to_begin=0) + jday = np.where(if_wrong_jday < 0, max(jday), jday) + + if_wrong_msec = np.where(msec < 1) + if_wrong_msec = if_wrong_msec[0] + if len(if_wrong_msec) > 0: + if if_wrong_msec[0] != 0: + msec = msec[0] + msec_lineno + else: + msec0 = np.median(msec - msec_lineno) + msec = msec0 + msec_lineno + + if_wrong_msec = np.ediff1d(msec, to_begin=0) + msec = np.where(np.logical_and(np.logical_or(if_wrong_msec < -1000, if_wrong_msec > 1000), if_wrong_jday != 1), + msec[0] + msec_lineno, msec) + + # checking if year value is out of valid range + if_wrong_year = np.where( + np.logical_or(year < 1978, year > datetime.datetime.now().year)) + if_wrong_year = if_wrong_year[0] + if len(if_wrong_year) > 0: + # if the first scanline has valid time stamp + if if_wrong_year[0] != 0: + year = year[0] + jday = jday[0] + msec = msec[0] + msec_lineno + # Otherwise use median time stamp + else: + year = np.median(year) + jday = np.median(jday) + msec0 = np.median(msec - msec_lineno) + msec = msec0 + msec_lineno + + return year, jday, msec + + def correct_scan_line_numbers(self): + """Remove scanlines with corrupted scanline numbers. + + This includes: + - Scanline numbers outside the valide range + - Scanline numbers deviating more than a certain threshold from the + ideal case (1,2,3,...N) + + Example files having corrupt scanline numbers: + - NSS.GHRR.NJ.D96144.S2000.E2148.B0720102.GC + - NSS.GHRR.NJ.D96064.S0043.E0236.B0606162.WI + - NSS.GHRR.NJ.D99286.S1818.E2001.B2466869.WI + + Returns: + Intermediate and final results (for plotting purpose) + + """ + along_track = np.arange(1, len(self.scans["scan_line_number"])+1) + results = {'along_track': along_track, + 'n_orig': self.scans['scan_line_number'].copy()} + + # Remove scanlines whose scanline number is outside the valid range + within_range = np.logical_and(self.scans["scan_line_number"] < 15000, + self.scans["scan_line_number"] >= 0) + self.scans = self.scans[within_range] + + # Remove scanlines deviating more than a certain threshold from the + # ideal case (1,2,3,...N). + ideal = np.arange(1, len(self.scans["scan_line_number"])+1) + + # ... Estimate possible offset (in case some scanlines are missing in + # the beginning of the scan) + offsets = self.scans["scan_line_number"] - ideal + med_offset = np.median(offsets) + + # ... Compute difference to ideal case (1,2,3,...N) + offset + diffs = np.abs(self.scans["scan_line_number"] - (ideal + med_offset)) + + # ... Remove those scanlines whose difference is larger than a certain + # threshold. For the threshold computation we only regard nonzero + # differences. + nz_diffs = diffs[diffs > 0] + if len(nz_diffs) < 50: + # Not enough differences for reliable statistics. Use fixed + # threshold. + thresh = 500 + else: + mean_nz_diffs = np.mean(nz_diffs) + std_nz_diffs = np.std(nz_diffs) + med_nz_diffs = np.median(nz_diffs) + mad_nz_diffs = np.median(np.abs(nz_diffs - med_nz_diffs)) + if mean_nz_diffs / float(med_nz_diffs) < 3: + # Relatively small variation, keep (almost) everything + thresh = mean_nz_diffs + 3*std_nz_diffs + else: + # Large variation, filter more agressively. Use median and + # median absolute deviation (MAD) as they are less sensitive to + # outliers. However, allow differences < 500 scanlines as they + # occur quite often. + thresh = max(500, med_nz_diffs + 3*mad_nz_diffs) + self.scans = self.scans[diffs <= thresh] + + LOG.debug('Removed %s scanline(s) with corrupt scanline numbers', + str(len(along_track) - len(self.scans))) + + results.update({'n_corr': self.scans['scan_line_number'], + 'within_range': within_range, + 'diffs': diffs, + 'thresh': thresh, + 'nz_diffs': nz_diffs}) + return results + + def correct_times_thresh(self, max_diff_from_t0_head=6*60*1000, + min_frac_near_t0_head=0.01, + max_diff_from_ideal_t=10*1000): + """Correct corrupted timestamps using a threshold approach. + + The threshold approach is based on the scanline number and the header + timestamp. It also works if the majority of scanlines has corrupted + timestamps. + + The header timestamp is used as a guideline to estimate the offset + between timestamps computed from the scanline number and the actual + scanline timestamps in the data. If header timestamp and scanline + timestamps do not match, no correction is applied. + + Once the offset has been estimated, one can calculate the ideal + timestamps based on the scanline number. Timestamps deviating more than + a certain threshold from the ideal timestamps are replaced by + the ideal timestamps. + + Example files having corrupt timestamps: + - NSS.GHRR.NA.D81193.S2329.E0116.B1061214.WI + - NSS.GHRR.NL.D01035.S2342.E0135.B0192627.WI + + Args: + max_diff_from_t0_head (int): Threshold for offset estimation: A + scanline timestamp matches the header timestamp t0_head if it is + within the interval + + [t0_head - max_diff_from_t0_head, + t0_head + max_diff_from_t0_head] + + around the header timestamp. + min_frac_near_t0_head (float): Specifies the minimum fraction of + scanline timestamps matching the header timestamp required for + applying the correction. + max_diff_from_ideal_t (float): Threshold for timestamp correction: + If a scanline timestamp deviates more than max_diff_from_ideal_t + from the ideal timestamp, it is regarded as corrupt and will be + replaced with the ideal timestamp. + Returns: + Intermediate and final results (for plotting purpose) + + """ + results = {} + dt64_msec = ">M8[ms]" + + # Check whether scanline number increases monotonically + nums = self.scans["scan_line_number"] + results.update({'t': self.utcs.copy(), 'n': nums}) + if np.any(np.diff(nums) < 0): + LOG.error("Cannot perform timestamp correction. Scanline number " + "does not increase monotonically.") + results['fail_reason'] = "Scanline number jumps backwards" + return results + + # Convert time to milliseconds since 1970-01-01 + t = self.utcs.astype("i8") + try: + t0_head = np.array([self.get_header_timestamp().isoformat()], + dtype="datetime64[ms]").astype("i8")[0] + except ValueError as err: + LOG.error("Cannot perform timestamp correction: %s", err) + return + + # Compute ideal timestamps based on the scanline number. Still + # without offset, i.e. scanline 0 has timestamp 1970-01-01 00:00 + tn = self.lineno2msec(nums) + + # Try to determine the timestamp t0 of the first scanline. Since none + # of the actual timestamps is trustworthy, use the header timestamp + # as a guideline. However, the header timestamp may also be corrupted, + # so we only apply corrections if there is a minimum fraction of + # scanlines whose timestamps match the header timestamp. + # + # 1) Compute offsets between actual timestamps and idealized timestamps + offsets = t - tn + + # 2) If the offsets of a certain minimum fraction of scanlines are + # within a certain interval around the header timestamp, estimate + # t0 by calculating the median offset among these timestamps. If not, + # we do not have reliable information and cannot proceed. + near_t0_head = np.where( + np.fabs(offsets - t0_head) <= max_diff_from_t0_head)[0] + results.update({'offsets': offsets, + 't0_head': t0_head, + 'max_diff_from_t0_head': max_diff_from_t0_head}) + if near_t0_head.size / float(nums.size) >= min_frac_near_t0_head: + t0 = np.median(offsets[near_t0_head]) + else: + LOG.error("Timestamp mismatch. Cannot perform correction.") + results['fail_reason'] = "Timestamp mismatch" + return results + + # Add estimated offset to the ideal timestamps + tn += t0 + + # Replace timestamps deviating more than a certain threshold from the + # ideal timestamp with the ideal timestamp. + corrupt_lines = np.where(np.fabs(t - tn) > max_diff_from_ideal_t) + self.utcs[corrupt_lines] = tn[corrupt_lines].astype(dt64_msec) + LOG.debug("Corrected %s timestamp(s)", str(len(corrupt_lines[0]))) + + results.update({'tn': tn, 'tcorr': self.utcs, 't0': t0}) + return results + + @abstractproperty + def tsm_affected_intervals(self): + """Specify time intervals being affected by the scan motor problem. + + Returns: + dict: Affected time intervals. A dictionary containing a list of + (start, end) tuples for each affected platform. Both start and + end must be datetime.datetime objects. + + """ + raise NotImplementedError + + def is_tsm_affected(self): + """Determine whether this orbit is affected by the scan motor problem. + + Returns: + bool: True if the orbit is affected, False otherwise. + + """ + self.get_times() + ts = self.times[0] + te = self.times[-1] + try: + for interval in self.tsm_affected_intervals[self.spacecraft_id]: + if ts >= interval[0] and te <= interval[1]: + # Found a matching interval + return True + + # No matching interval, orbit is not affected + return False + except KeyError: + # Platform is not affected at all + return False + + def get_midnight_scanline(self): + """Find the scanline where the UTC date increases by one day. + + Returns: + int: The midnight scanline if it exists and is unique. + None, else. + + """ + self.get_times() + d0 = np.datetime64(datetime.date(1970, 1, 1), 'D') + days = (self.utcs.astype('datetime64[D]') - d0).astype(int) + incr = np.where(np.diff(days) == 1)[0] + if len(incr) != 1: + if len(incr) > 1: + LOG.warning('Unable to determine midnight scanline: ' + 'UTC date increases more than once. ') + return None + else: + return incr[0] + + def get_miss_lines(self): + """Find missing scanlines. + + I.e. scanlines which were dropped for some reason or were never recorded. + + Returns: + Indices of missing scanlines + + """ + # Compare scanline number against the ideal case (1, 2, 3, ...) and + # find the missing line numbers. + ideal = set(range(1, self.scans['scan_line_number'][-1] + 1)) + missing = sorted(ideal.difference(set(self.scans['scan_line_number']))) + return np.array(missing, dtype=int) + + def mask_tsm_pixels(self, channels): + """Mask pixels affected by the scan motor issue.""" + idx = self.get_tsm_pixels(channels) + channels[idx] = np.nan + + @abstractmethod + def get_tsm_pixels(self, channels): + """Determine pixels affected by the scan motor issue. + + Channel selection is POD/KLM specific. + """ + raise NotImplementedError + + +def inherit_doc(cls): + """Make a class method inherit its docstring from the parent class. + + Copied from http://stackoverflow.com/a/8101598/5703449 . + """ + for name, func in vars(cls).items(): + if isinstance(func, types.FunctionType) and not func.__doc__: + for parent in cls.__bases__: + parfunc = getattr(parent, name, None) + if parfunc and getattr(parfunc, '__doc__', None): + func.__doc__ = parfunc.__doc__ + break + return cls diff --git a/pygac/tests/test_pod.py b/pygac/tests/test_pod.py index 3be438c3..a39e629d 100644 --- a/pygac/tests/test_pod.py +++ b/pygac/tests/test_pod.py @@ -17,6 +17,7 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . +"""Test module for the pod reading.""" import datetime as dt import unittest @@ -32,15 +33,16 @@ class TestPOD(unittest.TestCase): - """Test the POD GAC reader""" + """Test the POD GAC reader.""" longMessage = True def setUp(self): + """Set up the test.""" self.reader = GACPODReader() def test_decode_timestamps(self): - """Test POD timestamp decoding""" + """Test POD timestamp decoding.""" # Reference timestamps, one before 2000 one after 2000 t2000_ref = (2001, 335, 53644260) t1900_ref = (1983, 336, 35058207) @@ -67,6 +69,7 @@ def test_get_header_timestamp(self, decode_timestamps): @mock.patch('pygac.gac_pod.GACPODReader.decode_timestamps') def test_get_times(self, decode_timestamps): + """Test getting times.""" self.reader.scans = {'time_code': 123} self.reader._get_times() decode_timestamps.assert_called_with(123) @@ -86,7 +89,7 @@ def test_get_lonlat(self): numpy.testing.assert_array_equal(lons, lons_exp) numpy.testing.assert_array_equal(lats, lats_exp) - @mock.patch('pygac.gac_pod.get_tsm_idx') + @mock.patch('pygac.pod_reader.get_tsm_idx') def test_get_tsm_pixels(self, get_tsm_idx): """Test channel set used for TSM correction.""" ones = np.ones((409, 100)) @@ -104,7 +107,7 @@ def test_get_tsm_pixels(self, get_tsm_idx): def suite(): - """The suite for test_pod""" + """Test suite for test_pod.""" loader = unittest.TestLoader() mysuite = unittest.TestSuite() mysuite.addTest(loader.loadTestsFromTestCase(TestPOD)) diff --git a/pygac/tests/test_reader.py b/pygac/tests/test_reader.py index f1db636f..b430f7b6 100644 --- a/pygac/tests/test_reader.py +++ b/pygac/tests/test_reader.py @@ -62,7 +62,7 @@ def test_lineno2msec(self): @mock.patch('pygac.gac_reader.GACReader._get_lonlat') @mock.patch('pygac.gac_reader.GACReader._get_corrupt_mask') @mock.patch('pygac.gac_reader.GACReader._adjust_clock_drift') - @mock.patch('pygac.gac_reader.gtp.Gac_Lat_Lon_Interpolator') + @mock.patch('pygac.gac_reader.gtp.gac_lat_lon_interpolator') def test_get_lonlat(self, interpolator, adjust_clockdrift, get_corrupt_mask, get_lonlat): """Test common lon/lat computation.""" @@ -253,8 +253,8 @@ def test_get_angles(self, get_corrupt_mask): np.testing.assert_allclose(sun_zenith, expected_sun_zenith, atol=0.01) np.testing.assert_allclose(rel_azi, expected_rel_azi, atol=0.01) - @mock.patch('pygac.gac_reader.ConfigParser.ConfigParser.read') - @mock.patch('pygac.gac_reader.ConfigParser.ConfigParser.items') + @mock.patch('pygac.reader.ConfigParser.ConfigParser.read') + @mock.patch('pygac.reader.ConfigParser.ConfigParser.items') def test_get_tle_file(self, items, *mocks): """Test get_tle_file.""" # Use TLE name/dir from config file From fef8f075d7d8b1f7acfeec92ba80bff32ffd60ea Mon Sep 17 00:00:00 2001 From: stickler-ci Date: Tue, 3 Dec 2019 17:54:38 +0000 Subject: [PATCH 02/12] Fixing style errors. --- pygac/pod_reader.py | 9 ++++++--- pygac/reader.py | 11 +++++++---- 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/pygac/pod_reader.py b/pygac/pod_reader.py index ed597f10..37e03a4c 100644 --- a/pygac/pod_reader.py +++ b/pygac/pod_reader.py @@ -191,7 +191,8 @@ def correct_scan_line_numbers(self): super(PODReader, self).correct_scan_line_numbers() # cleaning up the data - min_scanline_number = np.amin(np.absolute(self.scans["scan_line_number"][:])) + min_scanline_number = np.amin( + np.absolute(self.scans["scan_line_number"][:])) if self.scans["scan_line_number"][0] == self.scans["scan_line_number"][-1] + 1: while self.scans["scan_line_number"][0] != min_scanline_number: self.scans = np.roll(self.scans, -1) @@ -231,7 +232,8 @@ def read(self, filename): head = np.fromfile(fd_, dtype=header0, count=1)[0] year, jday, _ = self.decode_timestamps(head["start_time"]) - start_date = (datetime.date(year, 1, 1) + datetime.timedelta(days=int(jday) - 1)) + start_date = (datetime.date(year, 1, 1) + + datetime.timedelta(days=int(jday) - 1)) if start_date < datetime.date(1992, 9, 8): header = header1 @@ -363,7 +365,8 @@ def _adjust_clock_drift(self): utcs=missed_utcs, clock_drift_adjust=True) except IndexError as err: - LOG.warning('Cannot perform clock drift correction: %s', str(err)) + LOG.warning( + 'Cannot perform clock drift correction: %s', str(err)) return complete_lons[missed - min_idx] = mlons diff --git a/pygac/reader.py b/pygac/reader.py index 4083fe29..ba8e6707 100644 --- a/pygac/reader.py +++ b/pygac/reader.py @@ -218,7 +218,7 @@ def to_datetime64(year, jday, msec): """ return (((year - 1970).astype('datetime64[Y]') - + (jday - 1).astype('timedelta64[D]')).astype('datetime64[ms]') + + (jday - 1).astype('timedelta64[D]')).astype('datetime64[ms]') + msec.astype('timedelta64[ms]')) @staticmethod @@ -299,7 +299,8 @@ def compute_lonlat(self, width, utcs=None, clock_drift_adjust=True): from pyorbital.geoloc_instrument_definitions import avhrr_gac from pyorbital.geoloc import compute_pixels, get_lonlatalt # TODO: Are we sure all satellites have this scan width in degrees ? - sgeom = avhrr_gac(utcs.astype(datetime.datetime), self.scan_points, 55.385) + sgeom = avhrr_gac(utcs.astype(datetime.datetime), + self.scan_points, 55.385) s_times = sgeom.times(t) tle1, tle2 = self.get_tle_lines() @@ -425,7 +426,8 @@ def tle2datetime64(times): doys = (times % 1000).astype('int') - 1 years = (times // 1000).astype('int') msecs = np.rint(24 * 3600 * 1000 * (times % 1)) - times64 = (years - 1970).astype('datetime64[Y]').astype('datetime64[ms]') + times64 = ( + years - 1970).astype('datetime64[Y]').astype('datetime64[ms]') times64 += doys.astype('timedelta64[D]') times64 += msecs.astype('timedelta64[ms]') @@ -580,7 +582,8 @@ def correct_times_median(self, year, jday, msec): # offset, e.g. the first scanline has timestamp 1970-01-01 00:00 msec_lineno = self.lineno2msec(self.scans["scan_line_number"]) - jday = np.where(np.logical_or(jday < 1, jday > 366), np.median(jday), jday) + jday = np.where(np.logical_or(jday < 1, jday > 366), + np.median(jday), jday) if_wrong_jday = np.ediff1d(jday, to_begin=0) jday = np.where(if_wrong_jday < 0, max(jday), jday) From 6e375cb8d15b3effbd63df46733cc20f727d63da Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Tue, 3 Dec 2019 18:58:25 +0100 Subject: [PATCH 03/12] Fix stickler line length --- .stickler.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.stickler.yml b/.stickler.yml index c266a863..4a6f4e14 100644 --- a/.stickler.yml +++ b/.stickler.yml @@ -1,5 +1,6 @@ linters: flake8: fixer: true + max-line-length: 120 fixers: enable: true From 068bb4b8da691f4bfec4dc224e72c94fae557085 Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Tue, 3 Dec 2019 20:20:46 +0100 Subject: [PATCH 04/12] Fix test_reader interpolator tests --- pygac/tests/test_reader.py | 29 +++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/pygac/tests/test_reader.py b/pygac/tests/test_reader.py index b430f7b6..6cf748e3 100644 --- a/pygac/tests/test_reader.py +++ b/pygac/tests/test_reader.py @@ -37,8 +37,10 @@ class TestGacReader(unittest.TestCase): @mock.patch.multiple('pygac.gac_reader.GACReader', __abstractmethods__=set()) - def setUp(self, *mocks): + @mock.patch('pygac.gac_reader.gtp.gac_lat_lon_interpolator') + def setUp(self, interpolator, *mocks): """Set up the tests.""" + self.interpolator = interpolator self.reader = GACReader() def test_to_datetime64(self): @@ -62,14 +64,13 @@ def test_lineno2msec(self): @mock.patch('pygac.gac_reader.GACReader._get_lonlat') @mock.patch('pygac.gac_reader.GACReader._get_corrupt_mask') @mock.patch('pygac.gac_reader.GACReader._adjust_clock_drift') - @mock.patch('pygac.gac_reader.gtp.gac_lat_lon_interpolator') - def test_get_lonlat(self, interpolator, adjust_clockdrift, + def test_get_lonlat(self, adjust_clockdrift, get_corrupt_mask, get_lonlat): """Test common lon/lat computation.""" lon_i = np.array([np.nan, 1, 2, 3, -180.1, 180.1]) lat_i = np.array([1, 2, 3, np.nan, -90.1, 90.1]) get_lonlat.return_value = lon_i, lat_i - interpolator.return_value = lon_i, lat_i + self.interpolator.return_value = lon_i, lat_i get_corrupt_mask.return_value = np.array( [0, 0, 1, 0, 0, 0], dtype=bool) @@ -84,7 +85,7 @@ def test_get_lonlat(self, interpolator, adjust_clockdrift, numpy.testing.assert_array_equal(lats, lats_exp) # Interpolation disabled - interpolator.reset_mock() + self.interpolator.reset_mock() adjust_clockdrift.reset_mock() self.reader.interpolate_coords = False self.reader.adjust_clock_drift = True @@ -92,11 +93,11 @@ def test_get_lonlat(self, interpolator, adjust_clockdrift, self.reader.get_lonlat() numpy.testing.assert_array_equal(lons, lons_exp) numpy.testing.assert_array_equal(lats, lats_exp) - interpolator.assert_not_called() + self.interpolator.assert_not_called() adjust_clockdrift.assert_called() # Clock drift adjustment disabled - interpolator.reset_mock() + self.interpolator.reset_mock() adjust_clockdrift.reset_mock() self.reader.interpolate_coords = True self.reader.adjust_clock_drift = False @@ -104,11 +105,11 @@ def test_get_lonlat(self, interpolator, adjust_clockdrift, self.reader.get_lonlat() numpy.testing.assert_array_equal(lons, lons_exp) numpy.testing.assert_array_equal(lats, lats_exp) - interpolator.assert_called() + self.interpolator.assert_called() adjust_clockdrift.assert_not_called() # Test caching - methods = [get_lonlat, interpolator, + methods = [get_lonlat, self.interpolator, adjust_clockdrift, get_corrupt_mask] for method in methods: method.reset_mock() @@ -124,12 +125,16 @@ def test_interpolate(self, _get_lonlat, _adjust_clock_drift, """Test interpolate method in get_lonlat.""" self.lons = None self.lats = None - lons = 90 * np.random.rand(17, 51) - lats = 90 * np.random.rand(17, 51) - _get_lonlat.return_value = lons, lats + lr_lons = 90 * np.random.rand(17, 51) + lr_lats = 90 * np.random.rand(17, 51) + _get_lonlat.return_value = lr_lons, lr_lats self.interpolate_coors = True + self.interpolator.reset_mock() + self.interpolator.return_value = (90 * np.random.rand(17, 409), + 90 * np.random.rand(17, 409)) lons, lats = self.reader.get_lonlat() self.assertEqual(lons.shape[1], 409) + self.interpolator.assert_called_once_with(lr_lons, lr_lats) @mock.patch('pygac.gac_reader.GACReader._get_corrupt_mask') def test_get_corrupt_mask(self, get_corrupt_mask): From 537ab85c2902779f3e0cca08247ab945b3b4447f Mon Sep 17 00:00:00 2001 From: Stephan Finkensieper Date: Wed, 4 Dec 2019 09:41:44 +0000 Subject: [PATCH 05/12] Fix main methods --- pygac/gac_klm.py | 36 ++++++++++++++++++++++++++++++++++++ pygac/gac_pod.py | 36 ++++++++++++++++++++++++++++++++++++ pygac/lac_klm.py | 10 ++++++---- pygac/lac_pod.py | 13 ++++++++----- pygac/pod_reader.py | 35 ----------------------------------- 5 files changed, 86 insertions(+), 44 deletions(-) diff --git a/pygac/gac_klm.py b/pygac/gac_klm.py index 47b49e9b..9ce9ca43 100644 --- a/pygac/gac_klm.py +++ b/pygac/gac_klm.py @@ -32,6 +32,7 @@ from __future__ import print_function +import datetime import logging import numpy as np @@ -195,3 +196,38 @@ def __init__(self, *args, **kwargs): GACReader.__init__(self, *args, **kwargs) self.scanline_type = scanline self.offset = 4608 + + +def main(filename, start_line, end_line): + """Generate a l1c file.""" + from pygac import gac_io + tic = datetime.datetime.now() + reader = GACKLMReader() + reader.read(filename) + reader.get_lonlat() + channels = reader.get_calibrated_channels() + sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles() + qual_flags = reader.get_qual_flags() + if (np.all(reader.mask)): + print("ERROR: All data is masked out. Stop processing") + raise ValueError("All data is masked out.") + + gac_io.save_gac(reader.spacecraft_name, + reader.utcs, + reader.lats, reader.lons, + channels[:, :, 0], channels[:, :, 1], + channels[:, :, 2], + channels[:, :, 3], + channels[:, :, 4], + channels[:, :, 5], + sun_zen, sat_zen, sun_azi, sat_azi, rel_azi, + qual_flags, start_line, end_line, + reader.filename, + reader.get_midnight_scanline(), + reader.get_miss_lines()) + LOG.info("pygac took: %s", str(datetime.datetime.now() - tic)) + + +if __name__ == "__main__": + import sys + main(sys.argv[1], sys.argv[2], sys.argv[3]) diff --git a/pygac/gac_pod.py b/pygac/gac_pod.py index 3c344781..e76607e6 100644 --- a/pygac/gac_pod.py +++ b/pygac/gac_pod.py @@ -31,6 +31,7 @@ http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/c3/sec3-1.htm """ +import datetime import logging import numpy as np @@ -64,3 +65,38 @@ def __init__(self, *args, **kwargs): self.scanline_type = scanline self.offset = 3220 self.scan_points = np.arange(3.5, 2048, 5) + + +def main(filename, start_line, end_line): + """Generate a l1c file.""" + from pygac import gac_io + tic = datetime.datetime.now() + reader = GACPODReader() + reader.read(filename) + reader.get_lonlat() + channels = reader.get_calibrated_channels() + sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles() + + qual_flags = reader.get_qual_flags() + if (np.all(reader.mask)): + print("ERROR: All data is masked out. Stop processing") + raise ValueError("All data is masked out.") + gac_io.save_gac(reader.spacecraft_name, + reader.utcs, + reader.lats, reader.lons, + channels[:, :, 0], channels[:, :, 1], + np.full_like(channels[:, :, 0], np.nan), + channels[:, :, 2], + channels[:, :, 3], + channels[:, :, 4], + sun_zen, sat_zen, sun_azi, sat_azi, rel_azi, + qual_flags, start_line, end_line, + reader.filename, + reader.get_midnight_scanline(), + reader.get_miss_lines()) + LOG.info("pygac took: %s", str(datetime.datetime.now() - tic)) + + +if __name__ == "__main__": + import sys + main(sys.argv[1], sys.argv[2], sys.argv[3]) diff --git a/pygac/lac_klm.py b/pygac/lac_klm.py index 56727a0e..3daab4f3 100644 --- a/pygac/lac_klm.py +++ b/pygac/lac_klm.py @@ -203,9 +203,8 @@ def main(filename, start_line, end_line): reader.get_lonlat() channels = reader.get_calibrated_channels() sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles() - - mask, qual_flags = reader.get_corrupt_mask() - if (np.all(mask)): + qual_flags = reader.get_qual_flags() + if (np.all(reader.mask)): print("ERROR: All data is masked out. Stop processing") raise ValueError("All data is masked out.") @@ -218,7 +217,10 @@ def main(filename, start_line, end_line): channels[:, :, 4], channels[:, :, 5], sun_zen, sat_zen, sun_azi, sat_azi, rel_azi, - mask, qual_flags, start_line, end_line, reader.get_ch3_switch()) + qual_flags, start_line, end_line, + reader.filename, + reader.get_midnight_scanline(), + reader.get_miss_lines()) LOG.info("pygac took: %s", str(datetime.datetime.now() - tic)) diff --git a/pygac/lac_pod.py b/pygac/lac_pod.py index 88f6ba32..73a32d8b 100644 --- a/pygac/lac_pod.py +++ b/pygac/lac_pod.py @@ -63,28 +63,31 @@ def __init__(self, *args, **kwargs): def main(filename, start_line, end_line): """Generate a l1c file.""" + from pygac import gac_io tic = datetime.datetime.now() reader = LACPODReader() reader.read(filename) reader.get_lonlat() - reader.adjust_clock_drift() channels = reader.get_calibrated_channels() sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles() - mask, qual_flags = reader.get_corrupt_mask() - if (np.all(mask)): + qual_flags = reader.get_qual_flags() + if (np.all(reader.mask)): print("ERROR: All data is masked out. Stop processing") raise ValueError("All data is masked out.") gac_io.save_gac(reader.spacecraft_name, reader.utcs, reader.lats, reader.lons, channels[:, :, 0], channels[:, :, 1], - np.ones_like(channels[:, :, 0]) * -1, + np.full_like(channels[:, :, 0], np.nan), channels[:, :, 2], channels[:, :, 3], channels[:, :, 4], sun_zen, sat_zen, sun_azi, sat_azi, rel_azi, - mask, qual_flags, start_line, end_line) + qual_flags, start_line, end_line, + reader.filename, + reader.get_midnight_scanline(), + reader.get_miss_lines()) LOG.info("pygac took: %s", str(datetime.datetime.now() - tic)) diff --git a/pygac/pod_reader.py b/pygac/pod_reader.py index 37e03a4c..b4f8d71e 100644 --- a/pygac/pod_reader.py +++ b/pygac/pod_reader.py @@ -457,38 +457,3 @@ def get_tsm_pixels(self, channels): """ return get_tsm_idx(channels[:, :, 0], channels[:, :, 1], channels[:, :, 3], channels[:, :, 4]) - - -def main(filename, start_line, end_line): - """Generate a l1c file.""" - from pygac import gac_io - tic = datetime.datetime.now() - reader = PODReader() - reader.read(filename) - reader.get_lonlat() - channels = reader.get_calibrated_channels() - sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles() - qual_flags = reader.get_qual_flags() - if np.all(reader.mask): - print("ERROR: All data is masked out. Stop processing") - raise ValueError("All data is masked out.") - - gac_io.save_gac(reader.spacecraft_name, - reader.utcs, - reader.lats, reader.lons, - channels[:, :, 0], channels[:, :, 1], - np.full_like(channels[:, :, 0], np.nan), - channels[:, :, 2], - channels[:, :, 3], - channels[:, :, 4], - sun_zen, sat_zen, sun_azi, sat_azi, rel_azi, - qual_flags, start_line, end_line, - reader.filename, - reader.get_midnight_scanline(), - reader.get_miss_lines()) - LOG.info("pygac took: %s", str(datetime.datetime.now() - tic)) - - -if __name__ == "__main__": - import sys - main(sys.argv[1], sys.argv[2], sys.argv[3]) From a36a388dd12587482877caf213d8fc05cf47da02 Mon Sep 17 00:00:00 2001 From: Stephan Finkensieper Date: Wed, 4 Dec 2019 09:47:42 +0000 Subject: [PATCH 06/12] Factorize main methods --- pygac/gac_klm.py | 30 ++---------------------------- pygac/gac_pod.py | 30 ++---------------------------- pygac/klm_reader.py | 30 ++++++++++++++++++++++++++++++ pygac/lac_klm.py | 28 +--------------------------- pygac/lac_pod.py | 30 ++---------------------------- pygac/pod_reader.py | 30 ++++++++++++++++++++++++++++++ 6 files changed, 67 insertions(+), 111 deletions(-) diff --git a/pygac/gac_klm.py b/pygac/gac_klm.py index 9ce9ca43..c7913e89 100644 --- a/pygac/gac_klm.py +++ b/pygac/gac_klm.py @@ -38,7 +38,7 @@ import numpy as np from pygac.gac_reader import GACReader -from pygac.klm_reader import KLMReader +from pygac.klm_reader import KLMReader, main_klm LOG = logging.getLogger(__name__) @@ -199,33 +199,7 @@ def __init__(self, *args, **kwargs): def main(filename, start_line, end_line): - """Generate a l1c file.""" - from pygac import gac_io - tic = datetime.datetime.now() - reader = GACKLMReader() - reader.read(filename) - reader.get_lonlat() - channels = reader.get_calibrated_channels() - sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles() - qual_flags = reader.get_qual_flags() - if (np.all(reader.mask)): - print("ERROR: All data is masked out. Stop processing") - raise ValueError("All data is masked out.") - - gac_io.save_gac(reader.spacecraft_name, - reader.utcs, - reader.lats, reader.lons, - channels[:, :, 0], channels[:, :, 1], - channels[:, :, 2], - channels[:, :, 3], - channels[:, :, 4], - channels[:, :, 5], - sun_zen, sat_zen, sun_azi, sat_azi, rel_azi, - qual_flags, start_line, end_line, - reader.filename, - reader.get_midnight_scanline(), - reader.get_miss_lines()) - LOG.info("pygac took: %s", str(datetime.datetime.now() - tic)) + return main_klm(GACKLMReader, filename, start_line, end_line) if __name__ == "__main__": diff --git a/pygac/gac_pod.py b/pygac/gac_pod.py index e76607e6..7bfcc903 100644 --- a/pygac/gac_pod.py +++ b/pygac/gac_pod.py @@ -36,7 +36,7 @@ import numpy as np -from pygac.pod_reader import PODReader +from pygac.pod_reader import PODReader, main_pod from pygac.gac_reader import GACReader LOG = logging.getLogger(__name__) @@ -68,33 +68,7 @@ def __init__(self, *args, **kwargs): def main(filename, start_line, end_line): - """Generate a l1c file.""" - from pygac import gac_io - tic = datetime.datetime.now() - reader = GACPODReader() - reader.read(filename) - reader.get_lonlat() - channels = reader.get_calibrated_channels() - sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles() - - qual_flags = reader.get_qual_flags() - if (np.all(reader.mask)): - print("ERROR: All data is masked out. Stop processing") - raise ValueError("All data is masked out.") - gac_io.save_gac(reader.spacecraft_name, - reader.utcs, - reader.lats, reader.lons, - channels[:, :, 0], channels[:, :, 1], - np.full_like(channels[:, :, 0], np.nan), - channels[:, :, 2], - channels[:, :, 3], - channels[:, :, 4], - sun_zen, sat_zen, sun_azi, sat_azi, rel_azi, - qual_flags, start_line, end_line, - reader.filename, - reader.get_midnight_scanline(), - reader.get_miss_lines()) - LOG.info("pygac took: %s", str(datetime.datetime.now() - tic)) + return main_pod(GACPODReader, filename, start_line, end_line) if __name__ == "__main__": diff --git a/pygac/klm_reader.py b/pygac/klm_reader.py index 247d1429..c8c173f3 100644 --- a/pygac/klm_reader.py +++ b/pygac/klm_reader.py @@ -730,3 +730,33 @@ def get_tsm_pixels(self, channels): """ return get_tsm_idx(channels[:, :, 0], channels[:, :, 1], channels[:, :, 4], channels[:, :, 5]) + + +def main_klm(reader_cls, filename, start_line, end_line): + """Generate a l1c file.""" + from pygac import gac_io + tic = datetime.datetime.now() + reader = reader_cls() + reader.read(filename) + reader.get_lonlat() + channels = reader.get_calibrated_channels() + sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles() + qual_flags = reader.get_qual_flags() + if (np.all(reader.mask)): + print("ERROR: All data is masked out. Stop processing") + raise ValueError("All data is masked out.") + + gac_io.save_gac(reader.spacecraft_name, + reader.utcs, + reader.lats, reader.lons, + channels[:, :, 0], channels[:, :, 1], + channels[:, :, 2], + channels[:, :, 3], + channels[:, :, 4], + channels[:, :, 5], + sun_zen, sat_zen, sun_azi, sat_azi, rel_azi, + qual_flags, start_line, end_line, + reader.filename, + reader.get_midnight_scanline(), + reader.get_miss_lines()) + LOG.info("pygac took: %s", str(datetime.datetime.now() - tic)) diff --git a/pygac/lac_klm.py b/pygac/lac_klm.py index 3daab4f3..c06b52dd 100644 --- a/pygac/lac_klm.py +++ b/pygac/lac_klm.py @@ -195,33 +195,7 @@ def __init__(self, *args, **kwargs): def main(filename, start_line, end_line): - """Generate a l1c file.""" - from pygac import gac_io - tic = datetime.datetime.now() - reader = LACKLMReader() - reader.read(filename) - reader.get_lonlat() - channels = reader.get_calibrated_channels() - sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles() - qual_flags = reader.get_qual_flags() - if (np.all(reader.mask)): - print("ERROR: All data is masked out. Stop processing") - raise ValueError("All data is masked out.") - - gac_io.save_gac(reader.spacecraft_name, - reader.utcs, - reader.lats, reader.lons, - channels[:, :, 0], channels[:, :, 1], - channels[:, :, 2], - channels[:, :, 3], - channels[:, :, 4], - channels[:, :, 5], - sun_zen, sat_zen, sun_azi, sat_azi, rel_azi, - qual_flags, start_line, end_line, - reader.filename, - reader.get_midnight_scanline(), - reader.get_miss_lines()) - LOG.info("pygac took: %s", str(datetime.datetime.now() - tic)) + return main_klm(LACKLMReader, filename, start_line, end_line) if __name__ == "__main__": diff --git a/pygac/lac_pod.py b/pygac/lac_pod.py index 73a32d8b..f28b3cb3 100644 --- a/pygac/lac_pod.py +++ b/pygac/lac_pod.py @@ -30,7 +30,7 @@ import numpy as np from pygac import gac_io -from pygac.pod_reader import PODReader +from pygac.pod_reader import PODReader, main_pod from pygac.lac_reader import LACReader LOG = logging.getLogger(__name__) @@ -62,33 +62,7 @@ def __init__(self, *args, **kwargs): def main(filename, start_line, end_line): - """Generate a l1c file.""" - from pygac import gac_io - tic = datetime.datetime.now() - reader = LACPODReader() - reader.read(filename) - reader.get_lonlat() - channels = reader.get_calibrated_channels() - sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles() - - qual_flags = reader.get_qual_flags() - if (np.all(reader.mask)): - print("ERROR: All data is masked out. Stop processing") - raise ValueError("All data is masked out.") - gac_io.save_gac(reader.spacecraft_name, - reader.utcs, - reader.lats, reader.lons, - channels[:, :, 0], channels[:, :, 1], - np.full_like(channels[:, :, 0], np.nan), - channels[:, :, 2], - channels[:, :, 3], - channels[:, :, 4], - sun_zen, sat_zen, sun_azi, sat_azi, rel_azi, - qual_flags, start_line, end_line, - reader.filename, - reader.get_midnight_scanline(), - reader.get_miss_lines()) - LOG.info("pygac took: %s", str(datetime.datetime.now() - tic)) + return main_pod(LACPODReader, filename, start_line, end_line) if __name__ == "__main__": diff --git a/pygac/pod_reader.py b/pygac/pod_reader.py index b4f8d71e..94032837 100644 --- a/pygac/pod_reader.py +++ b/pygac/pod_reader.py @@ -457,3 +457,33 @@ def get_tsm_pixels(self, channels): """ return get_tsm_idx(channels[:, :, 0], channels[:, :, 1], channels[:, :, 3], channels[:, :, 4]) + + +def main_pod(reader_cls, filename, start_line, end_line): + """Generate a l1c file.""" + from pygac import gac_io + tic = datetime.datetime.now() + reader = reader_cls() + reader.read(filename) + reader.get_lonlat() + channels = reader.get_calibrated_channels() + sat_azi, sat_zen, sun_azi, sun_zen, rel_azi = reader.get_angles() + + qual_flags = reader.get_qual_flags() + if (np.all(reader.mask)): + print("ERROR: All data is masked out. Stop processing") + raise ValueError("All data is masked out.") + gac_io.save_gac(reader.spacecraft_name, + reader.utcs, + reader.lats, reader.lons, + channels[:, :, 0], channels[:, :, 1], + np.full_like(channels[:, :, 0], np.nan), + channels[:, :, 2], + channels[:, :, 3], + channels[:, :, 4], + sun_zen, sat_zen, sun_azi, sat_azi, rel_azi, + qual_flags, start_line, end_line, + reader.filename, + reader.get_midnight_scanline(), + reader.get_miss_lines()) + LOG.info("pygac took: %s", str(datetime.datetime.now() - tic)) From 47223bc80f5c05498830803ac982bc40529087c2 Mon Sep 17 00:00:00 2001 From: Stephan Finkensieper Date: Wed, 4 Dec 2019 10:58:19 +0000 Subject: [PATCH 07/12] Make stickler happy --- pygac/klm_reader.py | 1 + pygac/lac_klm.py | 3 +-- pygac/lac_pod.py | 2 -- 3 files changed, 2 insertions(+), 4 deletions(-) diff --git a/pygac/klm_reader.py b/pygac/klm_reader.py index c8c173f3..bfa8e05e 100644 --- a/pygac/klm_reader.py +++ b/pygac/klm_reader.py @@ -30,6 +30,7 @@ http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c8/sec83142-1.htm """ +import datetime import logging import numpy as np diff --git a/pygac/lac_klm.py b/pygac/lac_klm.py index c06b52dd..43087a4d 100644 --- a/pygac/lac_klm.py +++ b/pygac/lac_klm.py @@ -24,12 +24,11 @@ """Reader for LAC KLM data.""" -import datetime import logging import numpy as np -from pygac.klm_reader import KLMReader +from pygac.klm_reader import KLMReader, main_klm from pygac.lac_reader import LACReader LOG = logging.getLogger(__name__) diff --git a/pygac/lac_pod.py b/pygac/lac_pod.py index f28b3cb3..59f478d5 100644 --- a/pygac/lac_pod.py +++ b/pygac/lac_pod.py @@ -24,12 +24,10 @@ """The LAC POD reader routines.""" -import datetime import logging import numpy as np -from pygac import gac_io from pygac.pod_reader import PODReader, main_pod from pygac.lac_reader import LACReader From c6319e42e0457d98f109e48c93b8ee551de488e9 Mon Sep 17 00:00:00 2001 From: Stephan Finkensieper Date: Wed, 4 Dec 2019 12:53:47 +0000 Subject: [PATCH 08/12] Fix TSM test --- pygac/tests/test_klm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygac/tests/test_klm.py b/pygac/tests/test_klm.py index d9da7761..ed94b856 100644 --- a/pygac/tests/test_klm.py +++ b/pygac/tests/test_klm.py @@ -131,7 +131,7 @@ def test_postproc(self, get_ch3_switch): self.reader.postproc(channels) # masks in-place numpy.testing.assert_array_equal(channels, masked_exp) - @mock.patch('pygac.gac_klm.get_tsm_idx') + @mock.patch('pygac.klm_reader.get_tsm_idx') def test_get_tsm_pixels(self, get_tsm_idx): """Test channel set used for TSM correction.""" ones = np.ones((409, 100)) From 8a7f2bd6f6f83c771850f50f6267fa704ee3c444 Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Wed, 4 Dec 2019 15:01:22 +0100 Subject: [PATCH 09/12] Disable the static attitude correction More thorough comparison needs to be done on this --- pygac/reader.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/pygac/reader.py b/pygac/reader.py index e91a5bc3..2d4b565f 100644 --- a/pygac/reader.py +++ b/pygac/reader.py @@ -287,11 +287,13 @@ def compute_lonlat(self, width, utcs=None, clock_drift_adjust=True): self.head["constant_yaw_attitude_error"] / 1e3]) else: try: - rpy_spacecraft = rpy_coeffs[self.spacecraft_name] - rpy = [rpy_spacecraft['roll'], - rpy_spacecraft['pitch'], - rpy_spacecraft['yaw']] - LOG.debug("Using static attitude correction") + # This needs to be checked thoroughly first + # rpy_spacecraft = rpy_coeffs[self.spacecraft_name] + # rpy = [rpy_spacecraft['roll'], + # rpy_spacecraft['pitch'], + # rpy_spacecraft['yaw']] + # LOG.debug("Using static attitude correction") + raise KeyError except KeyError: LOG.debug("Not applying attitude correction") rpy = [0, 0, 0] From 3af409391e624641429e1c18e58b8ad4a7a145e6 Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Wed, 4 Dec 2019 15:38:20 +0100 Subject: [PATCH 10/12] Document the `scan_points` attribute --- pygac/gac_pod.py | 8 ++++++-- pygac/lac_pod.py | 7 ++++++- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/pygac/gac_pod.py b/pygac/gac_pod.py index 7bfcc903..740788e4 100644 --- a/pygac/gac_pod.py +++ b/pygac/gac_pod.py @@ -31,7 +31,6 @@ http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/c3/sec3-1.htm """ -import datetime import logging import numpy as np @@ -57,7 +56,11 @@ class GACPODReader(GACReader, PODReader): - """The GAC POD reader class.""" + """The GAC POD reader class. + + The `scan_points` attributes provides the position of the longitude and latitude points to + compute relative to the full swath width. + """ def __init__(self, *args, **kwargs): """Init the GAC POD reader.""" @@ -68,6 +71,7 @@ def __init__(self, *args, **kwargs): def main(filename, start_line, end_line): + """Generate a l1c file.""" return main_pod(GACPODReader, filename, start_line, end_line) diff --git a/pygac/lac_pod.py b/pygac/lac_pod.py index 59f478d5..a848b31a 100644 --- a/pygac/lac_pod.py +++ b/pygac/lac_pod.py @@ -49,7 +49,11 @@ class LACPODReader(LACReader, PODReader): - """The LAC POD reader.""" + """The LAC POD reader. + + The `scan_points` attributes provides the position of the longitude and latitude points to + compute relative to the full swath width. + """ def __init__(self, *args, **kwargs): """Init the LAC POD reader.""" @@ -60,6 +64,7 @@ def __init__(self, *args, **kwargs): def main(filename, start_line, end_line): + """Generate a l1c file.""" return main_pod(LACPODReader, filename, start_line, end_line) From b14827150ceac7643b47158e07ad038322f65a73 Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Wed, 4 Dec 2019 15:42:55 +0100 Subject: [PATCH 11/12] Add documentation for the `offset` attribute --- pygac/gac_klm.py | 7 +++++-- pygac/gac_pod.py | 2 ++ pygac/lac_klm.py | 6 +++++- pygac/lac_pod.py | 2 ++ 4 files changed, 14 insertions(+), 3 deletions(-) diff --git a/pygac/gac_klm.py b/pygac/gac_klm.py index c7913e89..293ac1fc 100644 --- a/pygac/gac_klm.py +++ b/pygac/gac_klm.py @@ -32,7 +32,6 @@ from __future__ import print_function -import datetime import logging import numpy as np @@ -189,7 +188,10 @@ class GACKLMReader(GACReader, KLMReader): - """The GAC KLM reader class.""" + """The GAC KLM reader class. + + The offset attribute tells where in the file the scanline data starts. + """ def __init__(self, *args, **kwargs): """Init the GAC KLM reader.""" @@ -199,6 +201,7 @@ def __init__(self, *args, **kwargs): def main(filename, start_line, end_line): + """Generate a l1c file.""" return main_klm(GACKLMReader, filename, start_line, end_line) diff --git a/pygac/gac_pod.py b/pygac/gac_pod.py index 740788e4..9cabb56f 100644 --- a/pygac/gac_pod.py +++ b/pygac/gac_pod.py @@ -60,6 +60,8 @@ class GACPODReader(GACReader, PODReader): The `scan_points` attributes provides the position of the longitude and latitude points to compute relative to the full swath width. + + The offset attribute tells where in the file the scanline data starts. """ def __init__(self, *args, **kwargs): diff --git a/pygac/lac_klm.py b/pygac/lac_klm.py index 43087a4d..f2a50d15 100644 --- a/pygac/lac_klm.py +++ b/pygac/lac_klm.py @@ -181,7 +181,10 @@ class LACKLMReader(LACReader, KLMReader): - """The LAC KLM reader.""" + """The LAC KLM reader. + + The offset attribute tells where in the file the scanline data starts. + """ def __init__(self, *args, **kwargs): """Init the LAC KLM reader.""" @@ -194,6 +197,7 @@ def __init__(self, *args, **kwargs): def main(filename, start_line, end_line): + """Generate a l1c file.""" return main_klm(LACKLMReader, filename, start_line, end_line) diff --git a/pygac/lac_pod.py b/pygac/lac_pod.py index a848b31a..760e09ba 100644 --- a/pygac/lac_pod.py +++ b/pygac/lac_pod.py @@ -53,6 +53,8 @@ class LACPODReader(LACReader, PODReader): The `scan_points` attributes provides the position of the longitude and latitude points to compute relative to the full swath width. + + The offset attribute tells where in the file the scanline data starts. """ def __init__(self, *args, **kwargs): From fe07e305751db03993ea8a62b260683e93bcadcd Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Thu, 5 Dec 2019 10:56:53 +0100 Subject: [PATCH 12/12] Fix pygac description to include LAC --- README.md | 10 ++++------ setup.py | 13 ++++++++----- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index 28a7bb4f..37edd75f 100644 --- a/README.md +++ b/README.md @@ -1,12 +1,10 @@ pygac ===== -[![Build Status](https://travis-ci.org/adybbroe/pygac.png?branch=feature-clock)](https://travis-ci.org/adybbroe/pygac) -[![Coverage Status](https://coveralls.io/repos/adybbroe/pygac/badge.png?branch=feature-clock)](https://coveralls.io/r/adybbroe/pygac?branch=develop) -[![Code Health](https://landscape.io/github/pytroll/pygac/develop/landscape.png)](https://landscape.io/github/pytroll/pygac/develop) +[![Build Status](https://travis-ci.org/adybbroe/pygac.png?branch=feature-clock)](https://travis-ci.org/pytroll/pygac) +[![Coverage Status](https://coveralls.io/repos/adybbroe/pygac/badge.png?branch=feature-clock)](https://coveralls.io/r/pytroll/pygac?branch=develop) -A python package to read, calibrate and navigate NOAA AVHRR GAC data. +A python package to read, calibrate and navigate NOAA and Metop AVHRR GAC and LAC data. + -Abhay Devatshale, Martin Raspaud and Adam Dybbroe -April 2014, Norrkoping, Sweden diff --git a/setup.py b/setup.py index cefd25e1..22bd50ce 100644 --- a/setup.py +++ b/setup.py @@ -19,7 +19,7 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . - +"""The setup module.""" try: with open("./README", "r") as fd: @@ -39,19 +39,22 @@ def set_builtin(name, value): + """Set builtin.""" if isinstance(__builtins__, dict): __builtins__[name] = value else: setattr(__builtins__, name, value) -class build_ext(_build_ext): +class build_ext(_build_ext): # noqa """Work around to bootstrap numpy includes in to extensions. + Copied from: http://stackoverflow.com/questions/19919905/how-to-bootstrap-numpy-installation-in-setup-py """ def finalize_options(self): + """Finalize options.""" _build_ext.finalize_options(self) # Prevent numpy from thinking it is still in its setup process: set_builtin('__NUMPY_SETUP__', False) @@ -77,9 +80,9 @@ def finalize_options(self): setup(name='pygac', version=version.__version__, - description='NOAA AVHRR GAC reader and calibration', - author='Abhay Devasthale', - author_email='adam.dybbroe@smhi.se', + description='NOAA AVHRR GAC/LAC reader and calibration', + author='Abhay Devasthale, Martin Raspaud', + author_email='martin.raspaud@smhi.se', classifiers=["Development Status :: 4 - Beta", "Intended Audience :: Science/Research", "License :: OSI Approved :: GNU General Public License v3 " +