diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 4071a67..ec1d934 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -4,7 +4,7 @@ fail_fast: false repos: - repo: https://github.com/astral-sh/ruff-pre-commit # Ruff version. - rev: 'v0.3.2' + rev: 'v0.6.2' hooks: - id: ruff args: [--fix] diff --git a/continuous_integration/environment.yaml b/continuous_integration/environment.yaml index 5338027..89a879c 100644 --- a/continuous_integration/environment.yaml +++ b/continuous_integration/environment.yaml @@ -11,6 +11,7 @@ dependencies: - packaging - pytest - pytest-cov + - xarray - pip - pip: - docutils diff --git a/doc/source/conf.py b/doc/source/conf.py index 388725e..0d065c6 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -11,7 +11,6 @@ # All configuration values have a default; values that are commented out # serve to show the default. -import sys, os # If extensions (or modules to document with autodoc) are in another directory, # add these directories to sys.path here. If the directory is relative to the @@ -25,33 +24,33 @@ # Add any Sphinx extension module names here, as strings. They can be extensions # coming with Sphinx (named 'sphinx.ext.*') or your custom ones. -extensions = ['sphinx.ext.autodoc', 'sphinx.ext.doctest', 'sphinx.ext.todo', - 'sphinx.ext.inheritance_diagram'] +extensions = ["sphinx.ext.autodoc", "sphinx.ext.doctest", "sphinx.ext.todo", + "sphinx.ext.inheritance_diagram"] # Add any paths that contain templates here, relative to this directory. -templates_path = ['_templates'] +templates_path = ["_templates"] # The suffix of source filenames. -source_suffix = '.rst' +source_suffix = ".rst" # The encoding of source files. #source_encoding = 'utf-8-sig' # The master toctree document. -master_doc = 'index' +master_doc = "index" # General information about the project. -project = u'pygac' -copyright = u'2014, Abhay Devasthale, Martin Raspaud and Adam Dybbroe' +project = u"pygac" +copyright = u"2014, Abhay Devasthale, Martin Raspaud and Adam Dybbroe" # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the # built documents. # # The short X.Y version. -version = '0.1' +version = "0.1" # The full version, including alpha/beta/rc tags. -release = '0.1' +release = "0.1" # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. @@ -82,7 +81,7 @@ #show_authors = False # The name of the Pygments (syntax highlighting) style to use. -pygments_style = 'sphinx' +pygments_style = "sphinx" # A list of ignored prefixes for module index sorting. #modindex_common_prefix = [] @@ -92,7 +91,7 @@ # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. -html_theme = 'default' +html_theme = "default" # Theme options are theme-specific and customize the look and feel of a theme # further. For a list of options available for each theme, see the @@ -121,7 +120,7 @@ # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ['_static'] +html_static_path = ["_static"] # If not '', a 'Last updated on:' timestamp is inserted at every page bottom, # using the given strftime format. @@ -165,7 +164,7 @@ #html_file_suffix = None # Output file base name for HTML help builder. -htmlhelp_basename = 'pygacdoc' +htmlhelp_basename = "pygacdoc" # -- Options for LaTeX output -------------------------------------------------- @@ -179,8 +178,8 @@ # Grouping the document tree into LaTeX files. List of tuples # (source start file, target name, title, author, documentclass [howto/manual]). latex_documents = [ - ('index', 'pygac.tex', u'pygac Documentation', - u'Abhay Devasthale, Martin Raspaud and Adam Dybbroe', 'manual'), + ("index", "pygac.tex", u"pygac Documentation", + u"Abhay Devasthale, Martin Raspaud and Adam Dybbroe", "manual"), ] # The name of an image file (relative to this directory) to place at the top of @@ -212,6 +211,6 @@ # One entry per manual page. List of tuples # (source start file, name, description, authors, manual section). man_pages = [ - ('index', 'pygac', u'pygac Documentation', - [u'Abhay Devasthale, Martin Raspaud and Adam Dybbroe'], 1) + ("index", "pygac", u"pygac Documentation", + [u"Abhay Devasthale, Martin Raspaud and Adam Dybbroe"], 1) ] diff --git a/pygac/__init__.py b/pygac/__init__.py index 8c97fc8..0a984b2 100644 --- a/pygac/__init__.py +++ b/pygac/__init__.py @@ -30,10 +30,11 @@ "install').") import logging + from pygac.configuration import get_config, read_config_file # noqa from pygac.runner import get_reader_class, process_file # noqa # add a NullHandler to prevent messages in sys.stderr if the using application does # not use logging, but pygac makes logging calls of severity WARNING and greater. # See https://docs.python.org/3/howto/logging.html (Configuring Logging for a Library) -logging.getLogger('pygac').addHandler(logging.NullHandler()) +logging.getLogger("pygac").addHandler(logging.NullHandler()) diff --git a/pygac/calibration.py b/pygac/calibration/noaa.py similarity index 86% rename from pygac/calibration.py rename to pygac/calibration/noaa.py index 1f19cf6..4ca20d1 100644 --- a/pygac/calibration.py +++ b/pygac/calibration/noaa.py @@ -40,14 +40,78 @@ LOG = logging.getLogger(__name__) +def calibrate(ds, custom_coeffs=None, coeffs_file=None): + """Calibrate the dataset `ds`, possibly using custom_coeffs or a coeff file. + + Args: + ds: xarray Dataset object containing the data to calibrate. + custom_calibration: dictionary with a subset of user defined satellite specific + calibration coefficients + calibration_file: path to json file containing default calibrations + """ + channels = ds["channels"].data + times = ds.coords["times"] + scan_line_numbers = ds["scan_line_index"].data + + calibration_coeffs = Calibrator( + ds.attrs["spacecraft_name"], + custom_coeffs=custom_coeffs, + coeffs_file=coeffs_file + ) + + # `times` is in nanosecond resolution. However, convertion from datetime64[ns] to datetime objects + # does not work, and anyway only the date is needed. + start_time = times[0].dt + year = start_time.year.item() + jday = start_time.dayofyear.item() + + + corr = ds.attrs["sun_earth_distance_correction_factor"] + + # 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, + calibration_coeffs, + corr + ) + + prt = ds["prt_counts"].data + ict = ds["ict_counts"].data + space = ds["space_counts"].data + + ir_channels_to_calibrate = [3, 4, 5] + + for chan in ir_channels_to_calibrate: + channels[:, :, chan - 6] = calibrate_thermal( + channels[:, :, chan - 6], + prt, + ict[:, chan - 3], + space[:, chan - 3], + scan_line_numbers, + chan, + calibration_coeffs + ) + + new_ds = ds.copy() + new_ds["channels"].data = channels + + new_ds.attrs["calib_coeffs_version"] = calibration_coeffs.version + + return new_ds + + class CoeffStatus(Enum): """Indicates the status of calibration coefficients.""" - NOMINAL = 'nominal' - PROVISIONAL = 'provisional' - EXPERIMENTAL = 'experimental' + NOMINAL = "nominal" + PROVISIONAL = "provisional" + EXPERIMENTAL = "experimental" -class Calibrator(object): +class Calibrator: """Factory class to create namedtuples holding the calibration coefficients. Attributes: @@ -56,17 +120,17 @@ class Calibrator(object): default_coeffs: dictonary containing default values for all spacecrafts """ version_hashs = { - '963af9b66268475ed500ad7b37da33c5': { - 'name': 'PATMOS-x, v2017r1', - 'status': CoeffStatus.NOMINAL + "963af9b66268475ed500ad7b37da33c5": { + "name": "PATMOS-x, v2017r1", + "status": CoeffStatus.NOMINAL }, - '689386c822de18a07194ac7fd71652ea': { - 'name': 'PATMOS-x, v2017r1, with provisional coefficients for MetOp-C', - 'status': CoeffStatus.PROVISIONAL + "689386c822de18a07194ac7fd71652ea": { + "name": "PATMOS-x, v2017r1, with provisional coefficients for MetOp-C", + "status": CoeffStatus.PROVISIONAL }, - 'e8735ec394ecdb87b7edcd261e72d2eb': { - 'name': 'PATMOS-x, v2023', - 'status': CoeffStatus.PROVISIONAL + "e8735ec394ecdb87b7edcd261e72d2eb": { + "name": "PATMOS-x, v2023", + "status": CoeffStatus.PROVISIONAL }, } fields = [ @@ -75,7 +139,7 @@ class Calibrator(object): "to_eff_blackbody_slope", "date_of_launch", "d", "spacecraft", "version" ] - Calibrator = namedtuple('Calibrator', fields) + Calibrator = namedtuple("Calibrator", fields) default_coeffs = None default_file = None default_version = None @@ -112,21 +176,21 @@ def __new__(cls, spacecraft, custom_coeffs=None, coeffs_file=None): for key in ("dark_count", "gain_switch", "s0", "s1", "s2"): arraycoeffs[key] = np.array([ coeffs[channel][key] - for channel in ('channel_1', 'channel_2', 'channel_3a') + for channel in ("channel_1", "channel_2", "channel_3a") ], dtype=float) # thermal channels for key in ("centroid_wavenumber", "space_radiance", "to_eff_blackbody_intercept", "to_eff_blackbody_slope"): arraycoeffs[key] = np.array([ coeffs[channel][key] - for channel in ('channel_3b', 'channel_4', 'channel_5') + for channel in ("channel_3b", "channel_4", "channel_5") ], dtype=float) arraycoeffs["b"] = np.array([ [ coeffs[channel][key] for key in ("b0", "b1", "b2") ] - for channel in ('channel_3b', 'channel_4', 'channel_5') + for channel in ("channel_3b", "channel_4", "channel_5") ], dtype=float) # thermometers # Note, that "thermometer_0" does not exists, and is filled with zeros to @@ -139,7 +203,7 @@ def __new__(cls, spacecraft, custom_coeffs=None, coeffs_file=None): for d in range(5) ], dtype=float) # parse date of launch - date_of_launch_str = coeffs["date_of_launch"].replace('Z', '+00:00') + date_of_launch_str = coeffs["date_of_launch"].replace("Z", "+00:00") if sys.version_info < (3, 7): # Note that here any time information is lost import dateutil.parser @@ -174,7 +238,7 @@ def read_coeffs(cls, coeffs_file): else: LOG.debug("Read PyGAC internal calibration coefficients.") coeffs_file = files("pygac") / "data/calibration.json" - with open(coeffs_file, mode='rb') as json_file: + with open(coeffs_file, mode="rb") as json_file: content = json_file.read() coeffs = json.loads(content) version = cls._get_coeffs_version(content) @@ -187,10 +251,10 @@ def _get_coeffs_version(cls, coeff_file_content): digest = md5_hash.hexdigest() version_dict = cls.version_hashs.get( digest, - {'name': None, 'status': None} + {"name": None, "status": None} ) - version = version_dict['name'] - status = version_dict['status'] + version = version_dict["name"] + status = version_dict["status"] if version is None: warning = "Unknown calibration coefficients version!" warnings.warn(warning, RuntimeWarning) @@ -199,7 +263,7 @@ def _get_coeffs_version(cls, coeff_file_content): LOG.info('Identified calibration coefficients version "%s".', version) if status != CoeffStatus.NOMINAL: - warning = 'Using {} calibration coefficients'.format(status) + warning = "Using {} calibration coefficients".format(status) warnings.warn(warning, RuntimeWarning) LOG.warning(warning) return version @@ -393,14 +457,15 @@ def calibrate_thermal(counts, prt, ict, space, line_numbers, channel, cal): # Note that the prt values are the average value of the three readings from one of the four # PRTs. See reader.get_telemetry implementations. prt_threshold = 50 # empirically found and set by Abhay Devasthale - offset = 0 - for i, prt_val in enumerate(prt): + for offset in range(5): # According to the KLM Guide the fill value between PRT measurments is 0, but we search - # for the first measurment gap using the threshold. Is this on purpose? - if prt_val < prt_threshold: - offset = i + # for the first measurement gap using the threshold, because the fill value is in practice + # not always exactly 0. + if np.median(prt[(line_numbers - line_numbers[0]) % 5 == offset]) < prt_threshold: break + else: + raise IndexError("No PRT 0-index found!") # get the PRT index, iprt equals to 0 corresponds to the measurement gaps iprt = (line_numbers - line_numbers[0] + 5 - offset) % 5 @@ -450,16 +515,18 @@ def calibrate_thermal(counts, prt, ict, space, line_numbers, channel, cal): if channel == 3: zeros = ict < ict_threshold nonzeros = np.logical_not(zeros) - - ict[zeros] = np.interp((zeros).nonzero()[0], - (nonzeros).nonzero()[0], - ict[nonzeros]) + try: + ict[zeros] = np.interp((zeros).nonzero()[0], + (nonzeros).nonzero()[0], + ict[nonzeros]) + except ValueError: # 3b has no valid data + return counts zeros = space < space_threshold nonzeros = np.logical_not(zeros) space[zeros] = np.interp((zeros).nonzero()[0], - (nonzeros).nonzero()[0], - space[nonzeros]) + (nonzeros).nonzero()[0], + space[nonzeros]) # convolving and smoothing PRT, ICT and SPACE values if lines > 51: @@ -468,9 +535,9 @@ def calibrate_thermal(counts, prt, ict, space, line_numbers, channel, cal): wlength = 3 weighting_function = np.ones(wlength, dtype=float) / wlength - tprt_convolved = np.convolve(tprt, weighting_function, 'same') - ict_convolved = np.convolve(ict, weighting_function, 'same') - space_convolved = np.convolve(space, weighting_function, 'same') + tprt_convolved = np.convolve(tprt, weighting_function, "same") + ict_convolved = np.convolve(ict, weighting_function, "same") + space_convolved = np.convolve(space, weighting_function, "same") # take care of the beginning and end tprt_convolved[0:(wlength - 1) // 2] = tprt_convolved[(wlength - 1) // 2] @@ -491,7 +558,7 @@ def calibrate_thermal(counts, prt, ict, space, line_numbers, channel, cal): # TsBB = A + B*TBB (7.1.2.4-3) # NBB = c1*nu_e^3/(exp(c2*nu_e/TsBB) - 1) (7.1.2.4-3) # where the constants of the Planck function are defined as c1 = 2*h*c^2, c2 = h*c/k_B. - # constatns + # constants c1 = 1.1910427e-5 # mW/m^2/sr/cm^{-4} c2 = 1.4387752 # cm K # coefficients diff --git a/pygac/configuration.py b/pygac/configuration.py index 9ca301e..f808d3e 100644 --- a/pygac/configuration.py +++ b/pygac/configuration.py @@ -27,6 +27,7 @@ import logging import os import sys + try: import configparser except ImportError: @@ -42,7 +43,7 @@ class FileNotFoundError(OSError): class Configuration(configparser.ConfigParser, object): """Configuration container for pygac.""" - config_file = '' + config_file = "" def read(self, config_file): """Read and parse the configuration file @@ -69,8 +70,8 @@ def read(self, config_file): def get(self, *args, **kwargs): """python 2 compatibility for fallback attribute""" if sys.version_info.major < 3: - if 'fallback' in kwargs: - fallback = kwargs.pop('fallback') + if "fallback" in kwargs: + fallback = kwargs.pop("fallback") else: fallback = None try: @@ -99,7 +100,7 @@ def get_config(initialized=True): try: config_file = os.environ["PYGAC_CONFIG_FILE"] except KeyError: - LOG.error('Environment variable PYGAC_CONFIG_FILE not set!') + LOG.error("Environment variable PYGAC_CONFIG_FILE not set!") raise _config.read(config_file) return _config diff --git a/pygac/correct_tsm_issue.py b/pygac/correct_tsm_issue.py index f940125..394a215 100644 --- a/pygac/correct_tsm_issue.py +++ b/pygac/correct_tsm_issue.py @@ -23,10 +23,11 @@ an orbit contain corrupt data. This module identifies affected pixels and flags them with fill values.""" -import numpy as np -import bottleneck as bn import datetime +import bottleneck as bn +import numpy as np + TSM_AFFECTED_INTERVALS_POD = { 3: [(datetime.datetime(2001, 10, 19, 4, 50), datetime.datetime(2001, 10, 19, 13, 38)), (datetime.datetime(2001, 10, 19, 16, 58), @@ -409,7 +410,7 @@ def std_filter(data, box_size): # need to surround the data with NaNs to calculate values at the boundary padded_data = np.pad( data, (border, border), - mode='constant', + mode="constant", constant_values=np.nan ) windows = _rolling_window(padded_data, size) diff --git a/pygac/gac_io.py b/pygac/gac_io.py index 8b08b2d..69110bc 100644 --- a/pygac/gac_io.py +++ b/pygac/gac_io.py @@ -33,7 +33,7 @@ import h5py import numpy as np -from pygac.utils import slice_channel, strip_invalid_lat, check_user_scanlines +from pygac.utils import check_user_scanlines, slice_channel, strip_invalid_lat LOG = logging.getLogger(__name__) @@ -52,20 +52,20 @@ def save_gac(satellite_name, gac_file, meta_data, output_file_prefix, avhrr_dir, qual_dir, sunsatangles_dir): - midnight_scanline = meta_data['midnight_scanline'] - miss_lines = meta_data['missing_scanlines'] - corr = meta_data['sun_earth_distance_correction_factor'] + midnight_scanline = meta_data["midnight_scanline"] + miss_lines = meta_data["missing_scanlines"] + corr = meta_data["sun_earth_distance_correction_factor"] last_scan_line_number = qual_flags[-1, 0] # Strip invalid coordinates first_valid_lat, last_valid_lat = strip_invalid_lat(lats) if first_valid_lat > start_line: - LOG.info('New start_line chosen (due to invalid lat/lon ' - 'info) = ' + str(first_valid_lat)) + LOG.info("New start_line chosen (due to invalid lat/lon " + "info) = " + str(first_valid_lat)) if end_line > last_valid_lat: - LOG.info('New end_line chosen (due to invalid lat/lon ' - 'info) = ' + str(last_valid_lat)) + LOG.info("New end_line chosen (due to invalid lat/lon " + "info) = " + str(last_valid_lat)) # Check user-defined scanlines start_line, end_line = check_user_scanlines( @@ -209,35 +209,35 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, t_obj = time.strptime(enddate + endtime[0:6], "%Y%m%d%H%M%S") endtime_sec1970 = calendar.timegm(t_obj) - LOG.info('Output file prefix = ' + str(output_file_prefix)) - LOG.info('AVHRR data will be written to ' + str(avhrr_dir)) - ofn = os.path.join(avhrr_dir, (output_file_prefix + '_avhrr_' + - satellite_name + '_99999_' + - startdate + 'T' + starttime + 'Z_' + - enddate + 'T' + endtime + 'Z.h5')) + LOG.info("Output file prefix = " + str(output_file_prefix)) + LOG.info("AVHRR data will be written to " + str(avhrr_dir)) + ofn = os.path.join(avhrr_dir, (output_file_prefix + "_avhrr_" + + satellite_name + "_99999_" + + startdate + "T" + starttime + "Z_" + + enddate + "T" + endtime + "Z.h5")) - LOG.info('Filename: ' + str(os.path.basename(ofn))) + LOG.info("Filename: " + str(os.path.basename(ofn))) fout = h5py.File(ofn, "w") - dset1 = fout.create_dataset("/image1/data", dtype='int16', data=ref1) - dset2 = fout.create_dataset("/image2/data", dtype='int16', data=ref2) - dset3 = fout.create_dataset("/image3/data", dtype='int16', data=bt3) - dset4 = fout.create_dataset("/image4/data", dtype='int16', data=bt4) - dset5 = fout.create_dataset("/image5/data", dtype='int16', data=bt5) - dset6 = fout.create_dataset("/image6/data", dtype='int16', data=ref3) - dset7 = fout.create_dataset("/where/lat/data", dtype='int32', + dset1 = fout.create_dataset("/image1/data", dtype="int16", data=ref1) + dset2 = fout.create_dataset("/image2/data", dtype="int16", data=ref2) + dset3 = fout.create_dataset("/image3/data", dtype="int16", data=bt3) + dset4 = fout.create_dataset("/image4/data", dtype="int16", data=bt4) + dset5 = fout.create_dataset("/image5/data", dtype="int16", data=bt5) + dset6 = fout.create_dataset("/image6/data", dtype="int16", data=ref3) + dset7 = fout.create_dataset("/where/lat/data", dtype="int32", data=arrLat_full) - dset8 = fout.create_dataset("/where/lon/data", dtype='int32', + dset8 = fout.create_dataset("/where/lon/data", dtype="int32", data=arrLon_full) del dset8 channellist = [] - channellist.append("channel1".encode('utf8')) - channellist.append("channel2".encode('utf8')) - channellist.append("channel3b".encode('utf8')) - channellist.append("channel4".encode('utf8')) - channellist.append("channel5".encode('utf8')) - channellist.append("channel3a".encode('utf8')) + channellist.append("channel1".encode("utf8")) + channellist.append("channel2".encode("utf8")) + channellist.append("channel3b".encode("utf8")) + channellist.append("channel4".encode("utf8")) + channellist.append("channel5".encode("utf8")) + channellist.append("channel3a".encode("utf8")) dset10 = fout.create_dataset("/how/channel_list", data=channellist) del dset10 @@ -282,8 +282,8 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g1.attrs["product"] = np.string_("SATCH") g1.attrs["quantity"] = np.string_("REFL") - g1.attrs["dataset_name"] = np.string_('Channel 1 reflectance') - g1.attrs["units"] = np.string_('%') + g1.attrs["dataset_name"] = np.string_("Channel 1 reflectance") + g1.attrs["units"] = np.string_("%") g1.attrs["gain"] = np.float32(0.01) g1.attrs["offset"] = np.float32(0.0) g1.attrs["missingdata"] = np.int32(MISSING_DATA) @@ -295,8 +295,8 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g2.attrs["product"] = np.string_("SATCH") g2.attrs["quantity"] = np.string_("REFL") - g2.attrs["dataset_name"] = np.string_('Channel 2 reflectance') - g2.attrs["units"] = np.string_('%') + g2.attrs["dataset_name"] = np.string_("Channel 2 reflectance") + g2.attrs["units"] = np.string_("%") g2.attrs["gain"] = np.float32(0.01) g2.attrs["offset"] = np.float32(0.0) g2.attrs["missingdata"] = np.int32(MISSING_DATA) @@ -308,8 +308,8 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g6.attrs["product"] = np.string_("SATCH") g6.attrs["quantity"] = np.string_("REFL") - g6.attrs["dataset_name"] = np.string_('Channel 3a reflectance') - g6.attrs["units"] = np.string_('%') + g6.attrs["dataset_name"] = np.string_("Channel 3a reflectance") + g6.attrs["units"] = np.string_("%") g6.attrs["gain"] = np.float32(0.01) g6.attrs["offset"] = np.float32(0.0) g6.attrs["missingdata"] = np.int32(MISSING_DATA) @@ -321,8 +321,8 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g3.attrs["product"] = np.string_("SATCH") g3.attrs["quantity"] = np.string_("TB") - g3.attrs["dataset_name"] = np.string_('Channel 3b brightness temperature') - g3.attrs["units"] = np.string_('K') + g3.attrs["dataset_name"] = np.string_("Channel 3b brightness temperature") + g3.attrs["units"] = np.string_("K") g3.attrs["gain"] = np.float32(0.01) g3.attrs["offset"] = np.float32(273.15) g3.attrs["missingdata"] = np.int32(MISSING_DATA) @@ -334,8 +334,8 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g4.attrs["product"] = np.string_("SATCH") g4.attrs["quantity"] = np.string_("TB") - g4.attrs["dataset_name"] = np.string_('Channel 4 brightness temperature') - g4.attrs["units"] = np.string_('K') + g4.attrs["dataset_name"] = np.string_("Channel 4 brightness temperature") + g4.attrs["units"] = np.string_("K") g4.attrs["gain"] = np.float32(0.01) g4.attrs["offset"] = np.float32(273.15) g4.attrs["missingdata"] = np.int32(MISSING_DATA) @@ -347,8 +347,8 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g5.attrs["product"] = np.string_("SATCH") g5.attrs["quantity"] = np.string_("TB") - g5.attrs["dataset_name"] = np.string_('Channel 5 brightness temperature') - g5.attrs["units"] = np.string_('K') + g5.attrs["dataset_name"] = np.string_("Channel 5 brightness temperature") + g5.attrs["units"] = np.string_("K") g5.attrs["gain"] = np.float32(0.01) g5.attrs["offset"] = np.float32(273.15) g5.attrs["missingdata"] = np.int32(MISSING_DATA) @@ -358,8 +358,8 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g5.attrs["startdate"] = np.string_(startdate) g5.attrs["enddate"] = np.string_(enddate) - g7.attrs["dataset_name"] = np.string_('Latitude') - g7.attrs["units"] = np.string_('Deg') + g7.attrs["dataset_name"] = np.string_("Latitude") + g7.attrs["units"] = np.string_("Deg") g7.attrs["gain"] = np.float32(0.0010) g7.attrs["offset"] = np.float32(0.0) g7.attrs["missingdata"] = np.int32(MISSING_DATA_LATLON) @@ -369,8 +369,8 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g7.attrs["startdate"] = np.string_(startdate) g7.attrs["enddate"] = np.string_(enddate) - g8.attrs["dataset_name"] = np.string_('Longitude') - g8.attrs["units"] = np.string_('Deg') + g8.attrs["dataset_name"] = np.string_("Longitude") + g8.attrs["units"] = np.string_("Deg") g8.attrs["gain"] = np.float32(0.0010) g8.attrs["offset"] = np.float32(0.0) g8.attrs["missingdata"] = np.int32(MISSING_DATA_LATLON) @@ -419,25 +419,25 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, fout.close() - LOG.info('Sun and Satellite viewing angles will be ' + - 'written to ' + str(sunsatangles_dir)) + LOG.info("Sun and Satellite viewing angles will be " + + "written to " + str(sunsatangles_dir)) ofn = os.path.join(sunsatangles_dir, - (output_file_prefix + '_sunsatangles_' + - satellite_name + '_99999_' + startdate + - 'T' + starttime + 'Z_' + - enddate + 'T' + endtime + 'Z.h5')) + (output_file_prefix + "_sunsatangles_" + + satellite_name + "_99999_" + startdate + + "T" + starttime + "Z_" + + enddate + "T" + endtime + "Z.h5")) - LOG.info('Filename: ' + str(os.path.basename(ofn))) + LOG.info("Filename: " + str(os.path.basename(ofn))) fout = h5py.File(ofn, "w") - dset1 = fout.create_dataset("/image1/data", dtype='int16', data=arrSZA) - dset2 = fout.create_dataset("/image2/data", dtype='int16', data=arrSTZ) - dset3 = fout.create_dataset("/image3/data", dtype='int16', data=arrRAA) - dset4 = fout.create_dataset("/image4/data", dtype='int16', data=arrSAA) - dset5 = fout.create_dataset("/image5/data", dtype='int16', data=arrSTA) - dset6 = fout.create_dataset("/where/lat/data", dtype='int32', + dset1 = fout.create_dataset("/image1/data", dtype="int16", data=arrSZA) + dset2 = fout.create_dataset("/image2/data", dtype="int16", data=arrSTZ) + dset3 = fout.create_dataset("/image3/data", dtype="int16", data=arrRAA) + dset4 = fout.create_dataset("/image4/data", dtype="int16", data=arrSAA) + dset5 = fout.create_dataset("/image5/data", dtype="int16", data=arrSTA) + dset6 = fout.create_dataset("/where/lat/data", dtype="int32", data=arrLat_full) - dset7 = fout.create_dataset("/where/lon/data", dtype='int32', + dset7 = fout.create_dataset("/where/lon/data", dtype="int32", data=arrLon_full) del dset4, dset5, dset6, dset7 @@ -450,12 +450,12 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g5 = fout.require_group("/image5") g6 = fout.require_group("/where") - g1.attrs["description"] = np.string_('Solar zenith angle') - g2.attrs["description"] = np.string_('Satellite zenith angle') + g1.attrs["description"] = np.string_("Solar zenith angle") + g2.attrs["description"] = np.string_("Satellite zenith angle") g3.attrs["description"] = np.string_( - 'Relative satellite-sun azimuth angle') - g4.attrs["description"] = np.string_('Solar azimuth angle') - g5.attrs["description"] = np.string_('Satellite azimuth angle') + "Relative satellite-sun azimuth angle") + g4.attrs["description"] = np.string_("Solar azimuth angle") + g5.attrs["description"] = np.string_("Satellite azimuth angle") g6.attrs["num_of_pixels"] = np.int32(arrSZA.shape[1]) g6.attrs["num_of_lines"] = np.int32(arrSZA.shape[0]) g6.attrs["xscale"] = np.float32(0.0) @@ -476,8 +476,8 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g1.attrs["product"] = np.string_("SUNZ") g1.attrs["quantity"] = np.string_("DEG") - g1.attrs["dataset_name"] = np.string_('Solar zenith angle') - g1.attrs["units"] = np.string_('Deg') + g1.attrs["dataset_name"] = np.string_("Solar zenith angle") + g1.attrs["units"] = np.string_("Deg") g1.attrs["gain"] = np.float32(0.01) g1.attrs["offset"] = np.float32(0.0) g1.attrs["missingdata"] = np.int32(MISSING_DATA) @@ -489,8 +489,8 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g2.attrs["product"] = np.string_("SATZ") g2.attrs["quantity"] = np.string_("DEG") - g2.attrs["dataset_name"] = np.string_('Satellite zenith angle') - g2.attrs["units"] = np.string_('Deg') + g2.attrs["dataset_name"] = np.string_("Satellite zenith angle") + g2.attrs["units"] = np.string_("Deg") g2.attrs["gain"] = np.float32(0.01) g2.attrs["offset"] = np.float32(0.0) g2.attrs["missingdata"] = np.int32(MISSING_DATA) @@ -503,8 +503,8 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g3.attrs["product"] = np.string_("SSAZD") g3.attrs["quantity"] = np.string_("DEG") g3.attrs["dataset_name"] = np.string_( - 'Relative satellite-sun azimuth angle') - g3.attrs["units"] = np.string_('Deg') + "Relative satellite-sun azimuth angle") + g3.attrs["units"] = np.string_("Deg") g3.attrs["gain"] = np.float32(0.01) g3.attrs["offset"] = np.float32(0.0) g3.attrs["missingdata"] = np.int32(MISSING_DATA) @@ -516,8 +516,8 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g4.attrs["product"] = np.string_("SUNA") g4.attrs["quantity"] = np.string_("DEG") - g4.attrs["dataset_name"] = np.string_('Solar azimuth angle') - g4.attrs["units"] = np.string_('Deg') + g4.attrs["dataset_name"] = np.string_("Solar azimuth angle") + g4.attrs["units"] = np.string_("Deg") g4.attrs["gain"] = np.float32(0.01) g4.attrs["offset"] = np.float32(0.0) g4.attrs["missingdata"] = np.int32(MISSING_DATA) @@ -529,8 +529,8 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g5.attrs["product"] = np.string_("SATA") g5.attrs["quantity"] = np.string_("DEG") - g5.attrs["dataset_name"] = np.string_('Satellite azimuth angle') - g5.attrs["units"] = np.string_('Deg') + g5.attrs["dataset_name"] = np.string_("Satellite azimuth angle") + g5.attrs["units"] = np.string_("Deg") g5.attrs["gain"] = np.float32(0.01) g5.attrs["offset"] = np.float32(0.0) g5.attrs["missingdata"] = np.int32(MISSING_DATA) @@ -540,8 +540,8 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g5.attrs["startdate"] = np.string_(startdate) g5.attrs["enddate"] = np.string_(enddate) - g6.attrs["dataset_name"] = np.string_('Latitude') - g6.attrs["units"] = np.string_('Deg') + g6.attrs["dataset_name"] = np.string_("Latitude") + g6.attrs["units"] = np.string_("Deg") g6.attrs["gain"] = np.float32(0.0010) g6.attrs["offset"] = np.float32(0.0) g6.attrs["missingdata"] = np.int32(MISSING_DATA_LATLON) @@ -551,8 +551,8 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g6.attrs["startdate"] = np.string_(startdate) g6.attrs["enddate"] = np.string_(enddate) - g7.attrs["dataset_name"] = np.string_('Longitude') - g7.attrs["units"] = np.string_('Deg') + g7.attrs["dataset_name"] = np.string_("Longitude") + g7.attrs["units"] = np.string_("Deg") g7.attrs["gain"] = np.float32(0.0010) g7.attrs["offset"] = np.float32(0.0) g7.attrs["missingdata"] = np.int32(MISSING_DATA_LATLON) @@ -583,25 +583,25 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, fout.close() - LOG.info('Quality flags will be ' + - 'written to ' + str(qual_dir)) + LOG.info("Quality flags will be " + + "written to " + str(qual_dir)) ofn = os.path.join(qual_dir, - (output_file_prefix + '_qualflags_' + - satellite_name + '_99999_' + startdate + - 'T' + starttime + 'Z_' + - enddate + 'T' + endtime + 'Z.h5')) + (output_file_prefix + "_qualflags_" + + satellite_name + "_99999_" + startdate + + "T" + starttime + "Z_" + + enddate + "T" + endtime + "Z.h5")) - LOG.info('Filename: ' + str(os.path.basename(ofn))) + LOG.info("Filename: " + str(os.path.basename(ofn))) fout = h5py.File(ofn, "w") g1 = fout.require_group("/qual_flags") - dset1 = g1.create_dataset("data", dtype='int16', data=qual_flags) + dset1 = g1.create_dataset("data", dtype="int16", data=qual_flags) del dset1 g1.attrs["product"] = np.string_("QFLAG") g1.attrs["quantity"] = np.string_("INT") - g1.attrs["dataset_name"] = np.string_('Scanline quality flags') - g1.attrs["units"] = np.string_('None') + g1.attrs["dataset_name"] = np.string_("Scanline quality flags") + g1.attrs["units"] = np.string_("None") g1.attrs["gain"] = np.int32(1) g1.attrs["offset"] = np.int32(0) g1.attrs["missingdata"] = np.int32(MISSING_DATA) @@ -615,13 +615,13 @@ def avhrrGAC_io(satellite_name, xutcs, startdate, enddate, starttime, endtime, g1.attrs["last_scan_line_number"] = last_scan_line_number g2 = fout.require_group("/ancillary") - dset2 = g2.create_dataset("missing_scanlines", dtype='int16', + dset2 = g2.create_dataset("missing_scanlines", dtype="int16", data=miss_lines) del dset2 - dset3 = g2.create_dataset("scanline_timestamps", dtype='int64', - data=xutcs.astype('int64')) - dset3.attrs['units'] = 'Milliseconds since 1970-01-01 00:00:00 UTC' - dset3.attrs['calendar'] = 'standard' + dset3 = g2.create_dataset("scanline_timestamps", dtype="int64", + data=xutcs.astype("int64")) + dset3.attrs["units"] = "Milliseconds since 1970-01-01 00:00:00 UTC" + dset3.attrs["calendar"] = "standard" g2.attrs["midnight_scanline"] = np.string_(midnight_scanline) fout.close() diff --git a/pygac/gac_klm.py b/pygac/gac_klm.py index 081da8d..7d08697 100644 --- a/pygac/gac_klm.py +++ b/pygac/gac_klm.py @@ -127,16 +127,17 @@ ("spacecraft_altitude_above_reference_ellipsoid", ">u2"), ("angular_relationships", ">i2", (153, )), ("zero_fill3", ">i2", (3, )), - ("earth_location", ">i4", (102, )), + ("earth_location", [("lats", ">i2"), + ("lons", ">i2")], (51,)), ("zero_fill4", ">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'), ]), + ("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"), diff --git a/pygac/gac_pod.py b/pygac/gac_pod.py index 7700229..bc81416 100644 --- a/pygac/gac_pod.py +++ b/pygac/gac_pod.py @@ -30,8 +30,8 @@ import numpy as np -from pygac.pod_reader import PODReader, main_pod from pygac.gac_reader import GACReader +from pygac.pod_reader import PODReader, main_pod LOG = logging.getLogger(__name__) @@ -42,7 +42,8 @@ ("number_of_meaningful_zenith_angles_and_earth_location_appended", ">u1"), ("solar_zenith_angles", "i1", (51, )), - ("earth_location", ">i2", (102, )), + ("earth_location", [("lats", ">i2"), + ("lons", ">i2")], (51,)), ("telemetry", ">u4", (35, )), ("sensor_data", ">u4", (682, )), ("add_on_zenith", ">u2", (10, )), diff --git a/pygac/gac_reader.py b/pygac/gac_reader.py index 98f8706..52ec6cd 100644 --- a/pygac/gac_reader.py +++ b/pygac/gac_reader.py @@ -28,9 +28,9 @@ import logging -from pygac.reader import Reader, ReaderError import pygac.pygac_geotiepoints as gtp - +from pygac.pygac_geotiepoints import GAC_LONLAT_SAMPLE_POINTS +from pygac.reader import Reader, ReaderError LOG = logging.getLogger(__name__) @@ -42,10 +42,11 @@ class GACReader(Reader): scan_freq = 2.0 / 1000.0 # Max scanlines max_scanlines = 15000 + lonlat_sample_points = GAC_LONLAT_SAMPLE_POINTS def __init__(self, *args, **kwargs): """Init the GAC reader.""" - super(GACReader, self).__init__(*args, **kwargs) + super().__init__(*args, **kwargs) self.scan_width = 409 self.lonlat_interpolator = gtp.gac_lat_lon_interpolator @@ -55,9 +56,9 @@ def _validate_header(cls, header): # call super to enter the Method Resolution Order (MRO) super(GACReader, cls)._validate_header(header) LOG.debug("validate header") - data_set_name = header['data_set_name'].decode() + data_set_name = header["data_set_name"].decode() # split header into parts creation_site, transfer_mode, platform_id = ( - data_set_name.split('.')[:3]) - if transfer_mode != 'GHRR': + data_set_name.split(".")[:3]) + if transfer_mode != "GHRR": raise ReaderError('Improper transfer mode "%s"!' % transfer_mode) diff --git a/pygac/klm_reader.py b/pygac/klm_reader.py index da1f0e0..799056f 100644 --- a/pygac/klm_reader.py +++ b/pygac/klm_reader.py @@ -36,6 +36,7 @@ import datetime import logging + try: from enum import IntFlag except ImportError: @@ -579,60 +580,60 @@ class KLM_QualityIndicator(IntFlag): ("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'), +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, )), + ("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'), + ("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') + ("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', - 13: 'metopc', + spacecraft_names = {4: "noaa15", + 2: "noaa16", + 6: "noaa17", + 7: "noaa18", + 8: "noaa19", + 12: "metopa", + 11: "metopb", + 13: "metopc", } - spacecrafts_orbital = {4: 'noaa 15', - 2: 'noaa 16', - 6: 'noaa 17', - 7: 'noaa 18', - 8: 'noaa 19', - 12: 'metop 02', - 11: 'metop 01', - 13: 'metop 03', + spacecrafts_orbital = {4: "noaa 15", + 2: "noaa 16", + 6: "noaa 17", + 7: "noaa 18", + 8: "noaa 19", + 12: "metop 02", + 11: "metop 01", + 13: "metop 03", } tsm_affected_intervals = TSM_AFFECTED_INTERVALS_KLM @@ -659,7 +660,7 @@ def read(self, filename, fileobj=None): # file objects to (io.FileIO, io.BufferedReader, io.BufferedWriter) # see: numpy.compat.py3k.isfileobj self.filename = filename - LOG.info('Reading %s', self.filename) + LOG.info("Reading %s", self.filename) with file_opener(fileobj or filename) as fd_: self.ars_head, self.head = self.read_header( filename, fileobj=fd_) @@ -678,7 +679,7 @@ def read(self, filename, fileobj=None): fd_.read(analog_telemetry_v2.itemsize), dtype=analog_telemetry_v2, count=1) # LAC: 1, GAC: 2, ... - self.data_type = self.head['data_type_code'] + self.data_type = self.head["data_type_code"] # read until end of file fd_.seek(self.offset + ars_offset, 0) buffer = fd_.read() @@ -705,7 +706,7 @@ def read_header(cls, filename, fileobj=None): _ars_head, = np.frombuffer( fd_.read(ars_header.itemsize), dtype=ars_header, count=1) - if _ars_head['data_format'].startswith(b'NOAA Level 1b'): + if _ars_head["data_format"].startswith(b"NOAA Level 1b"): ars_head = _ars_head.copy() else: fd_.seek(0) @@ -724,11 +725,11 @@ def _validate_header(cls, header): # call super to enter the Method Resolution Order (MRO) super(KLMReader, cls)._validate_header(header) LOG.debug("validate header") - data_set_name = header['data_set_name'].decode() + data_set_name = header["data_set_name"].decode() # split header into parts creation_site, transfer_mode, platform_id = ( - data_set_name.split('.')[:3]) - allowed_ids = ['NK', 'NL', 'NM', 'NN', 'NP', 'M1', 'M2', 'M3'] + data_set_name.split(".")[:3]) + allowed_ids = ["NK", "NL", "NM", "NN", "NP", "M1", "M2", "M3"] if platform_id not in allowed_ids: raise ReaderError('Improper platform id "%s"!' % platform_id) @@ -741,7 +742,7 @@ def get_telemetry(self): space_counts: np.array """ - prt_counts = np.mean(self.scans["telemetry"]['PRT'], axis=1) + prt_counts = np.mean(self.scans["telemetry"]["PRT"], axis=1) # getting ICT counts @@ -759,10 +760,10 @@ def get_telemetry(self): return prt_counts, ict_counts, space_counts - def _get_lonlat(self): + def _get_lonlat_from_file(self): """Get the longitudes and latitudes.""" - lats = self.scans["earth_location"][:, 0::2] / 1e4 - lons = self.scans["earth_location"][:, 1::2] / 1e4 + lats = self.scans["earth_location"]["lats"] / np.float32(1e4) + lons = self.scans["earth_location"]["lons"] / np.float32(1e4) return lons, lats def get_header_timestamp(self): @@ -775,16 +776,16 @@ def get_header_timestamp(self): 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'] + 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)) + raise ValueError("Corrupt header timestamp: {0}".format(err)) - def _get_times(self): + def _get_times_from_file(self): """Get the times of the scanlines.""" year = self.scans["scan_line_year"] jday = self.scans["scan_line_day_of_year"] @@ -806,16 +807,15 @@ def _get_ir_channels_to_calibrate(self): ir_channels_to_calibrate = [4, 5] return ir_channels_to_calibrate - def postproc(self, channels): + def postproc(self, ds): """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 + ds["channels"].loc[dict(channel_name="3a", scan_line_index=((switch==0) | (switch==2)))] = np.nan + ds["channels"].loc[dict(channel_name="3b", scan_line_index=((switch==1) | (switch==2)))] = np.nan + def _adjust_clock_drift(self): """Adjust the geolocation to compensate for the clock error. diff --git a/pygac/lac_klm.py b/pygac/lac_klm.py index f2a50d1..d202f6d 100644 --- a/pygac/lac_klm.py +++ b/pygac/lac_klm.py @@ -126,16 +126,17 @@ ("spacecraft_altitude_above_reference_ellipsoid", ">u2"), ("angular_relationships", ">i2", (153, )), ("zero_fill2", ">i2", (3, )), - ("earth_location", ">i4", (102, )), + ("earth_location", [("lats", ">i2"), + ("lons", ">i2")], (51,)), ("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'), ]), + ("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"), diff --git a/pygac/lac_pod.py b/pygac/lac_pod.py index 49f6103..8d84dd9 100644 --- a/pygac/lac_pod.py +++ b/pygac/lac_pod.py @@ -28,8 +28,8 @@ import numpy as np -from pygac.pod_reader import PODReader, main_pod from pygac.lac_reader import LACReader +from pygac.pod_reader import PODReader, main_pod LOG = logging.getLogger(__name__) @@ -37,10 +37,10 @@ ("time_code", ">u2", (3, )), ("quality_indicators", ">u4"), ("calibration_coefficients", ">i4", (10, )), - ("number_of_meaningful_zenith_angles_and_earth_location_appended", - ">u1"), + ("number_of_meaningful_zenith_angles_and_earth_location_appended", ">u1"), ("solar_zenith_angles", "i1", (51, )), - ("earth_location", ">i2", (102, )), + ("earth_location", [("lats", ">i2"), + ("lons", ">i2")], (51,)), ("telemetry", ">u4", (35, )), ("sensor_data", ">u4", (3414, )), ("add_on_zenith", ">u2", (10, )), diff --git a/pygac/lac_reader.py b/pygac/lac_reader.py index e3690f4..2cb046b 100644 --- a/pygac/lac_reader.py +++ b/pygac/lac_reader.py @@ -26,9 +26,9 @@ import logging -from pygac.reader import Reader, ReaderError import pygac.pygac_geotiepoints as gtp - +from pygac.pygac_geotiepoints import LAC_LONLAT_SAMPLE_POINTS +from pygac.reader import Reader, ReaderError LOG = logging.getLogger(__name__) @@ -40,6 +40,7 @@ class LACReader(Reader): scan_freq = 6.0 / 1000.0 # Max scanlines max_scanlines = 65535 + lonlat_sample_points = LAC_LONLAT_SAMPLE_POINTS def __init__(self, *args, **kwargs): """Init the LAC reader.""" @@ -53,9 +54,9 @@ def _validate_header(cls, header): # call super to enter the Method Resolution Order (MRO) super(LACReader, cls)._validate_header(header) LOG.debug("validate header") - data_set_name = header['data_set_name'].decode() + data_set_name = header["data_set_name"].decode() # split header into parts creation_site, transfer_mode, platform_id = ( - data_set_name.split('.')[:3]) + data_set_name.split(".")[:3]) if transfer_mode not in ["LHRR", "HRPT", "FRAC"]: raise ReaderError('Improper transfer mode "%s"!' % transfer_mode) diff --git a/pygac/patmosx_coeff_reader.py b/pygac/patmosx_coeff_reader.py index 33c2258..b2cadfe 100644 --- a/pygac/patmosx_coeff_reader.py +++ b/pygac/patmosx_coeff_reader.py @@ -37,35 +37,35 @@ class PatmosxReader: """Read PATMOS-x coefficient files tarballs.""" # regular expression with named capturing groups to read an entire patmosx file regex = re.compile( - r'\s*(?P\w+)[^\n]*\n' - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' + r"\s*(?P\w+)[^\n]*\n" + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" r'\s*(?P[eE0-9\.-]+),?\s*(?P[eE0-9\.-]+),?\s*(?P[eE0-9\.-]+),?\s*(?P[eE0-9\.-]+)[^\n]*\n' # noqa r'\s*(?P[eE0-9\.-]+),?\s*(?P[eE0-9\.-]+),?\s*(?P[eE0-9\.-]+),?\s*(?P[eE0-9\.-]+)[^\n]*\n' # noqa r'\s*(?P[eE0-9\.-]+),?\s*(?P[eE0-9\.-]+),?\s*(?P[eE0-9\.-]+),?\s*(?P[eE0-9\.-]+)[^\n]*\n' # noqa - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' - r'(?:[a-z]+[^\n]*\n)?' + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" + r"(?:[a-z]+[^\n]*\n)?" r'\s*(?P[eE0-9\.-]+)\s*(?P[eE0-9\.-]+)\s*(?P[eE0-9\.-]+)[^\n]*\n' # noqa r'\s*(?P[eE0-9\.-]+)\s*(?P[eE0-9\.-]+)\s*(?P[eE0-9\.-]+)[^\n]*\n' # noqa r'\s*(?P[eE0-9\.-]+)\s*(?P[eE0-9\.-]+)\s*(?P[eE0-9\.-]+)[^\n]*\n' # noqa r'\s*(?P[eE0-9\.-]+)\s*(?P[eE0-9\.-]+)\s*(?P[eE0-9\.-]+)[^\n]*\n' # noqa r'\s*(?P[eE0-9\.-]+)\s*(?P[eE0-9\.-]+)\s*(?P[eE0-9\.-]+)[^\n]*\n' # noqa r'\s*(?P[eE0-9\.-]+)\s*(?P[eE0-9\.-]+)\s*(?P[eE0-9\.-]+)[^\n]*\n' # noqa - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' - r'\s*(?P[eE0-9\.-]+)[^\n]*\n' + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" + r"\s*(?P[eE0-9\.-]+)[^\n]*\n" r'\s*(?P[eE0-9\.-]+)\,*\s*(?P[eE0-9\.-]+)\,*\s*(?P[eE0-9\.-]+)\,*\s*(?P[eE0-9\.-]+)\,*\s*(?P[eE0-9\.-]+)[^\n]*\n' # noqa r'\s*(?P[eE0-9\.-]+)\,*\s*(?P[eE0-9\.-]+)\,*\s*(?P[eE0-9\.-]+)\,*\s*(?P[eE0-9\.-]+)\,*\s*(?P[eE0-9\.-]+)[^\n]*\n' # noqa r'\s*(?P[eE0-9\.-]+)\,*\s*(?P[eE0-9\.-]+)\,*\s*(?P[eE0-9\.-]+)\,*\s*(?P[eE0-9\.-]+)\,*\s*(?P[eE0-9\.-]+)[^\n]*\n' # noqa @@ -73,8 +73,8 @@ class PatmosxReader: r'\s*(?P[eE0-9\.-]+)\s*(?P[eE0-9\.-]+)\s*(?P[eE0-9\.-]+)\s*(?P[eE0-9\.-]+)[^\n]*\n' # noqa r'(?:\s*(?P[eE0-9\.-]+),\s*(?P[eE0-9\.-]+),\s*(?P[eE0-9\.-]+),\s*(?P[eE0-9\.-]+)[^\n]*\n)?' # noqa r'(?:\s*(?P[eE0-9\.-]+),\s*(?P[eE0-9\.-]+),\s*(?P[eE0-9\.-]+),\s*(?P[eE0-9\.-]+)[^\n]*\n)?' # noqa - r'(?:\![^v][^\n]*\n)*' - r'(?:\!(?Pv\w+))?' + r"(?:\![^v][^\n]*\n)*" + r"(?:\!(?Pv\w+))?" ) def __init__(self, tarball): @@ -115,9 +115,9 @@ class Translator: sat_names.update({"n{0:02d}".format(i): "noaa{0}".format(i) for i in range(6,20)}) description = { "visible": { - "channels": ['1', '2', '3a'], + "channels": ["1", "2", "3a"], "coefficients": { - 'dark_count': "instrument counts under dark conditions []", + "dark_count": "instrument counts under dark conditions []", "gain_switch": "dual-gain switch count, set to 'null' for single-gain instruments []", "s0": "single-gain calibration slope at launch date [%]", "s1": "linear single-gain calibration slope parameter [% years^{-1}]", @@ -127,9 +127,9 @@ class Translator: "method": 'Heidinger, A.K., W.C. Straka III, C.C. Molling, J.T. Sullivan, and X. Wu, 2010: Deriving an inter-sensor consistent calibration for the AVHRR solar reflectance data record. International Journal of Remote Sensing, 31:6493-6517' # noqa }, "thermal": { - "channels": ['3b', '4', '5'], + "channels": ["3b", "4", "5"], "coefficients": { - 'centroid_wavenumber': "centroid wavenumber [cm^{-1}]", + "centroid_wavenumber": "centroid wavenumber [cm^{-1}]", "b0": "constant non-linear radiance correction coefficient [mW m^{-2} sr cm^{-1}]", "b1": "linear non-linear radiance correction coefficient []", "b2": "quadratic non-linear radiance correction coefficient [(mW^{-1} m^2 sr^{-1} cm)]", @@ -158,34 +158,34 @@ def convert(cls, patmosx_sat_coeffs): pygac_sat_coeffs = {} # visible calibration for ch in ("1", "2", "3a"): - s0l = patmosx_sat_coeffs['ch{0}_low_gain_S0'.format(ch)] - s0h = patmosx_sat_coeffs['ch{0}_high_gain_S0'.format(ch)] + s0l = patmosx_sat_coeffs["ch{0}_low_gain_S0".format(ch)] + s0h = patmosx_sat_coeffs["ch{0}_high_gain_S0".format(ch)] if s0l == s0h: gain_switch = None s0 = s0l else: - gain_switch = patmosx_sat_coeffs['ch{0}_gain_switches_count'.format(ch)] + gain_switch = patmosx_sat_coeffs["ch{0}_gain_switches_count".format(ch)] s0 = cls.find_s0(s0l, s0h, ch) pygac_sat_coeffs["channel_{0}".format(ch)] = { - "dark_count": float(patmosx_sat_coeffs['ch{0}_dark_count'.format(ch)]), + "dark_count": float(patmosx_sat_coeffs["ch{0}_dark_count".format(ch)]), "gain_switch": gain_switch, "s0": s0, - "s1": patmosx_sat_coeffs['ch{0}_high_gain_S1'.format(ch)], - "s2": patmosx_sat_coeffs['ch{0}_high_gain_S2'.format(ch)] + "s1": patmosx_sat_coeffs["ch{0}_high_gain_S1".format(ch)], + "s2": patmosx_sat_coeffs["ch{0}_high_gain_S2".format(ch)] } - date_of_launch = cls.float2date(patmosx_sat_coeffs['date_of_launch']) - pygac_sat_coeffs['date_of_launch'] = date_of_launch.strftime("%Y-%m-%dT%H:%M:%S.%fZ") + date_of_launch = cls.float2date(patmosx_sat_coeffs["date_of_launch"]) + pygac_sat_coeffs["date_of_launch"] = date_of_launch.strftime("%Y-%m-%dT%H:%M:%S.%fZ") # thermal channels for ch in ("3b", "4", "5"): pygac_sat_coeffs["channel_{0}".format(ch)] = { - "b0": patmosx_sat_coeffs['ch{0}_b0'.format(ch)], - "b1": patmosx_sat_coeffs['ch{0}_b1'.format(ch)], - "b2": patmosx_sat_coeffs['ch{0}_b2'.format(ch)], - "centroid_wavenumber": patmosx_sat_coeffs['nu_{0}'.format(ch)], - "space_radiance": patmosx_sat_coeffs['ch{0}_Ns'.format(ch)], - "to_eff_blackbody_intercept": (-patmosx_sat_coeffs['a1_{0}'.format(ch)] - / patmosx_sat_coeffs['a2_{0}'.format(ch)]), - "to_eff_blackbody_slope": 1/patmosx_sat_coeffs['a2_{0}'.format(ch)] + "b0": patmosx_sat_coeffs["ch{0}_b0".format(ch)], + "b1": patmosx_sat_coeffs["ch{0}_b1".format(ch)], + "b2": patmosx_sat_coeffs["ch{0}_b2".format(ch)], + "centroid_wavenumber": patmosx_sat_coeffs["nu_{0}".format(ch)], + "space_radiance": patmosx_sat_coeffs["ch{0}_Ns".format(ch)], + "to_eff_blackbody_intercept": (-patmosx_sat_coeffs["a1_{0}".format(ch)] + / patmosx_sat_coeffs["a2_{0}".format(ch)]), + "to_eff_blackbody_slope": 1/patmosx_sat_coeffs["a2_{0}".format(ch)] } for t in range(1, 5): pygac_sat_coeffs["thermometer_{0}".format(t)] = { @@ -209,7 +209,7 @@ def find_s0(s0_low, s0_high, ch): if s0_low == s0_high: # single gain case return s0_low - if ch == '3a': + if ch == "3a": g_low, g_high = 0.25, 1.75 else: g_low, g_high = 0.5, 1.5 @@ -249,32 +249,32 @@ def float2date(date_float): def save(self, filepath): """Save coefficients as PyGAC json file.""" - with open(filepath, mode='w') as json_file: + with open(filepath, mode="w") as json_file: json.dump(self.coeffs, json_file, indent=4, sort_keys=True) def main(): """The main function.""" parser = argparse.ArgumentParser(description=__doc__) - parser.add_argument('tarball', type=str, help='path to PATMOS-x coefficients tarball') - parser.add_argument('-o', '--output', type=str, metavar="JSON", + parser.add_argument("tarball", type=str, help="path to PATMOS-x coefficients tarball") + parser.add_argument("-o", "--output", type=str, metavar="JSON", help='path to PyGAC json file, defaults to tarball path with suffix ".json"') - parser.add_argument('-v', '--verbose', action='store_true', help='explain what is being done') + parser.add_argument("-v", "--verbose", action="store_true", help="explain what is being done") args = parser.parse_args() if args.verbose: loglevel = logging.INFO else: loglevel = logging.WARNING - logging.basicConfig(level=loglevel, format='[%(asctime)s] %(message)s') + logging.basicConfig(level=loglevel, format="[%(asctime)s] %(message)s") tarball = pathlib.Path(args.tarball) logging.info('Read PATMOS-x tarball "%s".', tarball) patmosx_coeffs = PatmosxReader(tarball) - logging.info('Translate PATMOS-x coefficients to PyGAC format.') + logging.info("Translate PATMOS-x coefficients to PyGAC format.") pygac_coeffs = Translator(patmosx_coeffs) output = args.output or tarball.with_suffix(".json") logging.info('Write PyGAC calibration json file "%s".', output) pygac_coeffs.save(output) - logging.info('Done!') + logging.info("Done!") if __name__ == "__main__": main() diff --git a/pygac/pod_reader.py b/pygac/pod_reader.py index 39fc900..0990e09 100644 --- a/pygac/pod_reader.py +++ b/pygac/pod_reader.py @@ -198,44 +198,44 @@ class POD_QualityIndicator(IntFlag): ("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')]) +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', + 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', + 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 @@ -246,7 +246,7 @@ class PODReader(Reader): def correct_scan_line_numbers(self): """Correct the scan line numbers.""" # Perform common corrections first. - super(PODReader, self).correct_scan_line_numbers() + super().correct_scan_line_numbers() # cleaning up the data min_scanline_number = np.amin( @@ -275,7 +275,7 @@ def read(self, filename, fileobj=None): """ self.filename = filename - LOG.info('Reading %s', self.filename) + LOG.info("Reading %s", self.filename) # choose the right header depending on the date with file_opener(fileobj or filename) as fd_: self.tbm_head, self.head = self.read_header( @@ -341,8 +341,8 @@ def read_header(cls, filename, fileobj=None, header_date="auto"): @classmethod def _validate_tbm_header(cls, potential_tbm_header): - data_set_name = potential_tbm_header['data_set_name'] - allowed_empty = (42*b'\x00' + b' ') + data_set_name = potential_tbm_header["data_set_name"] + allowed_empty = (42*b"\x00" + b" ") if data_set_name == allowed_empty: return potential_tbm_header.copy() @@ -382,12 +382,12 @@ def _validate_header(cls, header): # call super to enter the Method Resolution Order (MRO) super(PODReader, cls)._validate_header(header) LOG.debug("validate header") - data_set_name = header['data_set_name'].decode() + data_set_name = header["data_set_name"].decode() # split header into parts creation_site, transfer_mode, platform_id = ( - data_set_name.split('.')[:3]) - allowed_ids = ['TN', 'NA', 'NB', 'NC', 'ND', 'NE', 'NF', 'NG', - 'NH', 'NI', 'NJ'] + data_set_name.split(".")[:3]) + allowed_ids = ["TN", "NA", "NB", "NC", "ND", "NE", "NF", "NG", + "NH", "NI", "NJ"] if platform_id not in allowed_ids: raise ReaderError('Improper platform id "%s"!' % platform_id) @@ -406,7 +406,7 @@ def get_header_timestamp(self): 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)) + raise ValueError("Corrupt header timestamp: {0}".format(err)) @staticmethod def decode_timestamps(encoded): @@ -430,7 +430,7 @@ def decode_timestamps(encoded): enc1 = encoded[:, 1] enc2 = encoded[:, 2] else: - raise ValueError('Invalid timestamp dimension') + raise ValueError("Invalid timestamp dimension") year = enc0 >> 9 year = np.where(year > 75, year + 1900, year + 2000) @@ -439,16 +439,15 @@ def decode_timestamps(encoded): return year, jday, msec - def _get_times(self): + def _get_times_from_file(self): return self.decode_timestamps(self.scans["time_code"]) def _compute_missing_lonlat(self, missed_utcs): """Compute lon lat values using pyorbital.""" tic = datetime.datetime.now() - scan_rate = datetime.timedelta(milliseconds=1/self.scan_freq).total_seconds() sgeom = avhrr_gac(missed_utcs.astype(datetime.datetime), - self.scan_points, frequency=scan_rate) + self.lonlat_sample_points, frequency=scan_rate) t0 = missed_utcs[0].astype(datetime.datetime) s_times = sgeom.times(t0) tle1, tle2 = self.get_tle_lines() @@ -479,10 +478,10 @@ def _adjust_clock_drift(self): self.spacecraft_name) return - error_utcs = np.array(error_utcs, dtype='datetime64[ms]') + error_utcs = np.array(error_utcs, dtype="datetime64[ms]") # interpolate to get the clock offsets at the scan line utcs # the clock_error is given in seconds, so offsets are in seconds, too. - offsets = np.interp(self.utcs.astype(np.uint64), + offsets = np.interp(self._times_as_np_datetime64.astype(np.uint64), error_utcs.astype(np.uint64), clock_error) LOG.info("Adjusting for clock drift of %s to %s seconds.", @@ -507,12 +506,12 @@ def _adjust_clock_drift(self): num_lines = max_line - min_line + 1 missed_lines = np.setdiff1d(np.arange(min_line, max_line+1), scan_lines) missed_utcs = ((missed_lines - scan_lines[0])*np.timedelta64(scan_rate, "ms") - + self.utcs[0]) + + self._times_as_np_datetime64[0]) # calculate the missing geo locations try: missed_lons, missed_lats = self._compute_missing_lonlat(missed_utcs) except NoTLEData as err: - LOG.warning('Cannot perform clock drift correction: %s', str(err)) + LOG.warning("Cannot perform clock drift correction: %s", str(err)) return # create arrays of lons and lats for interpolation. The locations @@ -540,14 +539,14 @@ def _adjust_clock_drift(self): # set corrected values self.lons = slerp_res[:, :, 0] self.lats = slerp_res[:, :, 1] - self.utcs -= (offsets * 1000).astype('timedelta64[ms]') + self._times_as_np_datetime64 -= (offsets * 1000).astype("timedelta64[ms]") 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 + def _get_lonlat_from_file(self): + lats = self.scans["earth_location"]["lats"] / np.float32(128.0) + lons = self.scans["earth_location"]["lons"] / np.float32(128.0) return lons, lats def get_telemetry(self): diff --git a/pygac/pygac_geotiepoints.py b/pygac/pygac_geotiepoints.py index 96e5ffb..ed41c60 100644 --- a/pygac/pygac_geotiepoints.py +++ b/pygac/pygac_geotiepoints.py @@ -25,9 +25,11 @@ import geotiepoints as gtp - import numpy as np +GAC_LONLAT_SAMPLE_POINTS = np.arange(4, 405, 8) +LAC_LONLAT_SAMPLE_POINTS = np.arange(24, 2048, 40) + def gac_lat_lon_interpolator(lons_subset, lats_subset): """Interpolate lat-lon values in the AVHRR GAC data. @@ -37,16 +39,15 @@ def gac_lat_lon_interpolator(lons_subset, lats_subset): ranging from 5 to 405. Interpolate to full resolution. """ - cols_subset = np.arange(4, 405, 8) cols_full = np.arange(409) - return lat_lon_interpolator(lons_subset, lats_subset, cols_subset, cols_full) + return lat_lon_interpolator(lons_subset, lats_subset, GAC_LONLAT_SAMPLE_POINTS, 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) + return lat_lon_interpolator(lons_subset, lats_subset, LAC_LONLAT_SAMPLE_POINTS, cols_full) def lat_lon_interpolator(lons_subset, lats_subset, cols_subset, cols_full): diff --git a/pygac/reader.py b/pygac/reader.py index 88e4c87..1883727 100644 --- a/pygac/reader.py +++ b/pygac/reader.py @@ -30,17 +30,18 @@ import re import types import warnings -from abc import ABCMeta, abstractmethod +from abc import ABC, abstractmethod +from contextlib import suppress +from importlib.metadata import entry_points import numpy as np import pyorbital -import six +import xarray as xr from packaging.version import Version from pyorbital import astronomy from pyorbital.orbital import Orbital from pygac import gac_io -from pygac.calibration import Calibrator, calibrate_solar, calibrate_thermal from pygac.utils import calculate_sun_earth_distance_correction, centered_modulus, get_absolute_azimuth_angle_diff LOG = logging.getLogger(__name__) @@ -48,29 +49,29 @@ # 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, + "noaa7": {"roll": 0.000, + "pitch": 0.000, + "yaw": 0.000, }, - 'noaa9': {'roll': 0.000, - 'pitch': 0.0025, - 'yaw': 0.000, + "noaa9": {"roll": 0.000, + "pitch": 0.0025, + "yaw": 0.000, }, - 'noaa10': {'roll': 0.000, - 'pitch': 0.000, - 'yaw': 0.000, + "noaa10": {"roll": 0.000, + "pitch": 0.000, + "yaw": 0.000, }, - 'noaa11': {'roll': -0.0019, - 'pitch': -0.0037, - 'yaw': 0.000, + "noaa11": {"roll": -0.0019, + "pitch": -0.0037, + "yaw": 0.000, }, - 'noaa12': {'roll': 0.000, - 'pitch': 0.000, - 'yaw': 0.000, + "noaa12": {"roll": 0.000, + "pitch": 0.000, + "yaw": 0.000, }, - 'noaa14': {'roll': 0.000, - 'pitch': 0.000, - 'yaw': 0.000, + "noaa14": {"roll": 0.000, + "pitch": 0.000, + "yaw": 0.000, }} @@ -88,16 +89,21 @@ class DecodingError(ValueError): """Raised when decoding of some value fails.""" -class Reader(six.with_metaclass(ABCMeta)): +class TimestampMismatch(IndexError): + """Raise when matching timestamps doesn't work.""" + + +class Reader(ABC): """Reader for GAC and LAC format, POD and KLM platforms.""" # data set header format, see _validate_header for more details data_set_pattern = re.compile( - r'\w{3}\.\w{4}\.\w{2}.D\d{5}\.S\d{4}\.E\d{4}\.B\d{7}\.\w{2}') + r"\w{3}\.\w{4}\.\w{2}.D\d{5}\.S\d{4}\.E\d{4}\.B\d{7}\.\w{2}") def __init__(self, interpolate_coords=True, adjust_clock_drift=True, tle_dir=None, tle_name=None, tle_thresh=7, creation_site=None, - custom_calibration=None, calibration_file=None, header_date="auto"): + custom_calibration=None, calibration_file=None, header_date="auto", + calibration_method="noaa", calibration_parameters=None): """Init the reader. Args: @@ -122,15 +128,13 @@ def __init__(self, interpolate_coords=True, adjust_clock_drift=True, 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.custom_calibration = custom_calibration - self.calibration_file = calibration_file + self.creation_site = (creation_site or "NSS").encode("utf-8") self.header_date = header_date self.head = None self.scans = None self.spacecraft_name = None self.spacecraft_id = None - self.utcs = None + self._times_as_np_datetime64 = None self.lats = None self.lons = None self.tle_lines = None @@ -138,10 +142,22 @@ def __init__(self, interpolate_coords=True, adjust_clock_drift=True, self._mask = None self._rpy = None + if calibration_method.lower() not in ["noaa"]: + raise ValueError(f"Unknown calibration method {calibration_method}.") + self.calibration_method = calibration_method + self.calibration_parameters = calibration_parameters or dict() + if custom_calibration is not None or calibration_file is not None: + warnings.warn("`custom_calibration` and `calibration_file` parameters will be deprecated soon. " + "Please use `calibration_parameters` dictionary instead.", + PendingDeprecationWarning) + if not self.calibration_parameters: + self.calibration_parameters = dict(custom_coeffs=custom_calibration, + coeffs_file=calibration_file) + @property def times(self): """Get the UTCs as datetime.datetime.""" - return self.to_datetime(self.utcs) + return self.to_datetime(self._times_as_np_datetime64) @property def filename(self): @@ -161,16 +177,6 @@ def filename(self, filepath): else: self._filename = os.path.basename(filepath) - @property - def calibration(self): - """Get the property 'calibration'.""" - calibration = Calibrator( - self.spacecraft_name, - custom_coeffs=self.custom_calibration, - coeffs_file=self.calibration_file - ) - return calibration - @abstractmethod def read(self, filename, fileobj=None): # pragma: no cover """Read the GAC/LAC data. @@ -209,7 +215,7 @@ def _correct_data_set_name(cls, header, filename): filename (str): path to file """ filename = str(filename) - data_set_name = header['data_set_name'] + data_set_name = header["data_set_name"] try: header["data_set_name"] = cls._decode_data_set_name(data_set_name) except DecodingError: @@ -219,10 +225,10 @@ def _correct_data_set_name(cls, header, filename): if match: data_set_name = match.group() LOG.debug(f"Set data_set_name, to filename {data_set_name}") - header['data_set_name'] = data_set_name.encode() + header["data_set_name"] = data_set_name.encode() else: LOG.debug(f"header['data_set_name']={header['data_set_name']}; filename='{filename}'") - raise ReaderError('Cannot determine data_set_name!') + raise ReaderError("Cannot determine data_set_name!") return header @classmethod @@ -239,10 +245,10 @@ def _decode_data_set_name(cls, data_set_name): @classmethod def _decode_data_set_name_for_encoding(cls, data_set_name, encoding): - data_set_name = data_set_name.decode(encoding, errors='ignore') + data_set_name = data_set_name.decode(encoding, errors="ignore") if not cls.data_set_pattern.match(data_set_name): - raise DecodingError(f'The data_set_name in header {data_set_name} ' - f'does not seem correct using encoding {encoding}.') + raise DecodingError(f"The data_set_name in header {data_set_name} " + f"does not seem correct using encoding {encoding}.") else: data_set_name = data_set_name.encode() return data_set_name @@ -280,10 +286,10 @@ def _validate_header(cls, header): # second use case "diamond diagrams". # Check if the data set name matches the pattern LOG.debug("validate header") - data_set_name = header['data_set_name'].decode(errors='ignore') + data_set_name = header["data_set_name"].decode(errors="ignore") if not cls.data_set_pattern.match(data_set_name): - raise ReaderError('Data set name %s does not match!' - % header['data_set_name']) + raise ReaderError("Data set name %s does not match!" + % header["data_set_name"]) def _read_scanlines(self, buffer, count): """Read the scanlines from the given buffer. @@ -371,7 +377,7 @@ def save(self, start_line, end_line, output_file_prefix="PyGAC", output_dir="./" raise ValueError("All data is masked out.") gac_io.save_gac( self.spacecraft_name, - self.utcs, self.lats, self.lons, + self._times_as_np_datetime64, self.lats, self.lons, channels[:, :, 0], channels[:, :, 1], channels[:, :, 2], channels[:, :, 3], channels[:, :, 4], channels[:, :, 5], @@ -428,8 +434,8 @@ def get_counts(self): return channels @abstractmethod - def _get_times(self): # pragma: no cover - """Specify how to read scanline timestamps from GAC data. + def _get_times_from_file(self): # pragma: no cover + """Specify how to read scanline timestamps from data. Returns: int: year @@ -442,24 +448,24 @@ def _get_times(self): # pragma: no cover 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: + if self._times_as_np_datetime64 is None: # Read timestamps - year, jday, msec = self._get_times() + year, jday, msec = self._get_times_from_file() # 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() + self._times_as_np_datetime64 = self.to_datetime64(year=year, jday=jday, msec=msec) + try: + self._times_as_np_datetime64 = self.correct_times_thresh() + except TimestampMismatch as err: + LOG.error(str(err)) - return self.utcs + return self._times_as_np_datetime64 @staticmethod def to_datetime64(year, jday, msec): @@ -474,9 +480,9 @@ def to_datetime64(year, jday, msec): numpy.datetime64: Converted timestamps """ - return (year.astype(str).astype('datetime64[Y]') - + (jday - 1).astype('timedelta64[D]') - + msec.astype('timedelta64[ms]')) + return (year.astype(str).astype("datetime64[Y]") + + (jday - 1).astype("timedelta64[D]") + + msec.astype("timedelta64[ms]")) @staticmethod def to_datetime(datetime64): @@ -509,72 +515,138 @@ def lineno2msec(self, scan_line_number): def get_sun_earth_distance_correction(self): """Get the julian day and the sun-earth distance correction.""" self.get_times() - jday = self.times[0].timetuple().tm_yday + # Shouldn't we use the start time from the header if available? + start_time = self._times_as_np_datetime64[0] + jday = (start_time - start_time.astype("datetime64[Y]")).astype("timedelta64[D]").astype(int) + 1 return calculate_sun_earth_distance_correction(jday) def update_meta_data(self): - """Add some metd data to the meta_data dicitonary.""" - if 'sun_earth_distance_correction_factor' not in self.meta_data: - self.meta_data['sun_earth_distance_correction_factor'] = ( + """Add some meta data to the meta_data dicitonary.""" + meta_data = self.meta_data + self._update_meta_data_object(meta_data) + if "gac_header" not in meta_data: + meta_data["gac_header"] = self.head.copy() + + def _update_meta_data_object(self, meta_data): + if "sun_earth_distance_correction_factor" not in meta_data: + meta_data["sun_earth_distance_correction_factor"] = ( self.get_sun_earth_distance_correction()) - if 'midnight_scanline' not in self.meta_data: - self.meta_data['midnight_scanline'] = self.get_midnight_scanline() - if 'missing_scanlines' not in self.meta_data: - self.meta_data['missing_scanlines'] = self.get_miss_lines() - if 'gac_header' not in self.meta_data: - self.meta_data['gac_header'] = self.head.copy() - self.meta_data['calib_coeffs_version'] = self.calibration.version + if "midnight_scanline" not in meta_data: + meta_data["midnight_scanline"] = self.get_midnight_scanline() + if "missing_scanlines" not in meta_data: + meta_data["missing_scanlines"] = self.get_miss_lines() + + def read_as_dataset(self, file_to_read): + self.read(file_to_read) + return self.create_counts_dataset() + + def create_counts_dataset(self): + head = dict(zip(self.head.dtype.names, self.head.item())) + scans = self.scans + + times = self.get_times() + line_numbers = scans["scan_line_number"] + counts = self.get_counts() + + scan_size = counts.shape[1] + + if counts.shape[-1] == 5: + channel_names = ["1", "2", "3", "4", "5"] + ir_channel_names = ["3", "4", "5"] + else: + channel_names = ["1", "2", "3a", "3b", "4", "5"] + ir_channel_names = ["3b", "4", "5"] + + columns = np.arange(scan_size) + channels = xr.DataArray(counts, dims=["scan_line_index", "columns", "channel_name"], + coords=dict(scan_line_index=line_numbers, + channel_name=channel_names, + columns=columns, + times=("scan_line_index", times))) + + prt, ict, space = self._get_telemetry_dataarrays(line_numbers, ir_channel_names) + + longitudes, latitudes = self._get_lonlat_dataarrays(line_numbers, columns) + + if self.interpolate_coords: + channels = channels.assign_coords(longitude=(("scan_line_index", "columns"), + longitudes.reindex_like(channels).data), + latitude=(("scan_line_index", "columns"), + latitudes.reindex_like(channels).data)) + + ds = xr.Dataset(dict(channels=channels, prt_counts=prt, ict_counts=ict, space_counts=space, + longitude=longitudes, latitude=latitudes), + attrs=head) + + ds.attrs["spacecraft_name"] = self.spacecraft_name + self._update_meta_data_object(ds.attrs) + return ds + + def _get_lonlat_dataarrays(self, line_numbers, columns): + lons, lats = self.get_lonlat() + lon_lat_columns = columns if self.interpolate_coords is True else self.lonlat_sample_points + + if self.interpolate_coords: + longitudes = xr.DataArray(lons, dims=["scan_line_index", "columns"], + coords={"scan_line_index": line_numbers, + "columns": lon_lat_columns}) + latitudes = xr.DataArray(lats, dims=["scan_line_index", "columns"], + coords={"scan_line_index": line_numbers, + "columns": lon_lat_columns}) + + else: + longitudes = xr.DataArray(lons, dims=["scan_line_index", "subsampled_columns"], + coords={"scan_line_index": line_numbers, + "subsampled_columns": lon_lat_columns}) + latitudes = xr.DataArray(lats, dims=["scan_line_index", "subsampled_columns"], + coords={"scan_line_index": line_numbers, + "subsampled_columns": lon_lat_columns}) + + return longitudes,latitudes + + def _get_telemetry_dataarrays(self, line_numbers, ir_channel_names): + prt, ict, space = self.get_telemetry() + + prt = xr.DataArray(prt, dims=["scan_line_index"], coords=dict(scan_line_index=line_numbers)) + ict = xr.DataArray(ict, dims=["scan_line_index", "ir_channel_name"], + coords=dict(ir_channel_name=ir_channel_names, scan_line_index=line_numbers)) + space = xr.DataArray(space, dims=["scan_line_index", "ir_channel_name"], + coords=dict(ir_channel_name=ir_channel_names, scan_line_index=line_numbers)) + + return prt,ict,space def get_calibrated_channels(self): """Calibrate and return the channels.""" - channels = self.get_counts() - self.get_times() - times = self.times - self.update_meta_data() - year = times[0].year - delta = times[0].date() - datetime.date(year, 1, 1) - jday = delta.days + 1 - - corr = self.meta_data['sun_earth_distance_correction_factor'] - calibration_coeffs = self.calibration - - # 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, - calibration_coeffs, - corr - ) - prt, ict, space = self.get_telemetry() + calibrated_ds = self.get_calibrated_dataset() - ir_channels_to_calibrate = self._get_ir_channels_to_calibrate() - - for chan in ir_channels_to_calibrate: - channels[:, :, chan - 6] = calibrate_thermal( - channels[:, :, chan - 6], - prt, - ict[:, chan - 3], - space[:, chan - 3], - self.scans["scan_line_number"], - chan, - calibration_coeffs - ) + channels = calibrated_ds["channels"].data + with suppress(KeyError): + self.meta_data["calib_coeffs_version"] = calibrated_ds.attrs["calib_coeffs_version"] + + return channels + + def get_calibrated_dataset(self): + """Create and calibrate the dataset for the pass.""" + ds = self.create_counts_dataset() + # calibration = {"1": "mitram", "2": "mitram", "4": {"method":"noaa", "coeff_file": "myfile.json"}} + + calibration_entrypoints = entry_points(group="pygac.calibration") + calibration_function = calibration_entrypoints[self.calibration_method].load() + calibrated_ds = calibration_function(ds, **self.calibration_parameters) # Mask out corrupt values - channels[self.mask] = np.nan + mask = xr.DataArray(self.mask==False, dims=["scan_line_index"]) # noqa + calibrated_ds = calibrated_ds.where(mask) + # Apply KLM/POD specific postprocessing - self.postproc(channels) + self.postproc(calibrated_ds) # 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 + LOG.info("Correcting for temporary scan motor issue") + self.mask_tsm_pixels(calibrated_ds) + return calibrated_ds @abstractmethod def get_telemetry(self): # pragma: no cover @@ -587,18 +659,18 @@ def get_lonlat(self): TODO: Switch to faster interpolator? """ if self.lons is None and self.lats is None: - self.lons, self.lats = self._get_lonlat() + self.lons, self.lats = self._get_lonlat_from_file() self.update_meta_data() + # Adjust clock drift + if self.adjust_clock_drift: + self._adjust_clock_drift() + # 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 @@ -610,7 +682,7 @@ def get_lonlat(self): return self.lons, self.lats @abstractmethod - def _get_lonlat(self): # pragma: no cover + def _get_lonlat_from_file(self): # pragma: no cover """KLM/POD specific readout of lat/lon coordinates.""" raise NotImplementedError @@ -662,7 +734,7 @@ def get_qual_flags(self): return qual_flags @abstractmethod - def postproc(self, channels): # pragma: no cover + def postproc(self, ds): # pragma: no cover """Apply KLM/POD specific postprocessing.""" raise NotImplementedError @@ -682,13 +754,13 @@ def tle2datetime64(times): times = np.where(times > 50000, times + 1900000, times + 2000000) # Convert float to datetime64 - doys = (times % 1000).astype('int') - 1 - years = (times // 1000).astype('int') + 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]') + years - 1970).astype("datetime64[Y]").astype("datetime64[ms]") + times64 += doys.astype("timedelta64[D]") + times64 += msecs.astype("timedelta64[ms]") return times64 @@ -701,13 +773,13 @@ def get_tle_file(self): raise RuntimeError("TLE name not specified!") values = {"satname": self.spacecraft_name, } tle_filename = os.path.join(tle_dir, tle_name % values) - LOG.info('TLE filename = ' + str(tle_filename)) + 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_: + with open(tle_filename, "r") as fp_: return fp_.readlines() def get_tle_lines(self): @@ -721,7 +793,7 @@ def get_tle_lines(self): return self.tle_lines self.get_times() tle_data = self.read_tle_file(self.get_tle_file()) - sdate = self.utcs[0] + sdate = self._times_as_np_datetime64[0] dates = self.tle2datetime64( np.array([float(line[18:32]) for line in tle_data[::2]])) @@ -740,7 +812,7 @@ def get_tle_lines(self): iindex -= 1 # Make sure the TLE we found is within the threshold - delta_days = abs(sdate - dates[iindex]) / np.timedelta64(1, 'D') + delta_days = abs(sdate - dates[iindex]) / np.timedelta64(1, "D") if delta_days > self.tle_thresh: raise NoTLEData( "Can't find tle data for %s within +/- %d days around %s" % @@ -769,8 +841,8 @@ def get_sat_angles(self): return self._get_sat_angles_with_tle() except NoTLEData: LOG.warning( - 'No TLE data available. Falling back to approximate ' - 'calculation of satellite angles.' + "No TLE data available. Falling back to approximate " + "calculation of satellite angles." ) return self._get_sat_angles_without_tle() @@ -778,25 +850,25 @@ def _get_sat_angles_with_tle(self): 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], + sat_azi, sat_elev = orb.get_observer_look(self._times_as_np_datetime64[:, np.newaxis], self.lons, self.lats, 0) return sat_azi, sat_elev def _get_sat_angles_without_tle(self): """Get satellite angles using lat/lon from data to approximate satellite postition instead of TLE.""" from pyorbital.orbital import get_observer_look as get_observer_look_no_tle - LOG.warning('Approximating satellite height to 850km (TIROS-N OSCAR)!') + LOG.warning("Approximating satellite height to 850km (TIROS-N OSCAR)!") sat_alt = 850.0 # km TIROS-N OSCAR mid_column = int(0.5*self.lons.shape[1]) sat_azi, sat_elev = get_observer_look_no_tle( self.lons[:, mid_column][:, np.newaxis], self.lats[:, mid_column][:, np.newaxis], # approximate satellite position sat_alt, # approximate satellite altitude - self.times[:, np.newaxis], + self._times_as_np_datetime64[:, np.newaxis], self.lons, self.lats, 0) # Sometimes (pyorbital <= 1.6.1) the get_observer_look_not_tle returns nodata instead of 90. # Problem solved with https://github.com/pytroll/pyorbital/pull/77 - if Version(pyorbital.__version__) <= Version('1.6.1'): + if Version(pyorbital.__version__) <= Version("1.6.1"): sat_elev[:, mid_column] = 90 return sat_azi, sat_elev @@ -824,7 +896,7 @@ def get_angles(self): """ self.get_times() self.get_lonlat() - times = self.times + times = self._times_as_np_datetime64 sat_azi, sat_elev = self.get_sat_angles() sat_zenith = 90 - sat_elev @@ -922,8 +994,8 @@ def correct_scan_line_numbers(self): """ 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()} + 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"] < self.max_scanlines, @@ -966,14 +1038,14 @@ def correct_scan_line_numbers(self): 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', + 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}) + 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, @@ -1024,15 +1096,15 @@ def correct_times_thresh(self, max_diff_from_t0_head=6*60*1000, # Check whether scanline number increases monotonically nums = self.scans["scan_line_number"] - results.update({'t': self.utcs.copy(), 'n': nums}) + results.update({"t": self._times_as_np_datetime64.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" + results["fail_reason"] = "Scanline number jumps backwards" return results # Convert time to milliseconds since 1970-01-01 - t = self.utcs.astype("i8") + t = self._times_as_np_datetime64.astype("i8") try: t0_head = np.array([self.get_header_timestamp().isoformat()], dtype="datetime64[ms]").astype("i8")[0] @@ -1059,15 +1131,14 @@ def correct_times_thresh(self, max_diff_from_t0_head=6*60*1000, # 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}) + 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 + results["fail_reason"] = "Timestamp mismatch" + raise TimestampMismatch("Timestamp mismatch. Cannot perform correction.") # Add estimated offset to the ideal timestamps tn += t0 @@ -1075,11 +1146,10 @@ def correct_times_thresh(self, max_diff_from_t0_head=6*60*1000, # 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) + self._times_as_np_datetime64[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 + return self._times_as_np_datetime64 @property @abstractmethod @@ -1102,7 +1172,7 @@ def is_tsm_affected(self): """ self.get_times() - ts, te = self.to_datetime(self.utcs[[0, -1]]) + ts, te = self.to_datetime(self._times_as_np_datetime64[[0, -1]]) try: for interval in self.tsm_affected_intervals[self.spacecraft_id]: if ts >= interval[0] and te <= interval[1]: @@ -1124,13 +1194,13 @@ def get_midnight_scanline(self): """ self.get_times() - d0 = np.datetime64(datetime.date(1970, 1, 1), 'D') - days = (self.utcs.astype('datetime64[D]') - d0).astype(int) + d0 = np.datetime64(datetime.date(1970, 1, 1), "D") + days = (self._times_as_np_datetime64.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. ') + LOG.warning("Unable to determine midnight scanline: " + "UTC date increases more than once. ") return None else: return incr[0] @@ -1146,14 +1216,15 @@ def get_miss_lines(self): """ # 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']))) + 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): + def mask_tsm_pixels(self, ds): """Mask pixels affected by the scan motor issue.""" - idx = self.get_tsm_pixels(channels) - channels[idx] = np.nan + idx = self.get_tsm_pixels(ds["channels"].values) + # This is because fancy/pixel indexing doesn't seem to work with xarray's `where` or `loc` + ds["channels"].data[idx] = np.nan @abstractmethod def get_tsm_pixels(self, channels): # pragma: no cover @@ -1196,7 +1267,7 @@ def inherit_doc(cls): 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): + if parfunc and getattr(parfunc, "__doc__", None): func.__doc__ = parfunc.__doc__ break return cls diff --git a/pygac/runner.py b/pygac/runner.py index 7fd7214..1029810 100644 --- a/pygac/runner.py +++ b/pygac/runner.py @@ -31,13 +31,12 @@ class for a given file. import logging import os +from pygac.configuration import get_config from pygac.gac_klm import GACKLMReader from pygac.gac_pod import GACPODReader from pygac.lac_klm import LACKLMReader from pygac.lac_pod import LACPODReader from pygac.utils import file_opener -from pygac.configuration import get_config - LOG = logging.getLogger(__name__) @@ -86,15 +85,15 @@ def process_file(filename, start_line, end_line, fileobj=None): # reader specific values config = get_config() - tle_dir = config.get('tle', 'tledir', raw=True) - tle_name = config.get('tle', 'tlename', raw=True) - coeffs_file = config.get("calibration", "coeffs_file", fallback='') + tle_dir = config.get("tle", "tledir", raw=True) + tle_name = config.get("tle", "tlename", raw=True) + coeffs_file = config.get("calibration", "coeffs_file", fallback="") # output specific values - output_dir = config.get('output', 'output_dir', raw=True) - output_file_prefix = config.get('output', 'output_file_prefix', raw=True) - avhrr_dir = os.environ.get('SM_AVHRR_DIR') - qual_dir = os.environ.get('SM_AVHRR_DIR') - sunsatangles_dir = os.environ.get('SM_SUNSATANGLES_DIR') + output_dir = config.get("output", "output_dir", raw=True) + output_file_prefix = config.get("output", "output_file_prefix", raw=True) + avhrr_dir = os.environ.get("SM_AVHRR_DIR") + qual_dir = os.environ.get("SM_AVHRR_DIR") + sunsatangles_dir = os.environ.get("SM_SUNSATANGLES_DIR") # Keep the file open while searching for the reader class and later # creation of the instance. diff --git a/pygac/tests/test_angles.py b/pygac/tests/test_angles.py index 97b5e6b..e3be986 100644 --- a/pygac/tests/test_angles.py +++ b/pygac/tests/test_angles.py @@ -28,8 +28,7 @@ import numpy as np -from pygac.utils import (get_absolute_azimuth_angle_diff, - centered_modulus) +from pygac.utils import centered_modulus, get_absolute_azimuth_angle_diff class TestAngles(unittest.TestCase): diff --git a/pygac/tests/test_calibrate_klm.py b/pygac/tests/test_calibrate_klm.py index f94ad55..70d595e 100644 --- a/pygac/tests/test_calibrate_klm.py +++ b/pygac/tests/test_calibrate_klm.py @@ -25,13 +25,10 @@ import unittest -try: - import mock -except ImportError: - from unittest import mock + import numpy as np -from pygac.calibration import Calibrator, calibrate_solar, calibrate_thermal +from pygac.calibration.noaa import Calibrator, calibrate_solar, calibrate_thermal class TestGenericCalibration(unittest.TestCase): diff --git a/pygac/tests/test_calibrate_pod.py b/pygac/tests/test_calibrate_pod.py index 67a4c72..51f4860 100644 --- a/pygac/tests/test_calibrate_pod.py +++ b/pygac/tests/test_calibrate_pod.py @@ -25,13 +25,10 @@ import unittest -try: - import mock -except ImportError: - from unittest import mock + import numpy as np -from pygac.calibration import Calibrator, calibrate_solar, calibrate_thermal +from pygac.calibration.noaa import Calibrator, calibrate_solar, calibrate_thermal class TestGenericCalibration(unittest.TestCase): diff --git a/pygac/tests/test_io.py b/pygac/tests/test_io.py index df3a4cb..6967785 100644 --- a/pygac/tests/test_io.py +++ b/pygac/tests/test_io.py @@ -20,12 +20,11 @@ """Test the I/O.""" import unittest -try: - import mock -except ImportError: - from unittest import mock +from unittest import mock + import numpy as np import numpy.testing + import pygac.gac_io as gac_io import pygac.utils as utils @@ -44,14 +43,14 @@ def test_strip_invalid_lat(self): def test_update_scanline(self): """Test updating the scanlines.""" - test_data = [{'new_start_line': 100, 'new_end_line': 200, - 'scanline': 110, 'scanline_exp': 10}, - {'new_start_line': 100, 'new_end_line': 200, - 'scanline': 90, 'scanline_exp': None}, - {'new_start_line': 100, 'new_end_line': 200, - 'scanline': 210, 'scanline_exp': None}] + test_data = [{"new_start_line": 100, "new_end_line": 200, + "scanline": 110, "scanline_exp": 10}, + {"new_start_line": 100, "new_end_line": 200, + "scanline": 90, "scanline_exp": None}, + {"new_start_line": 100, "new_end_line": 200, + "scanline": 210, "scanline_exp": None}] for t in test_data: - scanline_exp = t.pop('scanline_exp') + scanline_exp = t.pop("scanline_exp") scanline = utils._update_scanline(**t) self.assertEqual(scanline, scanline_exp) @@ -59,12 +58,12 @@ def test_update_missing_scanlines(self): """Test updating the missing scanlines.""" qual_flags = np.array([[1, 2, 4, 5, 6, 8, 9, 11, 12]]).transpose() miss_lines = np.array([3, 7, 10]) - test_data = [{'start_line': 0, 'end_line': 8, - 'miss_lines_exp': [3, 7, 10]}, - {'start_line': 3, 'end_line': 6, - 'miss_lines_exp': [1, 2, 3, 4, 7, 10, 11, 12]}] + test_data = [{"start_line": 0, "end_line": 8, + "miss_lines_exp": [3, 7, 10]}, + {"start_line": 3, "end_line": 6, + "miss_lines_exp": [1, 2, 3, 4, 7, 10, 11, 12]}] for t in test_data: - miss_lines_exp = t.pop('miss_lines_exp') + miss_lines_exp = t.pop("miss_lines_exp") miss_lines = utils._update_missing_scanlines( miss_lines=miss_lines, qual_flags=qual_flags, **t) numpy.testing.assert_array_equal(miss_lines, miss_lines_exp) @@ -185,10 +184,10 @@ def test_check_user_scanlines(self): self.assertRaises(ValueError, gac_io.check_user_scanlines, 110, 120, None, None, 100) - @mock.patch('pygac.gac_io.strip_invalid_lat') - @mock.patch('pygac.gac_io.avhrrGAC_io') - @mock.patch('pygac.gac_io.slice_channel') - @mock.patch('pygac.gac_io.check_user_scanlines') + @mock.patch("pygac.gac_io.strip_invalid_lat") + @mock.patch("pygac.gac_io.avhrrGAC_io") + @mock.patch("pygac.gac_io.slice_channel") + @mock.patch("pygac.gac_io.check_user_scanlines") def test_save_gac(self, check_user_scanlines, slice_channel, avhrr_gac_io, strip_invalid_lat): """Test saving.""" @@ -218,13 +217,13 @@ def test_save_gac(self, check_user_scanlines, slice_channel, avhrr_gac_io, qual_dir=mm, sunsatangles_dir=mm ) - slice_channel.return_value = mm, 'miss', 'midnight' + slice_channel.return_value = mm, "miss", "midnight" strip_invalid_lat.return_value = 0, 0 - check_user_scanlines.return_value = 'start', 'end' + check_user_scanlines.return_value = "start", "end" gac_io.save_gac(start_line=0, end_line=0, **kwargs) slice_channel.assert_called_with(mock.ANY, - start_line='start', end_line='end', + start_line="start", end_line="end", first_valid_lat=mock.ANY, last_valid_lat=mock.ANY ) @@ -255,8 +254,8 @@ def test_save_gac(self, check_user_scanlines, slice_channel, avhrr_gac_io, mock.ANY, mock.ANY, mock.ANY, - 'midnight', - 'miss', + "midnight", + "miss", mock.ANY, mock.ANY, mock.ANY, diff --git a/pygac/tests/test_klm.py b/pygac/tests/test_klm.py index d0bf9bf..bc66920 100644 --- a/pygac/tests/test_klm.py +++ b/pygac/tests/test_klm.py @@ -28,6 +28,7 @@ import numpy as np import numpy.testing +import xarray as xr from pygac.gac_klm import GACKLMReader from pygac.klm_reader import header @@ -44,40 +45,39 @@ def setup_method(self): def test_get_lonlat(self): """Test readout of lon/lat coordinates.""" - earth_loc = 1e4 * np.array([[1, 2, 3, 4], - [5, 6, 7, 8]]) - self.reader.scans = {'earth_location': earth_loc} - lons_exp = np.array([[2, 4], [6, 8]]) lats_exp = np.array([[1, 3], [5, 7]]) - lons, lats = self.reader._get_lonlat() + self.reader.scans = {"earth_location": {"lats": lats_exp * 1e4, + "lons": lons_exp * 1e4}} + + lons, lats = self.reader._get_lonlat_from_file() numpy.testing.assert_array_equal(lons, lons_exp) numpy.testing.assert_array_equal(lats, lats_exp) def test_get_header_timestamp(self): """Test readout of header timestamp.""" self.reader.head = { - 'start_of_data_set_year': np.array([2019]), - 'start_of_data_set_day_of_year': np.array([123]), - 'start_of_data_set_utc_time_of_day': np.array([123456]) + "start_of_data_set_year": np.array([2019]), + "start_of_data_set_day_of_year": np.array([123]), + "start_of_data_set_utc_time_of_day": np.array([123456]) } time = self.reader.get_header_timestamp() assert time == dt.datetime(2019, 5, 3, 0, 2, 3, 456000) def test_get_times(self): """Test readout of scanline timestamps.""" - self.reader.scans = {'scan_line_year': 1, - 'scan_line_day_of_year': 2, - 'scan_line_utc_time_of_day': 3} - assert self.reader._get_times() == (1, 2, 3) + self.reader.scans = {"scan_line_year": 1, + "scan_line_day_of_year": 2, + "scan_line_utc_time_of_day": 3} + assert self.reader._get_times_from_file() == (1, 2, 3) def test_get_ch3_switch(self): """Test channel 3 identification.""" self.reader.scans = { - 'scan_line_bit_field': np.array([1, 2, 3, 4, 5, 6])} + "scan_line_bit_field": np.array([1, 2, 3, 4, 5, 6])} switch_exp = np.array([1, 2, 3, 0, 1, 2]) numpy.testing.assert_array_equal( self.reader.get_ch3_switch(), switch_exp) @@ -85,7 +85,7 @@ def test_get_ch3_switch(self): def test_postproc(self): """Test KLM specific postprocessing.""" self.reader.scans = { - 'scan_line_bit_field': np.array([0, 1, 2])} + "scan_line_bit_field": np.array([0, 1, 2])} channels = np.array([[[1., 2., 3., 4.], [1., 2., 3., 4.]], [[1., 2., 3., 4.], @@ -99,16 +99,22 @@ def test_postproc(self): [1., 2., 3., np.nan]], [[1., 2., np.nan, np.nan], [1., 2., np.nan, np.nan]]]) - self.reader.postproc(channels) # masks in-place + channels = xr.DataArray(channels, dims=["scan_line_index", "columns", "channel_name"], + coords=dict(channel_name=["1", "2", "3a", "3b"] )) + ds = xr.Dataset(dict(channels=channels)) + self.reader.postproc(ds) # masks in-place numpy.testing.assert_array_equal(channels, masked_exp) def test_quality_indicators(self): """Test the quality indicator unpacking.""" reader = self.reader QFlag = reader.QFlag + + dtype = np.uint32 + quality_indicators = np.array([ 0, # nothing flagged - -1, # everything flagged + np.iinfo(dtype).max, # everything flagged QFlag.CALIBRATION | QFlag.NO_EARTH_LOCATION, QFlag.TIME_ERROR | QFlag.DATA_GAP, QFlag.FATAL_FLAG @@ -130,7 +136,7 @@ def setup_method(self): """Set up the tests.""" self.reader = GACKLMReader() - @mock.patch('pygac.klm_reader.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)) @@ -154,6 +160,8 @@ def setup_method(self): """Set up the tests.""" self.reader = LACKLMReader() self.reader.scans = np.ones(100, dtype=scanline) + self.reader.scans["scan_line_year"] = 2001 + self.reader.scans["sensor_data"] = 2047 self.reader.head = np.ones(1, dtype=header)[0] self.reader.spacecraft_id = 12 self.reader.head["noaa_spacecraft_identification_code"] = self.reader.spacecraft_id @@ -166,7 +174,7 @@ def setup_method(self): def test_get_ch3_switch(self): """Test channel 3 identification.""" self.reader.scans = { - 'scan_line_bit_field': np.array([1, 2, 3, 4, 5, 6])} + "scan_line_bit_field": np.array([1, 2, 3, 4, 5, 6])} switch_exp = np.array([1, 2, 3, 0, 1, 2]) numpy.testing.assert_array_equal( self.reader.get_ch3_switch(), switch_exp) diff --git a/pygac/tests/test_calibration_coefficients.py b/pygac/tests/test_noaa_calibration_coefficients.py similarity index 89% rename from pygac/tests/test_calibration_coefficients.py rename to pygac/tests/test_noaa_calibration_coefficients.py index e650afe..8a0c0a9 100644 --- a/pygac/tests/test_calibration_coefficients.py +++ b/pygac/tests/test_noaa_calibration_coefficients.py @@ -24,13 +24,11 @@ import sys import unittest -try: - import mock -except ImportError: - from unittest import mock +from unittest import mock + import numpy as np -from pygac.calibration import Calibrator, calibrate_solar, CoeffStatus +from pygac.calibration.noaa import Calibrator, CoeffStatus, calibrate_solar # dummy user json file including only noaa19 data with changed channel 1 coefficients user_json_file = b"""{ @@ -118,14 +116,14 @@ class TestCalibrationCoefficientsHandling(unittest.TestCase): - @mock.patch('pygac.calibration.open', mock.mock_open(read_data=user_json_file)) + @mock.patch("pygac.calibration.noaa.open", mock.mock_open(read_data=user_json_file)) def test_user_coefficients_file(self): if sys.version_info.major < 3: - cal = Calibrator('noaa19', coeffs_file="/path/to/unknow/defaults.json") + cal = Calibrator("noaa19", coeffs_file="/path/to/unknow/defaults.json") else: with self.assertWarnsRegex(RuntimeWarning, "Unknown calibration coefficients version!"): - cal = Calibrator('noaa19', coeffs_file="/path/to/unknow/defaults.json") + cal = Calibrator("noaa19", coeffs_file="/path/to/unknow/defaults.json") self.assertEqual(cal.dark_count[0], 0) self.assertEqual(cal.gain_switch[0], 1000) @@ -159,7 +157,6 @@ def test_custom_coefficients(self): if sys.version_info.major > 2: self.assertIsNone(cal.version) - @unittest.skipIf(sys.version_info.major < 3, "Skipped in python2!") def test_vis_deprecation_warning(self): counts = np.arange(10) year = 2010 @@ -182,18 +179,17 @@ def test_default_coeffs(self): _, version = Calibrator.read_coeffs(None) self.assertIsNotNone(version) - @unittest.skipIf(sys.version_info.major < 3, "Skipped in python2!") def test_read_coeffs_warnings(self): """Test warnings issued by Calibrator.read_coeffs.""" version_dicts = [ # Non-nominal coefficients - {'name': 'v123', - 'status': CoeffStatus.PROVISIONAL}, + {"name": "v123", + "status": CoeffStatus.PROVISIONAL}, # Unknown coefficients - {'name': None, - 'status': None} + {"name": None, + "status": None} ] - with mock.patch.object(Calibrator, 'version_hashs') as version_hashs: + with mock.patch.object(Calibrator, "version_hashs") as version_hashs: for version_dict in version_dicts: version_hashs.get.return_value = version_dict with self.assertWarns(RuntimeWarning): diff --git a/pygac/tests/test_pod.py b/pygac/tests/test_pod.py index 40e79b7..c006ab4 100644 --- a/pygac/tests/test_pod.py +++ b/pygac/tests/test_pod.py @@ -21,19 +21,17 @@ """Test module for the pod reading.""" import datetime as dt +import sys import unittest +from unittest import mock + import numpy as np import numpy.testing -import sys -try: - from unittest import mock -except ImportError: - import mock from pygac.clock_offsets_converter import txt as clock_offsets_txt -from pygac.reader import ReaderError, NoTLEData from pygac.gac_pod import GACPODReader from pygac.lac_pod import LACPODReader +from pygac.reader import NoTLEData, ReaderError from pygac.tests.utils import CalledWithArray @@ -51,31 +49,31 @@ def setUp(self): def test__validate_header(self): """Test the header validation""" - filename = b'NSS.GHRR.TN.D80001.S0332.E0526.B0627173.WI' - head = {'data_set_name': filename} + filename = b"NSS.GHRR.TN.D80001.S0332.E0526.B0627173.WI" + head = {"data_set_name": filename} GACPODReader._validate_header(head) # wrong name pattern with self.assertRaisesRegex(ReaderError, - 'Data set name .* does not match!'): - head = {'data_set_name': b'abc.txt'} + "Data set name .* does not match!"): + head = {"data_set_name": b"abc.txt"} GACPODReader._validate_header(head) # wrong platform - name = b'NSS.GHRR.NL.D02187.S1904.E2058.B0921517.GC' + name = b"NSS.GHRR.NL.D02187.S1904.E2058.B0921517.GC" with self.assertRaisesRegex(ReaderError, 'Improper platform id "NL"!'): - head = {'data_set_name': name} + head = {"data_set_name": name} GACPODReader._validate_header(head) # wrong transfer mode - name = filename.replace(b'GHRR', b'LHRR') + name = filename.replace(b"GHRR", b"LHRR") with self.assertRaisesRegex(ReaderError, 'Improper transfer mode "LHRR"!'): - head = {'data_set_name': name} + head = {"data_set_name": name} GACPODReader._validate_header(head) # other change reader - head = {'data_set_name': name} + head = {"data_set_name": name} LACPODReader._validate_header(head) - @mock.patch('pygac.reader.Reader.get_calibrated_channels') + @mock.patch("pygac.reader.Reader.get_calibrated_channels") def test__get_calibrated_channels_uniform_shape(self, get_channels): """Test the uniform shape as required by gac_io.save_gac.""" channels = np.arange(2*2*5, dtype=float).reshape((2, 2, 5)) @@ -97,43 +95,42 @@ def test_decode_timestamps(self): # Test whether PODReader decodes them correctly self.assertEqual(GACPODReader.decode_timestamps(t2000_enc), t2000_ref, - msg='Timestamp after 2000 was decoded incorrectly') + msg="Timestamp after 2000 was decoded incorrectly") self.assertEqual(GACPODReader.decode_timestamps(t1900_enc), t1900_ref, - msg='Timestamp before 2000 was decoded incorrectly') + msg="Timestamp before 2000 was decoded incorrectly") - @mock.patch('pygac.gac_pod.GACPODReader.decode_timestamps') + @mock.patch("pygac.gac_pod.GACPODReader.decode_timestamps") def test_get_header_timestamp(self, decode_timestamps): """Test readout of header timestamp.""" - self.reader.head = {'start_time': 123} + self.reader.head = {"start_time": 123} decode_timestamps.return_value = np.array( [2019]), np.array([123]), np.array([123456]) time = self.reader.get_header_timestamp() decode_timestamps.assert_called_with(123) self.assertEqual(time, dt.datetime(2019, 5, 3, 0, 2, 3, 456000)) - @mock.patch('pygac.gac_pod.GACPODReader.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() + self.reader.scans = {"time_code": 123} + self.reader._get_times_from_file() decode_timestamps.assert_called_with(123) def test_get_lonlat(self): """Test readout of lon/lat coordinates.""" - earth_loc = 128 * np.array([[1, 2, 3, 4], - [5, 6, 7, 8]]) - self.reader.scans = {'earth_location': earth_loc} - lons_exp = np.array([[2, 4], [6, 8]]) lats_exp = np.array([[1, 3], [5, 7]]) - lons, lats = self.reader._get_lonlat() + self.reader.scans = {"earth_location": {"lats": lats_exp * 128, + "lons": lons_exp * 128}} + + lons, lats = self.reader._get_lonlat_from_file() numpy.testing.assert_array_equal(lons, lons_exp) numpy.testing.assert_array_equal(lats, lats_exp) - @mock.patch('pygac.pod_reader.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)) @@ -158,7 +155,7 @@ def test_quality_indicators(self): QFlag.FATAL_FLAG, # 100...00 QFlag.CALIBRATION | QFlag.NO_EARTH_LOCATION, QFlag.TIME_ERROR | QFlag.DATA_GAP, - ], dtype='>u4') + ], dtype=">u4") # check if the bits look as expected bits = np.unpackbits(quality_indicators.view(np.uint8)).reshape((-1, 32)) # For a big endian integer, the number 1 fills only the last of the 32 bits @@ -189,10 +186,10 @@ def test_quality_indicators(self): expected_mask ) - @mock.patch('pygac.pod_reader.get_lonlatalt') - @mock.patch('pygac.pod_reader.compute_pixels') - @mock.patch('pygac.reader.Reader.get_tle_lines') - @mock.patch('pygac.pod_reader.avhrr_gac') + @mock.patch("pygac.pod_reader.get_lonlatalt") + @mock.patch("pygac.pod_reader.compute_pixels") + @mock.patch("pygac.reader.Reader.get_tle_lines") + @mock.patch("pygac.pod_reader.avhrr_gac") def test__adjust_clock_drift(self, avhrr_gac, get_tle_lines, compute_pixels, get_lonlatalt): """Test the clock drift adjustment.""" @@ -210,7 +207,7 @@ def test__adjust_clock_drift(self, avhrr_gac, get_tle_lines, # '1980-01-01T00:00:03.500', '1980-01-01T00:00:04.000', # '1980-01-01T00:00:04.500', '1980-01-01T00:00:05.000'] scan_utcs = ( - (1000 * scan_rate * (scan_lines - scan_lines[0])).astype('timedelta64[ms]') + (1000 * scan_rate * (scan_lines - scan_lines[0])).astype("timedelta64[ms]") + np.datetime64("1980", "ms") ) # For the geolocations, we assume an artificial swath of two pixels width @@ -230,7 +227,7 @@ def test__adjust_clock_drift(self, avhrr_gac, get_tle_lines, # '1980-01-01T00:00:00.750', '1980-01-01T00:00:01.250'] offset = 3.75 scan_offsets = offset*np.ones_like(scan_lines, dtype=float) # seconds - expected_utcs = scan_utcs - (1000*scan_offsets).astype('timedelta64[ms]') + expected_utcs = scan_utcs - (1000*scan_offsets).astype("timedelta64[ms]") # the adjustment of geolocations should keep the lons unchanged, # but should shift the lats by scan_angel * offset / scan_rate @@ -243,7 +240,7 @@ def test__adjust_clock_drift(self, avhrr_gac, get_tle_lines, # prepare the reader reader.scans = {"scan_line_number": scan_lines} - reader.utcs = scan_utcs + reader._times_as_np_datetime64 = scan_utcs reader.lons = scan_lons reader.lats = scan_lats reader.spacecraft_name = sat_name @@ -277,20 +274,20 @@ def test__adjust_clock_drift(self, avhrr_gac, get_tle_lines, # check output # use allclose for geolocations, because the slerp interpolation # includes a transormation to cartesian coordinates and back to lon, lats. - numpy.testing.assert_array_equal(reader.utcs, expected_utcs) + numpy.testing.assert_array_equal(reader._times_as_np_datetime64, expected_utcs) numpy.testing.assert_allclose(reader.lons, expected_lons) numpy.testing.assert_allclose(reader.lats, expected_lats) # undo changes to clock_offsets_txt clock_offsets_txt.pop(sat_name) - @mock.patch('pygac.pod_reader.get_offsets') - @mock.patch('pygac.reader.Reader.get_tle_lines') + @mock.patch("pygac.pod_reader.get_offsets") + @mock.patch("pygac.reader.Reader.get_tle_lines") def test__adjust_clock_drift_without_tle(self, get_tle_lines, get_offsets): """Test that clockdrift adjustment can handle missing TLE data.""" reader = self.reader - reader.utcs = np.zeros(10, dtype='datetime64[ms]') + reader._times_as_np_datetime64 = np.zeros(10, dtype="datetime64[ms]") reader.scans = {"scan_line_number": np.arange(10)} get_offsets.return_value = np.zeros(10), np.zeros(10) - get_tle_lines.side_effect = NoTLEData('No TLE data available') + get_tle_lines.side_effect = NoTLEData("No TLE data available") reader._adjust_clock_drift() # should pass without errors diff --git a/pygac/tests/test_reader.py b/pygac/tests/test_reader.py index fded545..0942552 100644 --- a/pygac/tests/test_reader.py +++ b/pygac/tests/test_reader.py @@ -29,6 +29,7 @@ import numpy as np import numpy.testing import pytest +import xarray as xr from pygac.gac_pod import scanline from pygac.gac_reader import GACReader, ReaderError @@ -40,7 +41,7 @@ from pygac.reader import NoTLEData -class TestPath(os.PathLike): +class FakePath(os.PathLike): """Fake path class.""" def __init__(self, path): @@ -59,20 +60,21 @@ class FakeGACReader(GACReader): _quality_indicators_key = "quality_indicators" tsm_affected_intervals = {None: []} along_track = 3 - across_track = 4 + across_track = 409 - def __init__(self): + def __init__(self, *args, **kwargs): """Initialize the fake reader.""" - super().__init__() + super().__init__(*args, **kwargs) self.scan_width = self.across_track scans = np.zeros(self.along_track, dtype=scanline) scans["scan_line_number"] = np.arange(self.along_track) scans["sensor_data"] = 128 self.scans = scans - self.head = {'foo': 'bar'} - self.spacecraft_name = 'noaa6' + self.head = {"foo": "bar"} + self.head = np.rec.fromrecords([("bar", ),],names="foo") + self.spacecraft_name = "noaa6" - def _get_times(self): + def _get_times_from_file(self): year = np.full(self.along_track, 1970, dtype=int) jday = np.full(self.along_track, 1, dtype=int) msec = 1000 * np.arange(1, self.along_track+1, dtype=int) @@ -85,6 +87,7 @@ def get_header_timestamp(self): def get_telemetry(self): """Get the telemetry.""" prt = 51 * np.ones(self.along_track) # prt threshold is 50 + prt[::5] = 0 ict = 101 * np.ones((self.along_track, 3)) # ict threshold is 100 space = 101 * np.ones((self.along_track, 3)) # space threshold is 100 return prt, ict, space @@ -92,8 +95,8 @@ def get_telemetry(self): def _adjust_clock_drift(self): pass - def _get_lonlat(self): - pass + def _get_lonlat_from_file(self): + return np.zeros((self.along_track, 51)), np.zeros((self.along_track, 51)) @staticmethod def _get_ir_channels_to_calibrate(): @@ -117,14 +120,23 @@ def get_tsm_pixels(self, channels): pass +class FakeGACReaderWithWrongPRTs(FakeGACReader): + + along_track = 10 + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self.scans["scan_line_number"] = np.array([1, 3, 4, 5, 6, 7, 8, 9, 10, 11]) + + class TestGacReader(unittest.TestCase): """Test the common GAC Reader.""" longMessage = True - @mock.patch.multiple('pygac.gac_reader.GACReader', + @mock.patch.multiple("pygac.gac_reader.GACReader", __abstractmethods__=set()) - @mock.patch('pygac.gac_reader.gtp.gac_lat_lon_interpolator') + @mock.patch("pygac.gac_reader.gtp.gac_lat_lon_interpolator") def setUp(self, interpolator, *mocks): """Set up the tests.""" self.interpolator = interpolator @@ -133,107 +145,107 @@ def setUp(self, interpolator, *mocks): def test_filename(self): """Test the setter of the filename property.""" # test path with .gz extension - filename = 'NSS.GHRR.TN.D80001.S0332.E0526.B0627173.WI' - filepath = '/path/to/' + filename + '.gz' + filename = "NSS.GHRR.TN.D80001.S0332.E0526.B0627173.WI" + filepath = "/path/to/" + filename + ".gz" self.reader.filename = filepath self.assertEqual(self.reader.filename, filename) self.reader.filename = None self.assertIsNone(self.reader.filename) - self.reader.filename = TestPath(filepath) + self.reader.filename = FakePath(filepath) self.assertEqual(self.reader.filename, filename) @unittest.skipIf(sys.version_info.major < 3, "Skipped in python2!") def test__read_scanlines(self): """Test the scanline extraction.""" self.reader.scanline_type = np.dtype([ - ('a', 'S2'), ('b', 'M8[ms]") + self.reader._times_as_np_datetime64 = msecs.astype(">M8[ms]") self.reader.scans = np.array(scan_line_numbers, dtype=[("scan_line_number", ">u2")]) # Test correction self.reader.correct_times_thresh() - numpy.testing.assert_array_equal(self.reader.utcs, utcs_expected) + numpy.testing.assert_array_equal(self.reader._times_as_np_datetime64, utcs_expected) def test_calculate_sun_earth_distance_correction(self): """Test the calculate sun earth distance correction method.""" - self.reader.utcs = np.array([315748035469, 315748359969, + self.reader._times_as_np_datetime64 = np.array([315748035469, 315748359969, 315751135469, 315754371969, - 315754371969]).astype('datetime64[ms]') + 315754371969]).astype("datetime64[ms]") corr = self.reader.get_sun_earth_distance_correction() numpy.testing.assert_almost_equal(corr, 0.96660494, decimal=7) - @mock.patch('pygac.reader.Reader.get_sun_earth_distance_correction') - @mock.patch('pygac.reader.Reader.get_midnight_scanline') - @mock.patch('pygac.reader.Reader.get_miss_lines') - @mock.patch('pygac.reader.Reader.calibration', - new_callable=mock.PropertyMock) + @mock.patch("pygac.reader.Reader.get_sun_earth_distance_correction") + @mock.patch("pygac.reader.Reader.get_midnight_scanline") + @mock.patch("pygac.reader.Reader.get_miss_lines") def test_update_metadata(self, - calibration, get_miss_lines, get_midnight_scanline, get_sun_earth_distance_correction): """Test updating the metadata.""" - get_miss_lines.return_value = 'miss_lines' - get_midnight_scanline.return_value = 'midn_line' - get_sun_earth_distance_correction.return_value = 'factor' - self.reader.head = {'foo': 'bar'} - calibration.return_value = mock.MagicMock(version='version') + get_miss_lines.return_value = "miss_lines" + get_midnight_scanline.return_value = "midn_line" + get_sun_earth_distance_correction.return_value = "factor" + self.reader.head = {"foo": "bar"} self.reader.update_meta_data() - mda_exp = {'midnight_scanline': 'midn_line', - 'missing_scanlines': 'miss_lines', - 'sun_earth_distance_correction_factor': 'factor', - 'gac_header': {'foo': 'bar'}, - 'calib_coeffs_version': 'version'} + mda_exp = {"midnight_scanline": "midn_line", + "missing_scanlines": "miss_lines", + "sun_earth_distance_correction_factor": "factor", + "gac_header": {"foo": "bar"}} self.assertDictEqual(self.reader.meta_data, mda_exp) @@ -662,7 +677,7 @@ class TestLacReader(unittest.TestCase): longMessage = True - @mock.patch.multiple('pygac.lac_reader.LACReader', + @mock.patch.multiple("pygac.lac_reader.LACReader", __abstractmethods__=set()) def setUp(self, *mocks): """Set up the tests.""" @@ -670,7 +685,7 @@ def setUp(self, *mocks): def test_lac_reader_accepts_FRAC(self): """Test the header validation.""" - head = {'data_set_name': b'NSS.FRAC.M1.D19115.S2352.E0050.B3425758.SV'} + head = {"data_set_name": b"NSS.FRAC.M1.D19115.S2352.E0050.B3425758.SV"} self.reader._validate_header(head) def test_correct_scan_line_numbers(self): @@ -678,7 +693,7 @@ def test_correct_scan_line_numbers(self): scans, expected = _get_scanline_numbers(22000) self.reader.scans = scans self.reader.correct_scan_line_numbers() - numpy.testing.assert_array_equal(self.reader.scans['scan_line_number'], + numpy.testing.assert_array_equal(self.reader.scans["scan_line_number"], expected) @@ -694,12 +709,12 @@ def pod_file_with_tbm_header(tmp_path): tbm_header["ending_latitude"] = b"+22" tbm_header["beginning_longitude"] = b"-004" tbm_header["ending_longitude"] = b"+032" - tbm_header["start_hour"] = b'AL' - tbm_header["start_minute"] = b'L ' - tbm_header["number_of_minutes"] = b'ALL' - tbm_header["appended_data_flag"] = b'Y' - tbm_header["channel_select_flag"][0, :5] = b'\x01' - tbm_header["sensor_data_word_size"] = b'10' + tbm_header["start_hour"] = b"AL" + tbm_header["start_minute"] = b"L " + tbm_header["number_of_minutes"] = b"ALL" + tbm_header["appended_data_flag"] = b"Y" + tbm_header["channel_select_flag"][0, :5] = b"\x01" + tbm_header["sensor_data_word_size"] = b"10" header = np.zeros(1, dtype=header3) header["noaa_spacecraft_identification_code"] = 3 @@ -707,7 +722,7 @@ def pod_file_with_tbm_header(tmp_path): header["start_time"] = [51522, 181, 62790] header["number_of_scans"] = number_of_scans header["end_time"] = [51522, 195, 42286] - header["processing_block_id"] = b'3031919' + header["processing_block_id"] = b"3031919" header["ramp_auto_calibration"] = 0 header["number_of_data_gaps"] = 0 header["dacs_quality"] = [21, 7, 0, 0, 0, 0] @@ -745,23 +760,34 @@ def pod_file_with_tbm_header(tmp_path): scanlines["calibration_coefficients"] = [149722896, -23983736, 173651248, -27815402, -1803673, 6985664, -176321840, 666680320, -196992576, 751880640] scanlines["number_of_meaningful_zenith_angles_and_earth_location_appended"] = 51 - scanlines["solar_zenith_angles"] = [-24, -26, -28, -30, -32, -33, -34, -35, -37, -37, -38, -39, -40, -41, -41, -42, - -43, -43, -44, -45, -45, -46, -46, -47, -48, -48, -49, -49, -50, -50, -51, -52, - -52, -53, -53, -54, -55, -55, -56, -57, -58, -59, -60, -61, -62, -63, -64, -66, - -68, -70, -73] - scanlines["earth_location"] = [9929, -36, 9966, 747, 9983, 1415, 9988, 1991, 9984, 2492, 9975, 2931, - 9963, 3320, 9948, 3667, 9931, 3979, 9913, 4262, 9895, 4519, 9875, 4755, - 9856, 4972, 9836, 5174, 9816, 5362, 9796, 5538, 9775, 5703, 9755, 5860, - 9734, 6009, 9713, 6150, 9692, 6286, 9671, 6416, 9650, 6542, 9628, 6663, - 9606, 6781, 9583, 6896, 9560, 7009, 9537, 7119, 9513, 7227, 9488, 7334, - 9463, 7440, 9437, 7545, 9409, 7650, 9381, 7755, 9351, 7860, 9320, - 7966, 9288, 8073, 9253, 8181, 9216, 8291, 9177, 8403, 9135, 8518, - 9090, 8637, 9040, 8759, 8986, 8886, 8926, 9019, 8859, 9159, 8784, - 9307, 8697, 9466, 8596, 9637, 8476, 9824, 8327, 10033] - scanlines["telemetry"] = 0 + + one_line_angles = np.array([-24, -26, -28, -30, -32, -33, -34, -35, -37, -37, -38, -39, -40, -41, -41, -42, + -43, -43, -44, -45, -45, -46, -46, -47, -48, -48, -49, -49, -50, -50, -51, -52, + -52, -53, -53, -54, -55, -55, -56, -57, -58, -59, -60, -61, -62, -63, -64, -66, + -68, -70, -73]) + scanlines["solar_zenith_angles"] = [one_line_angles, one_line_angles + 1, one_line_angles + 2] + + one_line_lats = np.array([9929, 9966, 9983, 9988, 9984, 9975, 9963, 9948, 9931, 9913, 9895, 9875, + 9856, 9836, 9816, 9796, 9775, 9755, 9734, 9713, 9692, 9671, 9650, 9628, + 9606, 9583, 9560, 9537, 9513, 9488, 9463, 9437, 9409, 9381, 9351, 9320, + 9288, 9253, 9216, 9177, 9135, 9090, 9040, 8986, 8926, 8859, 8784, 8697, + 8596, 8476, 8327]) + + scanlines["earth_location"]["lats"] = [one_line_lats - 1, one_line_lats, one_line_lats + 1] + + one_line_lons = np.array([-36, 747, 1415, 1991, 2492, 2931, 3320, 3667, 3979, 4262, 4519, 4755, 4972, + 5174, 5362, 5538, 5703, 5860, 6009, 6150, 6286, 6416, 6542, 6663, 6781, + 6896, 7009, 7119, 7227, 7334, 7440, 7545, 7650, 7755, 7860, 7966, 8073, + 8181, 8291, 8403, 8518, 8637, 8759, 8886, 9019, 9159, 9307, 9466, 9637, + 9824, 10033]) + + scanlines["earth_location"]["lons"] = [one_line_lons, one_line_lons + 1, one_line_lons + 2] + + scanlines["telemetry"] = 2047 scanlines["sensor_data"] = 99 scanlines["add_on_zenith"] = 0 scanlines["clock_drift_delta"] = 0 + scanlines["sensor_data"] = 2047 pod_filename = tmp_path / "image.l1b" offset = 14800 @@ -775,6 +801,18 @@ def pod_file_with_tbm_header(tmp_path): return pod_filename +@pytest.fixture() +def pod_tle(tmp_path): + lines = ("1 23455U 94089A 00322.04713399 .00000318 00000-0 19705-3 0 5298\n" + "2 23455 99.1591 303.5706 0010037 25.7760 334.3905 14.12496755303183\n" + "1 23455U 94089A 00322.96799836 .00000229 00000-0 14918-3 0 5303\n" + "2 23455 99.1590 304.5117 0009979 23.1101 337.0518 14.12496633303313\n") + tle_filename = tmp_path / "noaa14.tle" + with tle_filename.open("w") as fd: + fd.write(lines) + return tle_filename + + def test_podlac_eosip(pod_file_with_tbm_header): """Test reading a real podlac file.""" reader = LACPODReader(interpolate_coords=False) @@ -782,3 +820,61 @@ def test_podlac_eosip(pod_file_with_tbm_header): assert reader.head.itemsize == header3.itemsize # this is broken in eosip pod data, tbm data set name has start and end times reversed. # assert reader.head["data_set_name"] == reader.tbm_head["data_set_name"] + # todo: test that duplicate lines are removed and recomputed + + +def test_read_to_dataset_is_a_dataset_including_channels_and_telemetry(pod_file_with_tbm_header, pod_tle): + """Test creating an xr.Dataset from a gac file.""" + import xarray as xr + reader = LACPODReader(interpolate_coords=False, tle_dir=pod_tle.parent, tle_name=os.path.basename(pod_tle)) + dataset = reader.read_as_dataset(pod_file_with_tbm_header) + assert isinstance(dataset, xr.Dataset) + assert dataset["channels"].shape == (3, 2048, 5) + assert dataset.attrs["processing_block_id"] == b"3031919" + assert "times" in dataset.coords + assert "scan_line_index" in dataset.coords + assert "channel_name" in dataset.coords + assert dataset["prt_counts"].shape == (3,) + assert dataset["ict_counts"].shape == (3, 3) + assert dataset["space_counts"].shape == (3, 3) + + +def test_read_to_dataset_without_interpolation(pod_file_with_tbm_header, pod_tle): + """Test creating an xr.Dataset from a gac file.""" + reader = LACPODReader(interpolate_coords=False, tle_dir=pod_tle.parent, tle_name=os.path.basename(pod_tle)) + dataset = reader.read_as_dataset(pod_file_with_tbm_header) + + assert dataset["longitude"].shape == (3, 51) + assert dataset["latitude"].shape == (3, 51) + + +def test_read_to_dataset_with_interpolation(pod_file_with_tbm_header, pod_tle): + """Test creating an xr.Dataset from a gac file.""" + reader = LACPODReader(interpolate_coords=True, tle_dir=pod_tle.parent, tle_name=os.path.basename(pod_tle)) + dataset = reader.read_as_dataset(pod_file_with_tbm_header) + + assert dataset.coords["longitude"].shape == (3, 2048) + assert dataset.coords["latitude"].shape == (3, 2048) + + +def test_passing_calibration_coeffs_to_reader_init_is_deprecated(): + """Test that passing calibration-specific info to the reader's init is not allowed anymore.""" + with pytest.warns(PendingDeprecationWarning): + FakeGACReader(calibration_method="noaa", custom_calibration=dict(a="1")) + with pytest.warns(PendingDeprecationWarning): + FakeGACReader(calibration_method="noaa", calibration_file="somefile") + + +def test_passing_calibration_to_reader(): + """Test passing calibration info to `get_calibrated_channels`.""" + method = "InvalidMethod" + with pytest.raises(ValueError, match=method): + reader = FakeGACReader(calibration_method=method) + + reader = FakeGACReader(calibration_method="noaa") + + res = reader.get_calibrated_channels() + + np.testing.assert_allclose(res[:, 1, 0], 8.84714652) + np.testing.assert_allclose(res[:, 2, 1], 10.23511303) + assert reader.meta_data["calib_coeffs_version"] == "PATMOS-x, v2023" diff --git a/pygac/tests/test_slerp.py b/pygac/tests/test_slerp.py index f225631..f48201f 100644 --- a/pygac/tests/test_slerp.py +++ b/pygac/tests/test_slerp.py @@ -25,9 +25,11 @@ import unittest -from pygac.slerp import slerp + import numpy as np +from pygac.slerp import slerp + class TestSlerp(unittest.TestCase): diff --git a/pygac/tests/test_tsm.py b/pygac/tests/test_tsm.py index 5a21d10..45d3fe8 100644 --- a/pygac/tests/test_tsm.py +++ b/pygac/tests/test_tsm.py @@ -21,11 +21,12 @@ """Test TSM module.""" import unittest -import pygac.correct_tsm_issue as tsm import numpy as np import numpy.testing +import pygac.correct_tsm_issue as tsm + class TSMTest(unittest.TestCase): """Test TSM module.""" diff --git a/pygac/tests/test_utils.py b/pygac/tests/test_utils.py index 3361ca1..5fcf411 100644 --- a/pygac/tests/test_utils.py +++ b/pygac/tests/test_utils.py @@ -31,7 +31,7 @@ import numpy as np -from pygac.utils import file_opener, calculate_sun_earth_distance_correction +from pygac.utils import calculate_sun_earth_distance_correction, file_opener class TestUtils(unittest.TestCase): @@ -39,20 +39,20 @@ class TestUtils(unittest.TestCase): longMessage = True - @mock.patch('pygac.utils.open', mock.MagicMock(return_value=io.BytesIO(b'file content'))) + @mock.patch("pygac.utils.open", mock.MagicMock(return_value=io.BytesIO(b"file content"))) def test_file_opener_1(self): """Test if a file is redirected correctly through file_opener.""" - with file_opener('path/to/file') as f: + with file_opener("path/to/file") as f: content = f.read() - self.assertEqual(content, b'file content') + self.assertEqual(content, b"file content") def test_file_opener_2(self): """Test file_opener with file objects and compression""" # prepare test - normal_message = b'normal message' - gzip_message_decoded = b'gzip message' + normal_message = b"normal message" + gzip_message_decoded = b"gzip message" with io.BytesIO() as f: - with gzip.open(f, mode='wb') as g: + with gzip.open(f, mode="wb") as g: g.write(gzip_message_decoded) f.seek(0) gzip_message_encoded = f.read() @@ -68,7 +68,7 @@ def test_file_opener_2(self): message = g.read() self.assertEqual(message, gzip_message_decoded) - @mock.patch('pygac.utils.open', mock.MagicMock(side_effect=FileNotFoundError)) + @mock.patch("pygac.utils.open", mock.MagicMock(side_effect=FileNotFoundError)) def test_file_opener_3(self): """Test file_opener with PathLike object""" # prepare test @@ -83,8 +83,8 @@ def __fspath__(self): def open(self): return io.BytesIO(self.raw_bytes) - filename = '/path/to/file' - file_bytes = b'TestTestTest' + filename = "/path/to/file" + file_bytes = b"TestTestTest" test_pathlike = RawBytes(filename, file_bytes) with file_opener(test_pathlike) as f: content = f.read() diff --git a/pygac/utils.py b/pygac/utils.py index 76a7fd9..e35b564 100644 --- a/pygac/utils.py +++ b/pygac/utils.py @@ -22,7 +22,6 @@ import gzip import io import logging - from contextlib import contextmanager, nullcontext import numpy as np @@ -33,7 +32,7 @@ def gzip_inspected(open_file): """Try to gzip decompress the file object if applicable.""" try: - file_object = gzip.GzipFile(mode='rb', fileobj=open_file) + file_object = gzip.GzipFile(mode="rb", fileobj=open_file) file_object.read(1) except OSError: file_object = open_file @@ -47,13 +46,13 @@ def file_opener(file): if isinstance(file, io.IOBase) and file.seekable(): # avoid closing the file using nullcontext open_file = nullcontext(file) - elif hasattr(file, 'open'): + elif hasattr(file, "open"): try: - open_file = file.open(mode='rb') + open_file = file.open(mode="rb") except TypeError: open_file = file.open() else: - open_file = open(file, mode='rb') + open_file = open(file, mode="rb") # set open_file into context in case of lazy loading in __enter__ method. with open_file as file_object: yield gzip_inspected(file_object) @@ -108,7 +107,7 @@ def check_user_scanlines(start_line, end_line, first_valid_lat=None, num_valid_lines = last_valid_lat - first_valid_lat + 1 else: if along_track is None: - raise ValueError('Need along_track') + raise ValueError("Need along_track") num_valid_lines = along_track start_line = int(start_line) @@ -119,10 +118,10 @@ def check_user_scanlines(start_line, end_line, first_valid_lat=None, end_line = num_valid_lines - 1 elif end_line >= num_valid_lines: end_line = num_valid_lines - 1 - LOG.warning('Given end line exceeds scanline range, resetting ' - 'to {}'.format(end_line)) + LOG.warning("Given end line exceeds scanline range, resetting " + "to {}".format(end_line)) if start_line > num_valid_lines: - raise ValueError('Given start line {} exceeds scanline range {}' + raise ValueError("Given start line {} exceeds scanline range {}" .format(start_line, num_valid_lines)) return start_line, end_line @@ -199,9 +198,9 @@ def _slice(ch, start_line, end_line, update=None): ch_slc = ch[start_line:end_line + 1, :].copy() if update: - updated = [_update_scanline(l, start_line, end_line) - if l is not None else None - for l in update] + updated = [_update_scanline(line, start_line, end_line) + if line is not None else None + for line in update] return ch_slc, updated return ch_slc @@ -239,13 +238,13 @@ def plot_correct_times_thresh(res, filename=None): """Visualize results of GACReader.correct_times_thresh.""" import matplotlib.pyplot as plt - t = res['t'] - tcorr = res.get('tcorr') - n = res['n'] - offsets = res.get('offsets') - t0_head = res.get('t0_head') - max_diff_from_t0_head = res.get('max_diff_from_t0_head') - fail_reason = res.get('fail_reason', 'Failed for unknown reason') + t = res["t"] + tcorr = res.get("tcorr") + n = res["n"] + offsets = res.get("offsets") + t0_head = res.get("t0_head") + max_diff_from_t0_head = res.get("max_diff_from_t0_head") + fail_reason = res.get("fail_reason", "Failed for unknown reason") # Setup figure along_track = np.arange(t.size) @@ -295,13 +294,13 @@ def plot_correct_scanline_numbers(res, filename=None): """Visualize results of GACReader.correct_scanline_numbers.""" import matplotlib.pyplot as plt - along_track = res['along_track'] - n_orig = res['n_orig'] - n_corr = res['n_corr'] - within_range = res['within_range'] - thresh = res['thresh'] - diffs = res['diffs'] - nz_diffs = res['nz_diffs'] + along_track = res["along_track"] + n_orig = res["n_orig"] + n_corr = res["n_corr"] + within_range = res["within_range"] + thresh = res["thresh"] + diffs = res["diffs"] + nz_diffs = res["nz_diffs"] # Setup figure _, (ax0, ax1) = plt.subplots(nrows=2) @@ -327,6 +326,6 @@ def plot_correct_scanline_numbers(res, filename=None): plt.tight_layout() if filename: - plt.savefig(filename, bbox_inches='tight') + plt.savefig(filename, bbox_inches="tight") else: plt.show() diff --git a/pyproject.toml b/pyproject.toml index c10a1f8..28cbf74 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,7 +13,8 @@ dependencies = [ "pyorbital>=0.3.2", "python-geotiepoints>=1.1.8", "scipy>=0.8.0", - "packaging" + "packaging", + "xarray>=2024.1.0" ] readme = "README.md" requires-python = ">=3.8" @@ -43,6 +44,9 @@ dev = ['pytest', 'pre-commit', 'ruff'] [project.scripts] pygac-convert-patmosx-coefficients = "pygac.patmosx_coeff_reader:main" +[project.entry-points."pygac.calibration"] +noaa = "pygac.calibration.noaa:calibrate" + [build-system] requires = ["hatchling", "hatch-vcs"] build-backend = "hatchling.build" @@ -63,4 +67,4 @@ version-file = "pygac/version.py" line-length = 120 [tool.ruff.lint] -select = ["E", "W", "F", "I"] +select = ["E", "W", "F", "I", "Q"] diff --git a/requirements.txt b/requirements.txt index 82c43ee..e9c71f2 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,5 +3,4 @@ h5py>=2.0.1 scipy>=0.8.0 python-geotiepoints>=1.1.8 packaging>=14.0 - - +xarray>=2024.1.0