diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 644d3aa6ae..13444de46c 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -103,7 +103,8 @@ jobs: git+https://github.com/Unidata/cftime \ git+https://github.com/mapbox/rasterio \ git+https://github.com/pydata/bottleneck \ - git+https://github.com/pydata/xarray; + git+https://github.com/pydata/xarray \ + git+https://github.com/astropy/astropy; - name: Install satpy shell: bash -l {0} diff --git a/AUTHORS.md b/AUTHORS.md index dd2b24750d..2bb496e52e 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -25,13 +25,13 @@ The following people have made contributions to this project: - [Adam Dybbroe (adybbroe)](https://github.com/adybbroe) - [Ulrik Egede (egede)](https://github.com/egede) - [Joleen Feltz (joleenf)](https://github.com/joleenf) -- [Stephan Finkensieper (sfinkens)](https://github.com/sfinkens) +- [Stephan Finkensieper (sfinkens)](https://github.com/sfinkens) - Deutscher Wetterdienst - [Andrea Grillini (AppLEaDaY)](https://github.com/AppLEaDaY) - [Blanka Gvozdikova (gvozdikb)](https://github.com/gvozdikb) - [Nina Håkansson (ninahakansson)](https://github.com/ninahakansson) - [Ulrich Hamann](https://github.com/) - [Mitch Herbertson (mherbertson)](https://github.com/mherbertson) -- [Gerrit Holl (gerritholl)](https://github.com/gerritholl) +- [Gerrit Holl (gerritholl)](https://github.com/gerritholl) - Deutscher Wetterdienst - [David Hoese (djhoese)](https://github.com/djhoese) - [Marc Honnorat (honnorat)](https://github.com/honnorat) - [Mikhail Itkin (mitkin)](https://github.com/mitkin) @@ -42,7 +42,8 @@ The following people have made contributions to this project: - [Janne Kotro (jkotro)](https://github.com/jkotro) - [Ralph Kuehn (ralphk11)](https://github.com/ralphk11) - [Panu Lahtinen (pnuu)](https://github.com/pnuu) -- [Thomas Leppelt (m4sth0)](https://github.com/m4sth0) +- [Jussi Leinonen (jleinonen)](https://github.com/jleinonen) - meteoswiss +- [Thomas Leppelt (m4sth0)](https://github.com/m4sth0) - Deutscher Wetterdienst - [Lu Liu (yukaribbba)](https://github.com/yukaribbba) - [Andrea Meraner (ameraner)](https://github.com/ameraner) - [Aronne Merrelli (aronnem)](https://github.com/aronnem) diff --git a/continuous_integration/environment.yaml b/continuous_integration/environment.yaml index ad768209b3..b86bc1e072 100644 --- a/continuous_integration/environment.yaml +++ b/continuous_integration/environment.yaml @@ -17,7 +17,7 @@ dependencies: - scipy - pyyaml - pyproj - - pyresample + - pyresample>=1.24 - coveralls - coverage - codecov @@ -49,6 +49,8 @@ dependencies: - python-geotiepoints - pooch - pip + - skyfield + - astropy - pip: - trollsift - trollimage diff --git a/doc/source/modifiers.rst b/doc/source/modifiers.rst index 4669aca145..4869fe9091 100644 --- a/doc/source/modifiers.rst +++ b/doc/source/modifiers.rst @@ -55,3 +55,84 @@ on those modifiers can be found in the linked API documentation. A complete list can be found in the `etc/composites `_ source code and in the :mod:`~satpy.modifiers` module documentation. + +Parallax correction +------------------- + +.. warning:: + + The Satpy parallax correction is experimental and subject to change. + +Since version 0.37 (mid 2022), Satpy has included a +modifier for parallax correction, implemented in the +:class:`~satpy.modifiers.parallax.ParallaxCorrectionModifier` class. +This modifier is important for some applications, but not applied +by default to any Satpy datasets or composites, because it can be +applied to any input dataset and used with any source of (cloud top) +height. Therefore, users wishing to apply the parallax correction +semi-automagically have to define their own modifier and then apply +that modifier for their datasets. An example is included +with the :class:`~satpy.modifiers.parallax.ParallaxCorrectionModifier` +API documentation. Note that Satpy cannot apply modifiers to +composites, so users wishing to apply parallax correction to a composite +will have to use a lower level API or duplicate an existing composite +recipe to use modified inputs. + +The parallax correction is directly calculated from the cloud top height. +Information on satellite position is obtained from cloud top height +metadata. If no orbital parameters are present in the cloud top height +metadata, Satpy will attempt to calculate orbital parameters from the +platform name and start time. The backup calculation requires skyfield +and astropy to be installed. If the metadata include neither orbital +parameters nor platform name and start time, parallax calculation will +fail. Because the cloud top height metadata are used, it is essential +that the cloud top height data are derived from the same platform as +the measurements to be corrected are taken by. + +The parallax error moves clouds away from the observer. Therefore, the +parallax correction shifts clouds in the direction of the observer. The +space left behind by the cloud will be filled with fill values. As the +cloud is shifted toward the observer, it may occupy less pixels than before, +because pixels closer to the observer have a smaller surface area. It can +also be deformed (a "rectangular" cloud may get the shape of a parallelogram). + +.. figure:: https://figshare.com/ndownloader/files/36422616/preview/36422616/preview.jpg + :width: 512 + :height: 512 + :alt: Satellite image without parallax correction. + + SEVIRI view of southern Sweden, 2021-11-30 12:15Z, without parallax correction. + This is the ``natural_color`` composite as built into Satpy. + + +.. figure:: https://figshare.com/ndownloader/files/36422613/preview/36422613/preview.jpg + :width: 512 + :height: 512 + :alt: Satellite image with parallax correction. + + The same satellite view with parallax correction. The most obvious change + are the gaps left behind by the parallax correction, shown as black pixels. + Otherwise it shows that clouds have "moved" south-south-west in the direction + of the satellite. To view the images side-by-side or alternating, look at + `the figshare page `_ + +The utility function :func:`~satpy.modifiers.parallax.get_surface_parallax_displacement` allows to calculate the magnitude of the parallax error. For a cloud with a cloud top height of 10 km: + +.. figure:: https://figshare.com/ndownloader/files/36462435/preview/36462435/preview.jpg + :width: 512 + :height: 512 + :alt: Figure showing magnitude of parallax effect. + + Magnitude of the parallax error for a fictitious cloud with a cloud top + height of 10 km for the GOES-East (GOES-16) full disc. + +The parallax correction is currently experimental and subject to change. +Although it is covered by tests, there may be cases that yield unexpected +or incorrect results. It does not yet perform any checks that the +provided (cloud top) height covers the area of the dataset for which +the parallax correction shall be applied. + +For more general background information and web routines related to the +parallax effect, see also `this collection at the CIMSS website _`. + +.. versionadded:: 0.37 diff --git a/satpy/modifiers/parallax.py b/satpy/modifiers/parallax.py new file mode 100644 index 0000000000..11185f2edb --- /dev/null +++ b/satpy/modifiers/parallax.py @@ -0,0 +1,551 @@ +# Copyright (c) 2021-2022 Satpy developers +# +# This file is part of satpy. +# +# satpy is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# satpy is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +# A PARTICULAR PURPOSE. See the GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License along with +# satpy. If not, see . +"""Parallax correction. + +Routines related to parallax correction using datasets involving height, such +as cloud top height. + +The geolocation of (geostationary) satellite imagery is calculated by +agencies or in satpy readers with the assumption of a clear view from +the satellite to the geoid. When a cloud blocks the view of the Earth +surface or the surface is above sea level, the geolocation is not accurate +for the cloud or mountain top. This module contains routines to correct +imagery such that pixels are shifted or interpolated to correct for this +parallax effect. + +Parallax correction is currently only supported for (cloud top) height +that arrives on an :class:`~pyresample.geometry.AreaDefinition`, such +as is standard for geostationary satellites. Parallax correction with +data described by a :class:`~pyresample.geometry.SwathDefinition`, +such as is common for polar satellites, is not (yet) supported. + +See also the :doc:`../modifiers` page in the documentation for an introduction to +parallax correction as a modifier in Satpy. +""" + +import datetime +import inspect +import logging +import warnings + +import dask.array as da +import numpy as np +import xarray as xr +from pyorbital.orbital import A as EARTH_RADIUS +from pyorbital.orbital import get_observer_look +from pyproj import Geod +from pyresample.bucket import BucketResampler +from pyresample.geometry import SwathDefinition + +from satpy.modifiers import ModifierBase +from satpy.resample import resample_dataset +from satpy.utils import get_satpos, lonlat2xyz, xyz2lonlat + +logger = logging.getLogger(__name__) + + +class MissingHeightError(ValueError): + """Raised when heights do not overlap with area to be corrected.""" + + +class IncompleteHeightWarning(UserWarning): + """Raised when heights only partially overlap with area to be corrected.""" + + +def get_parallax_corrected_lonlats(sat_lon, sat_lat, sat_alt, lon, lat, height): + """Calculate parallax corrected lon/lats. + + Satellite geolocation generally assumes an unobstructed view of a smooth + Earth surface. In reality, this view may be obstructed by clouds or + mountains. + + If the view of a pixel at location (lat, lon) is blocked by a cloud + at height h, this function calculates the (lat, lon) coordinates + of the cloud above/in front of the invisible surface. + + For scenes that are only partly cloudy, the user might set the cloud top + height for clear-sky pixels to NaN. This function will return a corrected + lat/lon as NaN as well. The user can use the original lat/lon for those + pixels or use the higher level :class:`ParallaxCorrection` class. + + This function assumes a spherical Earth. + + .. note:: + + Be careful with units! This code expects ``sat_alt`` and + ``height`` to be in meter above the Earth's surface. You may + have to convert your input correspondingly. Cloud Top Height + is usually reported in meters above the Earth's surface, rarely + in km. Satellite altitude may be reported in either m or km, but + orbital parameters are usually in relation to the Earth's centre. + The Earth radius from pyresample is reported in km. + + Args: + sat_lon (number): Satellite longitude in geodetic coordinates [degrees] + sat_lat (number): Satellite latitude in geodetic coordinates [degrees] + sat_alt (number): Satellite altitude above the Earth surface [m] + lon (array or number): Longitudes of pixel or pixels to be corrected, + in geodetic coordinates [degrees] + lat (array or number): Latitudes of pixel/pixels to be corrected, in + geodetic coordinates [degrees] + height (array or number): Heights of pixels on which the correction + will be based. Typically this is the cloud top height. [m] + + Returns: + tuple[float, float]: Corrected geolocation + Corrected geolocation ``(lon, lat)`` in geodetic coordinates for + the pixel(s) to be corrected. [degrees] + """ + elevation = _get_satellite_elevation(sat_lon, sat_lat, sat_alt, lon, lat) + parallax_distance = _calculate_slant_cloud_distance(height, elevation) + shifted_xyz = _get_parallax_shift_xyz( + sat_lon, sat_lat, sat_alt, lon, lat, parallax_distance) + + return xyz2lonlat( + shifted_xyz[..., 0], shifted_xyz[..., 1], shifted_xyz[..., 2]) + + +def get_surface_parallax_displacement( + sat_lon, sat_lat, sat_alt, lon, lat, height): + """Calculate surface parallax displacement. + + Calculate the displacement due to parallax error. Input parameters are + identical to :func:`get_parallax_corrected_lonlats`. + + Returns: + number or array: parallax displacement in meter + """ + (corr_lon, corr_lat) = get_parallax_corrected_lonlats(sat_lon, sat_lat, sat_alt, lon, lat, height) + # Get parallax displacement + geod = Geod(ellps="sphere") + _, _, parallax_dist = geod.inv(corr_lon, corr_lat, lon, lat) + return parallax_dist + + +def _get_parallax_shift_xyz(sat_lon, sat_lat, sat_alt, lon, lat, parallax_distance): + """Calculate the parallax shift in cartesian coordinates. + + From satellite position and cloud position, get the parallax shift in + cartesian coordinates: + + Args: + sat_lon (number): Satellite longitude in geodetic coordinates [degrees] + sat_lat (number): Satellite latitude in geodetic coordinates [degrees] + sat_alt (number): Satellite altitude above the Earth surface [m] + lon (array or number): Longitudes of pixel or pixels to be corrected, + in geodetic coordinates [degrees] + lat (array or number): Latitudes of pixel/pixels to be corrected, in + geodetic coordinates [degrees] + parallax_distance (array or number): Cloud to ground distance with parallax + effect [m]. + + Returns: + Parallax shift in cartesian coordinates in meter. + """ + sat_xyz = np.hstack(lonlat2xyz(sat_lon, sat_lat)) * sat_alt + cth_xyz = np.stack(lonlat2xyz(lon, lat), axis=-1) * EARTH_RADIUS*1e3 # km → m + delta_xyz = cth_xyz - sat_xyz + sat_distance = np.sqrt((delta_xyz*delta_xyz).sum(axis=-1)) + dist_shape = delta_xyz.shape[:-1] + (1,) # force correct array broadcasting + return cth_xyz - delta_xyz*(parallax_distance/sat_distance).reshape(dist_shape) + + +def _get_satellite_elevation(sat_lon, sat_lat, sat_alt, lon, lat): + """Get satellite elevation. + + Get the satellite elevation from satellite lon/lat/alt for positions + lon/lat. + """ + placeholder_date = datetime.datetime(2000, 1, 1) # no impact on get_observer_look? + (_, elevation) = get_observer_look( + sat_lon, sat_lat, sat_alt/1e3, # m → km (wanted by get_observer_look) + placeholder_date, lon, lat, 0) + return elevation + + +def _calculate_slant_cloud_distance(height, elevation): + """Calculate slant cloud to ground distance. + + From (cloud top) height and satellite elevation, calculate the + slant cloud-to-ground distance along the line of sight of the satellite. + """ + if np.isscalar(elevation) and elevation == 0: + raise NotImplementedError( + "Parallax correction not implemented for " + "satellite elevation 0") + if np.isscalar(elevation) and elevation < 0: + raise ValueError( + "Satellite is below the horizon. Cannot calculate parallax " + "correction.") + return height / np.sin(np.deg2rad(elevation)) + + +class ParallaxCorrection: + """Parallax correction calculations. + + This class contains higher-level functionality to wrap the parallax + correction calculations in :func:`get_parallax_corrected_lonlats`. The class is + initialised using a base area, which is the area for which a corrected + geolocation will be calculated. The resulting object is a callable. + Calling the object with an array of (cloud top) heights returns a + :class:`~pyresample.geometry.SwathDefinition` describing the new , + corrected geolocation. The cloud top height should cover at least the + area for which the corrected geolocation will be calculated. + + Note that the ``ctth`` dataset must contain satellite location + metadata, such as set in the ``orbital_parameters`` dataset attribute + that is set by many Satpy readers. It is essential that the datasets to be + corrected are coming from the same platform as the provided cloud top + height. + + A note on the algorithm and the implementation. Parallax correction + is inherently an inverse problem. The reported geolocation in + satellite data files is the true location plus the parallax error. + Therefore, this class first calculates the true geolocation (using + :func:`get_parallax_corrected_lonlats`), which gives a shifted longitude and + shifted latitude on an irregular grid. The difference between + the original and the shifted grid is the parallax error or shift. + The magnitude of this error can be estimated with + :func:`get_surface_parallax_displacement`. + With this difference, we need to invert the parallax correction to + calculate the corrected geolocation. Due to parallax correction, + high clouds shift a lot, low clouds shift a little, and cloud-free + pixels shift not at all. The shift may result in zero, one, + two, or more source pixel onto a destination pixel. Physically, + this corresponds to the situation where a narrow but high cloud is + viewed at a large angle. The cloud may occupy two or more pixels when + viewed at a large angle, but only one when viewed straight from above. + To accurately reproduce this perspective, the parallax correction uses + the :class:`~pyresample.bucket.BucketResampler` class, specifically + the :meth:`~pyresample.bucket.BucketResampler.get_abs_max` method, to + retain only the largest absolute shift (corresponding to the highest + cloud) within each pixel. Any other resampling method at this step + would yield incorrect results. When cloud moves over clear-sky, the + clear-sky pixel is unshifted and the shift is located exactly in the + centre of the grid box, so nearest-neighbour resampling would lead to + such shifts being deselected. Other resampling methods would average + large shifts with small shifts, leading to unpredictable results. + Now the reprojected shifts can be applied to the original lat/lon, + returning a new :class:`~pyresample.geometry.SwathDefinition`. + This is is the object returned by :meth:`corrected_area`. + + This procedure can be configured as a modifier using the + :class:`ParallaxCorrectionModifier` class. However, the modifier can only + be applied to one dataset at the time, which may not provide optimal + performance, although dask should reuse identical calculations between + multiple channels. + + """ + + def __init__(self, base_area, + debug_mode=False): + """Initialise parallax correction class. + + Args: + base_area (:class:`~pyresample.AreaDefinition`): Area for which calculated + geolocation will be calculated. + debug_mode (bool): Store diagnostic information in + self.diagnostics. This attribute always apply to the most + recently applied operation only. + """ + self.base_area = base_area + self.debug_mode = debug_mode + self.diagnostics = {} + + def __call__(self, cth_dataset, **kwargs): + """Apply parallax correction to dataset. + + Args: + cth_dataset: Dataset containing cloud top heights (or other heights + to be corrected). + + Returns: + :class:'~pyresample.geometry.SwathDefinition`: Swathdefinition with corrected + lat/lon coordinates. + """ + self.diagnostics.clear() + return self.corrected_area(cth_dataset, **kwargs) + + def corrected_area(self, cth_dataset, + cth_resampler="nearest", + cth_radius_of_influence=50000, + lonlat_chunks=1024): + """Return the parallax corrected SwathDefinition. + + Using the cloud top heights provided in ``cth_dataset``, calculate the + :class:`pyresample.geometry.SwathDefinition` that estimates the + geolocation for each pixel if it had been viewed from straight above + (without parallax error). The cloud top height will first be resampled + onto the area passed upon class initialisation in :meth:`__init__`. + Pixels that are invisible after parallax correction are not retained + but get geolocation NaN. + + Args: + cth_dataset (:class:`~xarray.DataArray`): Cloud top height in + meters. The variable attributes must contain an ``area`` + attribute describing the geolocation in a pyresample-aware way, + and they must contain satellite orbital parameters. The + dimensions must be ``(y, x)``. For best performance, this + should be a dask-based :class:`~xarray.DataArray`. + cth_resampler (string, optional): Resampler to use when resampling the + (cloud top) height to the base area. Defaults to "nearest". + cth_radius_of_influence (number, optional): Radius of influence to use when + resampling the (cloud top) height to the base area. Defaults + to 50000. + lonlat_chunks (int, optional): Chunking to use when calculating lon/lats. + Probably the default (1024) should be fine. + + Returns: + :class:`~pyresample.geometry.SwathDefinition` describing parallax + corrected geolocation. + """ + logger.debug("Calculating parallax correction using heights from " + f"{cth_dataset.attrs.get('name', cth_dataset.name)!s}, " + f"with base area {self.base_area.name!s}.") + (sat_lon, sat_lat, sat_alt_m) = _get_satpos_from_cth(cth_dataset) + self._check_overlap(cth_dataset) + + cth_dataset = self._prepare_cth_dataset( + cth_dataset, resampler=cth_resampler, + radius_of_influence=cth_radius_of_influence, + lonlat_chunks=lonlat_chunks) + + (base_lon, base_lat) = self.base_area.get_lonlats(chunks=lonlat_chunks) + # calculate the shift/error due to the parallax effect + (corrected_lon, corrected_lat) = get_parallax_corrected_lonlats( + sat_lon, sat_lat, sat_alt_m, + base_lon, base_lat, cth_dataset.data) + + shifted_area = self._get_swathdef_from_lon_lat(corrected_lon, corrected_lat) + + # But we are not actually moving pixels, rather we want a + # coordinate transformation. With this transformation we approximately + # invert the pixel coordinate transformation, giving the lon and lat + # where we should retrieve a value for a given pixel. + (proj_lon, proj_lat) = self._get_corrected_lon_lat( + base_lon, base_lat, shifted_area) + + return self._get_swathdef_from_lon_lat(proj_lon, proj_lat) + + @staticmethod + def _get_swathdef_from_lon_lat(lon, lat): + """Return a SwathDefinition from lon/lat. + + Turn ndarrays describing lon/lat into xarray with dimensions y, x, then + use these to create a :class:`~pyresample.geometry.SwathDefinition`. + """ + # lons and lats passed to SwathDefinition must be data-arrays with + # dimensions, see https://github.com/pytroll/satpy/issues/1434 + # and https://github.com/pytroll/satpy/issues/1997 + return SwathDefinition( + xr.DataArray(lon, dims=("y", "x")), + xr.DataArray(lat, dims=("y", "x"))) + + def _prepare_cth_dataset( + self, cth_dataset, resampler="nearest", radius_of_influence=50000, + lonlat_chunks=1024): + """Prepare CTH dataset. + + Set cloud top height to zero wherever lat/lon are valid but CTH is + undefined. Then resample onto the base area. + """ + # for calculating the parallax effect, set cth to 0 where it is + # undefined, unless pixels have no valid lat/lon + # NB: 0 may be below the surface... could be a problem for high + # resolution imagery in mountainous or high elevation terrain + # NB: how tolerant of xarray & dask is this? + resampled_cth_dataset = resample_dataset( + cth_dataset, self.base_area, resampler=resampler, + radius_of_influence=radius_of_influence) + (pixel_lon, pixel_lat) = resampled_cth_dataset.attrs["area"].get_lonlats( + chunks=lonlat_chunks) + masked_resampled_cth_dataset = resampled_cth_dataset.where( + np.isfinite(pixel_lon) & np.isfinite(pixel_lat)) + masked_resampled_cth_dataset = masked_resampled_cth_dataset.where( + masked_resampled_cth_dataset.notnull(), 0) + return masked_resampled_cth_dataset + + def _check_overlap(self, cth_dataset): + """Ensure cth_dataset is usable for parallax correction. + + Checks the coverage of ``cth_dataset`` compared to the ``base_area``. If + the entirety of ``base_area`` is covered by ``cth_dataset``, do + nothing. If only part of ``base_area`` is covered by ``cth_dataset``, + raise a `IncompleteHeightWarning`. If none of ``base_area`` is covered + by ``cth_dataset``, raise a `MissingHeightError`. + """ + warnings.warn( + "Overlap checking not impelemented. Waiting for " + "fix for https://github.com/pytroll/pyresample/issues/329") + + def _get_corrected_lon_lat(self, base_lon, base_lat, shifted_area): + """Calculate the corrected lon/lat based from the shifted area. + + After calculating the shifted area based on + :func:`get_parallax_corrected_lonlats`, + we invert the parallax error and estimate where those pixels came from. + For details on the algorithm, see the class docstring. + """ + (corrected_lon, corrected_lat) = shifted_area.get_lonlats(chunks=1024) + lon_diff = corrected_lon - base_lon + lat_diff = corrected_lat - base_lat + # We use the bucket resampler here, because parallax correction + # inevitably means there will be 2 source pixels ending up in the same + # destination pixel. We want to choose the biggest shift (max abs in + # lat_diff and lon_diff), because the biggest shift corresponds to the + # highest clouds, and if we move a 10 km cloud over a 2 km one, we + # should retain the 10 km. + # + # some things to keep in mind: + # - even with a constant cloud height, 3 source pixels may end up in + # the same destination pixel, because pixels get larger in the + # direction of the satellite. This means clouds may shrink as they + # approach the satellite. + # - the x-shift is a function of y and the y-shift is a function of x, + # so a cloud that was rectangular at the start may no longer be + # rectangular at the end + bur = BucketResampler(self.base_area, + da.array(corrected_lon), da.array(corrected_lat)) + inv_lat_diff = bur.get_abs_max(lat_diff) + inv_lon_diff = bur.get_abs_max(lon_diff) + + inv_lon = base_lon - inv_lon_diff + inv_lat = base_lat - inv_lat_diff + if self.debug_mode: + self.diagnostics["corrected_lon"] = corrected_lon + self.diagnostics["corrected_lat"] = corrected_lat + self.diagnostics["inv_lon"] = inv_lon + self.diagnostics["inv_lat"] = inv_lat + self.diagnostics["inv_lon_diff"] = inv_lon_diff + self.diagnostics["inv_lat_diff"] = inv_lat_diff + self.diagnostics["base_lon"] = base_lon + self.diagnostics["base_lat"] = base_lat + self.diagnostics["lon_diff"] = lon_diff + self.diagnostics["lat_diff"] = lat_diff + self.diagnostics["shifted_area"] = shifted_area + self.diagnostics["count"] = xr.DataArray( + bur.get_count(), dims=("y", "x"), attrs={"area": self.base_area}) + return (inv_lon, inv_lat) + + +class ParallaxCorrectionModifier(ModifierBase): + """Modifier for parallax correction. + + Apply parallax correction as a modifier. Uses the + :class:`ParallaxCorrection` class, which in turn uses the + :func:`get_parallax_corrected_lonlats` function. See the documentation there for + details on the behaviour. + + To use this, add to ``composites/visir.yaml`` within ``SATPY_CONFIG_PATH`` + something like:: + + sensor_name: visir + + modifiers: + parallax_corrected: + modifier: !!python/name:satpy.modifiers.parallax.ParallaxCorrectionModifier + prerequisites: + - "ctth_alti" + dataset_radius_of_influence: 50000 + + composites: + + parallax_corrected_VIS006: + compositor: !!python/name:satpy.composites.SingleBandCompositor + prerequisites: + - name: VIS006 + modifiers: [parallax_corrected] + + Here, ``ctth_alti`` is CTH provided by the ``nwcsaf-geo`` reader, so to use it + one would have to pass both on scene creation:: + + sc = Scene({"seviri_l1b_hrit": files_l1b, "nwcsaf-geo": files_l2}) + sc.load(["parallax_corrected_VIS006"]) + + The modifier takes optional global parameters, all of which are optional. + They affect various steps in the algorithm. Setting them may impact + performance:: + + cth_resampler + Resampler to use when resampling (cloud top) height to the base area. + Defaults to "nearest". + cth_radius_of_influence + Radius of influence to use when resampling the (cloud top) height to + the base area. Defaults to 50000. + lonlat_chunks + Chunk size to use when obtaining longitudes and latitudes from the area + definition. Defaults to 1024. If you set this to None, then parallax + correction will involve premature calculation. Changing this may or + may not make parallax correction slower or faster. + dataset_radius_of_influence + Radius of influence to use when resampling the dataset onto the + swathdefinition describing the parallax-corrected area. Defaults to + 50000. This always uses nearest neighbour resampling. + + Alternately, you can use the lower-level API directly with the + :class:`ParallaxCorrection` class, which may be more efficient if multiple + datasets need to be corrected. RGB Composites cannot be modified in this way + (i.e. you can't replace "VIS006" by "natural_color"). To get a parallax + corrected RGB composite, create a new composite where each input has the + modifier applied. The parallax calculation should only occur once, because + calculations are happening via dask and dask should reuse the calculation. + """ + + def __call__(self, projectables, optional_datasets=None, **info): + """Apply parallax correction. + + The argument ``projectables`` needs to contain the dataset to be + projected and the height to use for the correction. + """ + (to_be_corrected, cth) = projectables + base_area = to_be_corrected.attrs["area"] + corrector = self._get_corrector(base_area) + plax_corr_area = corrector( + cth, + cth_resampler=self.attrs.get("cth_resampler", "nearest"), + cth_radius_of_influence=self.attrs.get("cth_radius_of_influence", 50_000), + lonlat_chunks=self.attrs.get("lonlat_chunks", 1024), + ) + res = resample_dataset( + to_be_corrected, plax_corr_area, + radius_of_influence=self.attrs.get("dataset_radius_of_influence", 50_000), + fill_value=np.nan) + res.attrs["area"] = to_be_corrected.attrs["area"] + self.apply_modifier_info(to_be_corrected, res) + + return res + + def _get_corrector(self, base_area): + # only pass on those attributes that are arguments by + # ParallaxCorrection.__init__ + sig = inspect.signature(ParallaxCorrection.__init__) + kwargs = {} + for k in sig.parameters.keys() & self.attrs.keys(): + kwargs[k] = self.attrs[k] + corrector = ParallaxCorrection(base_area, **kwargs) + return corrector + + +def _get_satpos_from_cth(cth_dataset): + """Obtain satellite position from CTH dataset, height in meter. + + From a CTH dataset, obtain the satellite position lon, lat, altitude/m, + either directly from orbital parameters, or, when missing, from the + platform name using pyorbital and skyfield. + """ + (sat_lon, sat_lat, sat_alt_km) = get_satpos( + cth_dataset, use_tle=True) + return (sat_lon, sat_lat, sat_alt_km * 1000) diff --git a/satpy/readers/gpm_imerg.py b/satpy/readers/gpm_imerg.py index a511e404b8..3a68f8a9bb 100644 --- a/satpy/readers/gpm_imerg.py +++ b/satpy/readers/gpm_imerg.py @@ -83,9 +83,8 @@ def get_dataset(self, dataset_id, ds_info): val = data.attrs[key] if isinstance(val, h5py.h5r.Reference): del data.attrs[key] - if isinstance(val, np.ndarray): - if isinstance(val[0][0], h5py.h5r.Reference): - del data.attrs[key] + if isinstance(val, np.ndarray) and isinstance(val[0][0], h5py.h5r.Reference): + del data.attrs[key] return data def get_area_def(self, dsid): diff --git a/satpy/tests/modifier_tests/test_parallax.py b/satpy/tests/modifier_tests/test_parallax.py new file mode 100644 index 0000000000..70d51b49e9 --- /dev/null +++ b/satpy/tests/modifier_tests/test_parallax.py @@ -0,0 +1,798 @@ +# Copyright (c) 2021 Satpy developers +# +# This file is part of satpy. +# +# satpy is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# satpy is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +# A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +"""Tests related to parallax correction.""" + +import datetime +import logging +import math +import os +import unittest.mock + +import dask.array as da +import dask.config +import numpy as np +import pyorbital.tlefile +import pyresample.kd_tree +import pytest +import xarray as xr +from pyproj import Geod +from pyresample import create_area_def + +import satpy.resample + +from ...writers import get_enhanced_image + + +@pytest.fixture +def fake_tle(): + """Produce fake Two Line Element (TLE) object from pyorbital.""" + return pyorbital.tlefile.Tle( + "Meteosat-42", + line1="1 40732U 15034A 22011.84285506 .00000004 00000+0 00000+0 0 9995", + line2="2 40732 0.2533 325.0106 0000976 118.8734 330.4058 1.00272123 23817") + + +def _get_fake_areas(center, sizes, resolution, code=4326): + """Get multiple square areas with the same center. + + Returns multiple square areas centered at the same location + + Args: + center (Tuple[float, float]): Center of all areass + sizes (List[int]): Sizes of areas + resolution (float): Resolution of fake area. + + Returns: + List of areas. + """ + return [create_area_def( + "fribullus_xax", + code, + units="degrees", + resolution=resolution, + center=center, + shape=(size, size)) + for size in sizes] + + +def _get_attrs(lat, lon, height=35_000): + """Get attributes for datasets in fake scene.""" + return { + "orbital_parameters": { + "satellite_actual_altitude": height, # in km above surface + "satellite_actual_longitude": lon, + "satellite_actual_latitude": lat}, + "units": "m" # does not apply to orbital parameters, I think! + } + + +class TestForwardParallax: + """Test the forward parallax function with various inputs.""" + + def test_get_parallax_corrected_lonlats_ssp(self): + """Test that at SSP, parallax correction does nothing.""" + from ...modifiers.parallax import get_parallax_corrected_lonlats + sat_lat = sat_lon = lon = lat = 0. + height = 5000. # m + sat_alt = 30_000_000. # m + corr_lon, corr_lat = get_parallax_corrected_lonlats( + sat_lon, sat_lat, sat_alt, lon, lat, height) + assert corr_lon == corr_lat == 0 + + def test_get_parallax_corrected_lonlats_clearsky(self): + """Test parallax correction for clearsky case (returns NaN).""" + from ...modifiers.parallax import get_parallax_corrected_lonlats + sat_lat = sat_lon = 0 + lat = np.linspace(-20, 20, 25).reshape(5, 5) + lon = np.linspace(-20, 20, 25).reshape(5, 5).T + height = np.full((5, 5), np.nan) # no CTH --> clearsky + sat_alt = 35_000_000. # m above surface + (corr_lon, corr_lat) = get_parallax_corrected_lonlats( + sat_lon, sat_lat, sat_alt, lon, lat, height) + # clearsky becomes NaN + assert np.isnan(corr_lon).all() + assert np.isnan(corr_lat).all() + + @pytest.mark.parametrize("lat,lon", [(0, 0), (0, 40), (0, 179.9)]) + @pytest.mark.parametrize("resolution", [0.01, 0.5, 10]) + def test_get_parallax_corrected_lonlats_cloudy_ssp(self, lat, lon, resolution): + """Test parallax correction for fully cloudy scene at SSP.""" + from ...modifiers.parallax import get_parallax_corrected_lonlats + + N = 5 + lats = np.linspace(lat-N*resolution, lat+N*resolution, 25).reshape(N, N) + lons = np.linspace(lon-N*resolution, lon+N*resolution, 25).reshape(N, N).T + height = np.full((N, N), 10_000) # constant high clouds at 10 km + sat_alt = 35_000_000. # satellite at 35 Mm + (corr_lon, corr_lat) = get_parallax_corrected_lonlats( + lon, lat, sat_alt, lons, lats, height) + # confirm movements behave as expected + geod = Geod(ellps="sphere") + # need to use np.tile here as geod.inv doesn't seem to broadcast (not + # when turning lon/lat in arrays of size (1, 1) either) + corr_dist = geod.inv(np.tile(lon, [N, N]), np.tile(lat, [N, N]), corr_lon, corr_lat)[2] + corr_delta = geod.inv(corr_lon, corr_lat, lons, lats)[2] + uncorr_dist = geod.inv(np.tile(lon, [N, N]), np.tile(lat, [N, N]), lons, lats)[2] + # should be equal at SSP and nowhere else + np.testing.assert_allclose(corr_delta[2, 2], 0, atol=1e-9) + assert np.isclose(corr_delta, 0, atol=1e-9).sum() == 1 + # should always get closer to SSP + assert (uncorr_dist - corr_dist >= -1e-8).all() + # should be larger the further we get from SSP + assert (np.diff(corr_delta[N//2, :N//2+1]) < 0).all() + assert (np.diff(corr_delta[N//2, N//2:]) > 0).all() + assert (np.diff(corr_delta[N//2:, N//2]) > 0).all() + assert (np.diff(corr_delta[:N//2+1, N//2]) < 0).all() + assert (np.diff(np.diag(corr_delta)[:N//2+1]) < 0).all() + assert (np.diff(np.diag(corr_delta)[N//2:]) > 0).all() + + def test_get_parallax_corrected_lonlats_cloudy_slant(self): + """Test parallax correction for fully cloudy scene (not SSP).""" + from ...modifiers.parallax import get_parallax_corrected_lonlats + sat_lat = sat_lon = 0 + lat = np.linspace(-20, 20, 25).reshape(5, 5) + lon = np.linspace(-20, 20, 25).reshape(5, 5).T + height = np.full((5, 5), 10_000) # constant high clouds at 10 km + sat_alt = 35_000_000. # satellite at 35 Mm + (corr_lon, corr_lat) = get_parallax_corrected_lonlats( + sat_lon, sat_lat, sat_alt, lon, lat, height) + # reference value from Simon Proud + np.testing.assert_allclose( + corr_lat[4, 4], 19.955, rtol=5e-4) + np.testing.assert_allclose( + corr_lon[4, 4], 19.960, rtol=5e-4) + + def test_get_parallax_corrected_lonlats_mixed(self): + """Test parallax correction for mixed cloudy case.""" + from ...modifiers.parallax import get_parallax_corrected_lonlats + + sat_lon = sat_lat = 0 + sat_alt = 35_785_831.0 # m + lon = da.array([[-20, -10, 0, 10, 20]]*5) + lat = da.array([[-20, -10, 0, 10, 20]]*5).T + alt = da.array([ + [np.nan, np.nan, 5000., 6000., np.nan], + [np.nan, 6000., 7000., 7000., 7000.], + [np.nan, 7000., 8000., 9000., np.nan], + [np.nan, 7000., 7000., 7000., np.nan], + [np.nan, 4000., 3000., np.nan, np.nan]]) + (corrected_lon, corrected_lat) = get_parallax_corrected_lonlats( + sat_lon, sat_lat, sat_alt, lon, lat, alt) + assert corrected_lon.shape == lon.shape + assert corrected_lat.shape == lat.shape + # lon/lat should be nan for clear-sky pixels + assert np.isnan(corrected_lon[np.isnan(alt)]).all() + assert np.isnan(corrected_lat[np.isnan(alt)]).all() + # otherwise no nans + assert np.isfinite(corrected_lon[~np.isnan(alt)]).all() + assert np.isfinite(corrected_lat[~np.isnan(alt)]).all() + + def test_get_parallax_corrected_lonlats_horizon(self): + """Test that exception is raised if satellites exactly at the horizon. + + Test the rather unlikely case of a satellite elevation of exactly 0 + """ + from ...modifiers.parallax import get_parallax_corrected_lonlats + sat_lat = sat_lon = lon = lat = 0. + height = 5000. + sat_alt = 30_000_000. + with unittest.mock.patch("satpy.modifiers.parallax.get_observer_look") as smpg: + smpg.return_value = (0, 0) + with pytest.raises(NotImplementedError): + get_parallax_corrected_lonlats(sat_lon, sat_lat, sat_alt, lon, lat, height) + + def test_get_surface_parallax_displacement(self): + """Test surface parallax displacement.""" + from ...modifiers.parallax import get_surface_parallax_displacement + + val = get_surface_parallax_displacement( + 0, 0, 36_000_000, 0, 10, 10_000) + np.testing.assert_allclose(val, 2141.2404451757875) + + +class TestParallaxCorrectionClass: + """Test that the ParallaxCorrection class is behaving sensibly.""" + + @pytest.mark.parametrize("center", [(0, 0), (80, -10), (-180, 5)]) + @pytest.mark.parametrize("sizes", [[5, 9]]) + @pytest.mark.parametrize("resolution", [0.05, 1, 10]) + def test_init_parallaxcorrection(self, center, sizes, resolution): + """Test that ParallaxCorrection class can be instantiated.""" + from ...modifiers.parallax import ParallaxCorrection + fake_area = _get_fake_areas(center, sizes, resolution)[0] + pc = ParallaxCorrection(fake_area) + assert pc.base_area == fake_area + + @pytest.mark.parametrize("sat_pos,ar_pos", + [((0, 0), (0, 0)), ((0, 0), (40, 0))]) + @pytest.mark.parametrize("resolution", [0.01, 0.5, 10]) + def test_correct_area_clearsky(self, sat_pos, ar_pos, resolution, caplog): + """Test that ParallaxCorrection doesn't change clearsky geolocation.""" + from ...modifiers.parallax import ParallaxCorrection + from ..utils import make_fake_scene + (sat_lat, sat_lon) = sat_pos + (ar_lat, ar_lon) = ar_pos + small = 5 + large = 9 + (fake_area_small, fake_area_large) = _get_fake_areas( + (ar_lon, ar_lat), [small, large], resolution) + corrector = ParallaxCorrection(fake_area_small) + + sc = make_fake_scene( + {"CTH_clear": np.full((large, large), np.nan)}, + daskify=False, + area=fake_area_large, + common_attrs=_get_attrs(sat_lat, sat_lon, 35_000)) + + with caplog.at_level(logging.DEBUG): + new_area = corrector(sc["CTH_clear"]) + assert "Calculating parallax correction using heights from CTH_clear" in caplog.text + np.testing.assert_allclose( + new_area.get_lonlats(), + fake_area_small.get_lonlats()) + + @pytest.mark.parametrize("lat,lon", + [(0, 0), (0, 40), (0, 180), + (90, 0)]) # relevant for Арктика satellites + @pytest.mark.parametrize("resolution", [0.01, 0.5, 10]) + def test_correct_area_ssp(self, lat, lon, resolution): + """Test that ParallaxCorrection doesn't touch SSP.""" + from ...modifiers.parallax import ParallaxCorrection + from ..utils import make_fake_scene + codes = { + (0, 0): 4326, + (0, 40): 4326, + (0, 180): 3575, + (90, 0): 3575} + small = 5 + large = 9 + (fake_area_small, fake_area_large) = _get_fake_areas( + (lon, lat), [small, large], resolution, + code=codes[(lat, lon)]) + corrector = ParallaxCorrection(fake_area_small) + + sc = make_fake_scene( + {"CTH_constant": np.full((large, large), 10000)}, + daskify=False, + area=fake_area_large, + common_attrs=_get_attrs(lat, lon, 35_000)) + new_area = corrector(sc["CTH_constant"]) + assert new_area.shape == fake_area_small.shape + old_lonlats = fake_area_small.get_lonlats() + new_lonlats = new_area.get_lonlats() + if lat != 90: # don't check SSP longitude if lat=90 + np.testing.assert_allclose( + old_lonlats[0][2, 2], + new_lonlats[0][2, 2], + atol=1e-9) + np.testing.assert_allclose( + old_lonlats[0][2, 2], + lon, + atol=1e-9) + np.testing.assert_allclose( + old_lonlats[1][2, 2], + new_lonlats[1][2, 2], + atol=1e-9) + np.testing.assert_allclose( + old_lonlats[1][2, 2], + lat, + atol=1e-9) + + @pytest.mark.parametrize("daskify", [False, True]) + def test_correct_area_partlycloudy(self, daskify): + """Test ParallaxCorrection for partly cloudy situation.""" + from ...modifiers.parallax import ParallaxCorrection + from ..utils import make_fake_scene + small = 5 + large = 9 + (fake_area_small, fake_area_large) = _get_fake_areas( + (0, 50), [small, large], 0.1) + (fake_area_lons, fake_area_lats) = fake_area_small.get_lonlats() + corrector = ParallaxCorrection(fake_area_small) + + sc = make_fake_scene( + {"CTH": np.array([ + [np.nan, np.nan, 5000., 6000., 7000., 6000., 5000., np.nan, np.nan], + [np.nan, 6000., 7000., 7000., 7000., np.nan, np.nan, np.nan, np.nan], + [np.nan, 7000., 8000., 9000., np.nan, np.nan, np.nan, np.nan, np.nan], + [np.nan, 7000., 7000., 7000., np.nan, np.nan, np.nan, np.nan, np.nan], + [np.nan, 4000., 3000., np.nan, np.nan, np.nan, np.nan, np.nan, np.nan], + [np.nan, np.nan, 5000., 8000., 8000., 8000., 6000., np.nan, np.nan], + [np.nan, 9000., 9000., 9000., 9000., 9000., 9000., 9000., np.nan], + [np.nan, 9000., 9000., 9000., 9000., 9000., 9000., 9000., np.nan], + [np.nan, 9000., 9000., 9000., 9000., 9000., 9000., 9000., np.nan], + ])}, + daskify=daskify, + area=fake_area_large, + common_attrs=_get_attrs(0, 0, 40_000)) + new_area = corrector(sc["CTH"]) + assert new_area.shape == fake_area_small.shape + (new_lons, new_lats) = new_area.get_lonlats() + assert fake_area_lons[3, 4] != new_lons[3, 4] + + np.testing.assert_allclose( + new_lons, + np.array([ + [np.nan, np.nan, 0.0, 0.1, 0.2], + [-0.20078652, -0.10044222, 0.0, 0.1, 0.2], + [-0.20068529, -0.10034264, 0.0, 0.1, 0.2], + [np.nan, np.nan, np.nan, np.nan, np.nan], + [-0.20048537, -0.10038778, 0., 0.10038778, 0.20058219]]), + rtol=1e-5) + np.testing.assert_allclose( + new_lats, + np.array([ + [np.nan, np.nan, 50.2, 50.2, 50.2], + [50.2110675, 50.22493181, 50.1, 50.1, 50.1], + [50.09680357, 50.09680346, 50.0, 50.0, 50.0], + [np.nan, np.nan, np.nan, np.nan, np.nan], + [49.86860622, 49.9097198, 49.90971976, 49.9097198, 49.88231496]]), + rtol=1e-6) + + @pytest.mark.parametrize("res1,res2", [(0.08, 0.3), (0.3, 0.08)]) + def test_correct_area_clearsky_different_resolutions(self, res1, res2): + """Test clearsky correction when areas have different resolutions.""" + from ...modifiers.parallax import ParallaxCorrection + from ..utils import make_fake_scene + + # areas with different resolutions, but same coverage + + area1 = create_area_def( + "fribullus_xax", + 4326, + units="degrees", + resolution=res1, + area_extent=[-1, -1, 1, 1]) + + area2 = create_area_def( + "fribullus_xax", + 4326, + units="degrees", + resolution=res2, + area_extent=[-1, -1, 1, 1]) + + with pytest.warns(None) as record: + sc = make_fake_scene( + {"CTH_clear": np.full(area1.shape, np.nan)}, + daskify=False, + area=area1, + common_attrs=_get_attrs(0, 0, 35_000)) + assert len(record) == 0 + + corrector = ParallaxCorrection(area2) + new_area = corrector(sc["CTH_clear"]) + np.testing.assert_allclose( + new_area.get_lonlats(), + area2.get_lonlats()) + + @pytest.mark.xfail(reason="awaiting pyresample fixes") + def test_correct_area_cloudy_no_overlap(self, ): + """Test cloudy correction when areas have no overlap.""" + from ...modifiers.parallax import MissingHeightError, ParallaxCorrection + from ..utils import make_fake_scene + areas_00 = _get_fake_areas((0, 40), [5, 9], 0.1) + areas_shift = _get_fake_areas((90, 20), [5, 9], 0.1) + fake_area_small = areas_00[0] + fake_area_large = areas_shift[1] + + sc = make_fake_scene( + {"CTH_constant": np.full((9, 9), 10000)}, + daskify=False, + area=fake_area_large, + common_attrs=_get_attrs(0, 0, 35_000)) + + corrector = ParallaxCorrection(fake_area_small) + with pytest.raises(MissingHeightError): + corrector(sc["CTH_constant"]) + + @pytest.mark.xfail(reason="awaiting pyresample fixes") + def test_correct_area_cloudy_partly_shifted(self, ): + """Test cloudy correction when areas overlap only partly.""" + from ...modifiers.parallax import IncompleteHeightWarning, ParallaxCorrection + from ..utils import make_fake_scene + areas_00 = _get_fake_areas((0, 40), [5, 9], 0.1) + areas_shift = _get_fake_areas((0.5, 40), [5, 9], 0.1) + fake_area_small = areas_00[0] + fake_area_large = areas_shift[1] + + sc = make_fake_scene( + {"CTH_constant": np.full((9, 9), 10000)}, + daskify=False, + area=fake_area_large, + common_attrs=_get_attrs(0, 0, 35_000)) + + corrector = ParallaxCorrection(fake_area_small) + + with pytest.warns(IncompleteHeightWarning): + new_area = corrector(sc["CTH_constant"]) + assert new_area.shape == fake_area_small.shape + + def test_correct_area_cloudy_same_area(self, ): + """Test cloudy correction when areas are the same.""" + from ...modifiers.parallax import ParallaxCorrection + from ..utils import make_fake_scene + area = _get_fake_areas((0, 0), [9], 0.1)[0] + + sc = make_fake_scene( + {"CTH_constant": np.full((9, 9), 10000)}, + daskify=False, + area=area, + common_attrs=_get_attrs(0, 0, 35_000)) + + corrector = ParallaxCorrection(area) + corrector(sc["CTH_constant"]) + + def test_correct_area_no_orbital_parameters(self, caplog, fake_tle): + """Test ParallaxCorrection when CTH has no orbital parameters. + + Some CTH products, such as NWCSAF-GEO, do not include information + on satellite location directly. Rather, they include platform name, + sensor, start time, and end time, that we have to use instead. + """ + from ...modifiers.parallax import ParallaxCorrection + from ..utils import make_fake_scene + small = 5 + large = 9 + (fake_area_small, fake_area_large) = _get_fake_areas( + (0, 0), [small, large], 0.05) + corrector = ParallaxCorrection(fake_area_small) + + sc = make_fake_scene( + {"CTH_clear": np.full((large, large), np.nan)}, + daskify=False, + area=fake_area_large, + common_attrs={ + "platform_name": "Meteosat-42", + "sensor": "irives", + "start_time": datetime.datetime(3021, 11, 30, 12, 24, 17), + "end_time": datetime.datetime(3021, 11, 30, 12, 27, 22)}) + with unittest.mock.patch("pyorbital.tlefile.read") as plr: + plr.return_value = fake_tle + with caplog.at_level(logging.WARNING): + new_area = corrector(sc["CTH_clear"]) + assert "Orbital parameters missing from metadata." in caplog.text + np.testing.assert_allclose( + new_area.get_lonlats(), + fake_area_small.get_lonlats()) + + +class TestParallaxCorrectionModifier: + """Test that the parallax correction modifier works correctly.""" + + def test_parallax_modifier_interface(self): + """Test the modifier interface.""" + from ...modifiers.parallax import ParallaxCorrectionModifier + (area_small, area_large) = _get_fake_areas((0, 0), [5, 9], 0.1) + fake_bt = xr.DataArray( + np.linspace(220, 230, 25).reshape(5, 5), + dims=("y", "x"), + attrs={"area": area_small, **_get_attrs(0, 0, 35_000)}) + cth_clear = xr.DataArray( + np.full((9, 9), np.nan), + dims=("y", "x"), + attrs={"area": area_large, **_get_attrs(0, 0, 35_000)}) + modif = ParallaxCorrectionModifier( + name="parallax_corrected_dataset", + prerequisites=[fake_bt, cth_clear], + optional_prerequisites=[], + cth_radius_of_influence=48_000, + dataset_radius_of_influence=49_000) + res = modif([fake_bt, cth_clear], optional_datasets=[]) + np.testing.assert_allclose(res, fake_bt) + with unittest.mock.patch("satpy.modifiers.parallax.resample_dataset") as smp: + smp.side_effect = satpy.resample.resample_dataset + modif([fake_bt, cth_clear], optional_datasets=[]) + assert smp.call_args_list[0].kwargs["radius_of_influence"] == 48_000 + assert smp.call_args_list[1].kwargs["radius_of_influence"] == 49_000 + + def test_parallax_modifier_interface_with_cloud(self): + """Test the modifier interface with a cloud. + + Test corresponds to a real bug encountered when using CTH data + from NWCSAF-GEO, which created strange speckles in Africa (see + https://github.com/pytroll/satpy/pull/1904#issuecomment-1011161623 + for an example). Create fake CTH corresponding to NWCSAF-GEO area and + BT corresponding to full disk SEVIRI, and test that no strange speckles + occur. + """ + from ...modifiers.parallax import ParallaxCorrectionModifier + + w_cth = 25 + h_cth = 15 + proj_dict = {'a': '6378137', 'h': '35785863', 'proj': 'geos', 'units': 'm'} + fake_area_cth = pyresample.create_area_def( + area_id="test-area", + projection=proj_dict, + area_extent=(-2296808.75, 2785874.75, 2293808.25, 5570249.0), + shape=(h_cth, w_cth)) + + sz_bt = 20 + fake_area_bt = pyresample.create_area_def( + "test-area-2", + projection=proj_dict, + area_extent=(-5567248.0742, -5513240.8172, 5513240.8172, 5567248.0742), + shape=(sz_bt, sz_bt)) + + (lons_cth, lats_cth) = fake_area_cth.get_lonlats() + fake_cth_data = np.where( + np.isfinite(lons_cth) & np.isfinite(lats_cth), + 15000, + np.nan) + + (lons_bt, lats_bt) = fake_area_bt.get_lonlats() + fake_bt_data = np.where( + np.isfinite(lons_bt) & np.isfinite(lats_bt), + np.linspace(200, 300, lons_bt.size).reshape(lons_bt.shape), + np.nan) + + attrs = _get_attrs(0, 0) + + fake_bt = xr.DataArray( + fake_bt_data, + dims=("y", "x"), + attrs={**attrs, "area": fake_area_bt}) + fake_cth = xr.DataArray( + fake_cth_data, + dims=("y", "x"), + attrs={**attrs, "area": fake_area_cth}) + + modif = ParallaxCorrectionModifier( + name="parallax_corrected_dataset", + prerequisites=[fake_bt, fake_cth], + optional_prerequisites=[], + search_radius=25_000) + + res = modif([fake_bt, fake_cth], optional_datasets=[]) + + # with a constant cloud, a monotonically increasing BT should still + # do so after parallax correction + assert not (res.diff("x") < 0).any() + + @pytest.fixture + def test_area(self, request): + """Produce test area for parallax correction unit tests. + + Produce test area for the modifier-interface parallax correction unit + tests. + """ + extents = { + "foroyar": [-861785.8867075047, 6820719.391005835, -686309.8124887547, 6954386.383193335], + "ouagadougou": [-232482.90622750926, 1328206.360136668, + -114074.70310250926, 1422810.852324168], + } + where = request.param + return pyresample.create_area_def(where, 4087, area_extent=extents[where], resolution=500) + + def _get_fake_cloud_datasets(self, test_area, cth, use_dask): + """Return datasets for BT and CTH for fake cloud.""" + w_cloud = 20 + h_cloud = 5 + + # location of cloud in uncorrected data + lat_min_i = 155 + lat_max_i = lat_min_i + h_cloud + lon_min_i = 140 + lon_max_i = lon_min_i + w_cloud + + fake_bt_data = np.linspace( + 270, 330, math.prod(test_area.shape), dtype="f8").reshape( + test_area.shape).round(2) + fake_cth_data = np.full(test_area.shape, np.nan, dtype="f8") + fake_bt_data[lat_min_i:lat_max_i, lon_min_i:lon_max_i] = np.linspace( + 180, 220, w_cloud*h_cloud).reshape(h_cloud, w_cloud).round(2) + fake_cth_data[lat_min_i:lat_max_i, lon_min_i:lon_max_i] = cth + + if use_dask: + fake_bt_data = da.array(fake_bt_data) + fake_cth_data = da.array(fake_cth_data) + + attrs = _get_attrs(0, 0) + + fake_bt = xr.DataArray( + fake_bt_data, + dims=("y", "x"), + attrs={**attrs, "area": test_area}) + + fake_cth = xr.DataArray( + fake_cth_data, + dims=("y", "x"), + attrs={**attrs, "area": test_area}) + + cma = np.zeros(shape=fake_bt.shape, dtype="?") + cma[lat_min_i:lat_max_i, lon_min_i:lon_max_i] = True + + return (fake_bt, fake_cth, cma) + + @pytest.mark.parametrize("test_area", ["foroyar", "ouagadougou"], indirect=["test_area"]) + def test_modifier_interface_fog_no_shift(self, test_area): + """Test that fog isn't masked or shifted.""" + from ...modifiers.parallax import ParallaxCorrectionModifier + + (fake_bt, fake_cth, _) = self._get_fake_cloud_datasets(test_area, 50, use_dask=False) + + modif = ParallaxCorrectionModifier( + name="parallax_corrected_dataset", + prerequisites=[fake_bt, fake_cth], + optional_prerequisites=[], + debug_mode=True) + + res = modif([fake_bt, fake_cth], optional_datasets=[]) + + assert np.isfinite(res).all() + np.testing.assert_allclose(res, fake_bt) + + @pytest.mark.parametrize("cth", [7500, 15000]) + @pytest.mark.parametrize("use_dask", [True, False]) + @pytest.mark.parametrize("test_area", ["foroyar", "ouagadougou"], indirect=["test_area"]) + def test_modifier_interface_cloud_moves_to_observer(self, cth, use_dask, test_area): + """Test that a cloud moves to the observer. + + With the modifier interface, use a high resolution area and test that + pixels are moved in the direction of the observer and not away from it. + """ + from ...modifiers.parallax import ParallaxCorrectionModifier + + (fake_bt, fake_cth, cma) = self._get_fake_cloud_datasets(test_area, cth, use_dask=use_dask) + + # location of cloud in corrected data + # this may no longer be rectangular! + dest_mask = np.zeros(shape=test_area.shape, dtype="?") + cloud_location = { + "foroyar": { + 7500: (197, 202, 152, 172), + 15000: (239, 244, 165, 184)}, + "ouagadougou": { + 7500: (159, 164, 140, 160), + 15000: (163, 168, 141, 161)}} + (x_lo, x_hi, y_lo, y_hi) = cloud_location[test_area.name][cth] + dest_mask[x_lo:x_hi, y_lo:y_hi] = True + + modif = ParallaxCorrectionModifier( + name="parallax_corrected_dataset", + prerequisites=[fake_bt, fake_cth], + optional_prerequisites=[], + debug_mode=True) + + res = modif([fake_bt, fake_cth], optional_datasets=[]) + + assert fake_bt.attrs["area"] == test_area # should not be changed + assert res.attrs["area"] == fake_bt.attrs["area"] + # confirm old cloud area now fill value + # except where it overlaps with new cloud + assert np.isnan(res.data[cma & (~dest_mask)]).all() + # confirm rest of the area does not have fill values + assert np.isfinite(res.data[~cma]).all() + # confirm that rest of area pixel values did not change, except where + # cloud arrived or originated + delta = res - fake_bt + assert (delta.data[~(cma | dest_mask)] == 0).all() + # verify that cloud moved south. Pointwise comparison might not work because + # cloud may shrink. + assert ((res.attrs["area"].get_lonlats()[1][dest_mask]).mean() < + fake_bt.attrs["area"].get_lonlats()[1][cma].mean()) + # verify that all pixels at the new cloud location are indeed cloudy + assert (res.data[dest_mask] < 250).all() + + +_test_yaml_code = """ +sensor_name: visir + +modifiers: + parallax_corrected: + modifier: !!python/name:satpy.modifiers.parallax.ParallaxCorrectionModifier + prerequisites: + - name: "ctth_alti" + +composites: + parallax_corrected_VIS006: + compositor: !!python/name:satpy.composites.SingleBandCompositor + prerequisites: + - name: VIS006 + modifiers: [parallax_corrected] +""" + + +class TestParallaxCorrectionSceneLoad: + """Test that scene load interface works as expected.""" + + @pytest.fixture + def yaml_code(self): + """Return YAML code for parallax_corrected_VIS006.""" + return _test_yaml_code + + @pytest.fixture + def conf_file(self, yaml_code, tmp_path): + """Produce a fake configuration file.""" + conf_file = tmp_path / "test.yaml" + with conf_file.open(mode="wt", encoding="ascii") as fp: + fp.write(yaml_code) + return conf_file + + @pytest.fixture + def fake_scene(self, yaml_code): + """Produce fake scene and prepare fake composite config.""" + from satpy import Scene + from satpy.dataset.dataid import WavelengthRange + from satpy.tests.utils import make_dataid + + area = _get_fake_areas((0, 0), [5], 1)[0] + sc = Scene() + sc["VIS006"] = xr.DataArray( + np.linspace(0, 99, 25).reshape(5, 5), + dims=("y", "x"), + attrs={ + "_satpy_id": make_dataid( + name="VIS006", + wavelength=WavelengthRange(min=0.56, central=0.635, max=0.71, unit="µm"), + resolution=3000, + calibration="reflectance", + modifiers=()), + "modifiers": (), + "sensor": "seviri", + "area": area}) + sc["ctth_alti"] = xr.DataArray( + np.linspace(0, 99, 25).reshape(5, 5), + dims=("y", "x"), + attrs={ + "_satpy_id": make_dataid( + name="ctth_alti", + resolution=3000, + modifiers=()), + "modifiers": (), + "sensor": {"seviri"}, + "platform_name": "Meteosat-11", + "start_time": datetime.datetime(2022, 4, 12, 9, 0), + "area": area}) + return sc + + def test_double_load(self, fake_scene, conf_file, fake_tle): + """Test that loading corrected and uncorrected works correctly. + + When the modifier ``__call__`` method fails to call + ``self.apply_modifier_info(new, old)`` and the original and + parallax-corrected dataset are requested at the same time, the + DataArrays differ but the underlying dask arrays have object identity, + which in turn leads to both being parallax corrected. This unit test + confirms that there is no such object identity. + """ + with unittest.mock.patch( + "satpy.composites.config_loader.config_search_paths") as sccc, \ + unittest.mock.patch("pyorbital.tlefile.read") as plr: + sccc.return_value = [os.fspath(conf_file)] + plr.return_value = fake_tle + fake_scene.load(["parallax_corrected_VIS006", "VIS006"]) + assert fake_scene["VIS006"] is not fake_scene["parallax_corrected_VIS006"] + assert fake_scene["VIS006"].data is not fake_scene["parallax_corrected_VIS006"].data + + @pytest.mark.xfail(reason="awaiting pyresample fixes") + def test_no_compute(self, fake_scene, conf_file): + """Test that no computation occurs.""" + from satpy.tests.utils import CustomScheduler + with unittest.mock.patch( + "satpy.composites.config_loader.config_search_paths") as sccc, \ + dask.config.set(scheduler=CustomScheduler(max_computes=0)): + sccc.return_value = [os.fspath(conf_file)] + fake_scene.load(["parallax_corrected_VIS006"]) + + def test_enhanced_image(self, fake_scene, conf_file, fake_tle): + """Test that image enhancement is the same.""" + with unittest.mock.patch( + "satpy.composites.config_loader.config_search_paths") as sccc, \ + unittest.mock.patch("pyorbital.tlefile.read") as plr: + sccc.return_value = [os.fspath(conf_file)] + plr.return_value = fake_tle + fake_scene.load(["parallax_corrected_VIS006", "VIS006"]) + im1 = get_enhanced_image(fake_scene["VIS006"]) + im2 = get_enhanced_image(fake_scene["parallax_corrected_VIS006"]) + assert im1.data.attrs["enhancement_history"] == im2.data.attrs["enhancement_history"] diff --git a/satpy/tests/test_utils.py b/satpy/tests/test_utils.py index b2e53ebd46..354c1886bb 100644 --- a/satpy/tests/test_utils.py +++ b/satpy/tests/test_utils.py @@ -18,6 +18,7 @@ """Testing of utils.""" from __future__ import annotations +import datetime import logging import typing import unittest @@ -271,6 +272,28 @@ def test_get_satpos_fails_with_informative_error(self, attrs): with pytest.raises(KeyError, match="Unable to determine satellite position.*"): get_satpos(data_arr) + def test_get_satpos_from_satname(self, caplog): + """Test getting satellite position from satellite name only.""" + import pyorbital.tlefile + + data_arr = xr.DataArray( + (), + attrs={ + "platform_name": "Meteosat-42", + "sensor": "irives", + "start_time": datetime.datetime(2031, 11, 20, 19, 18, 17)}) + with mock.patch("pyorbital.tlefile.read") as plr: + plr.return_value = pyorbital.tlefile.Tle( + "Meteosat-42", + line1="1 40732U 15034A 22011.84285506 .00000004 00000+0 00000+0 0 9995", + line2="2 40732 0.2533 325.0106 0000976 118.8734 330.4058 1.00272123 23817") + with caplog.at_level(logging.WARNING): + (lon, lat, alt) = get_satpos(data_arr, use_tle=True) + assert "Orbital parameters missing from metadata" in caplog.text + np.testing.assert_allclose( + (lon, lat, alt), + (119.39533705010592, -1.1491628298731498, 35803.19986408156)) + def test_make_fake_scene(): """Test the make_fake_scene utility. diff --git a/satpy/utils.py b/satpy/utils.py index 95383e6f6a..e1489a9c19 100644 --- a/satpy/utils.py +++ b/satpy/utils.py @@ -21,6 +21,7 @@ from __future__ import annotations import contextlib +import datetime import logging import os import warnings @@ -42,6 +43,8 @@ _is_logging_on = False TRACE_LEVEL = 5 +logger = logging.getLogger(__name__) + class PerformanceWarning(Warning): """Warning raised when there is a possible performance impact.""" @@ -183,7 +186,18 @@ def in_ipynb(): def lonlat2xyz(lon, lat): - """Convert lon lat to cartesian.""" + """Convert lon lat to cartesian. + + For a sphere with unit radius, convert the spherical coordinates + longitude and latitude to cartesian coordinates. + + Args: + lon (number or array of numbers): Longitude in °. + lat (number or array of numbers): Latitude in °. + + Returns: + (x, y, z) Cartesian coordinates [1] + """ lat = np.deg2rad(lat) lon = np.deg2rad(lon) x = np.cos(lat) * np.cos(lon) @@ -193,7 +207,21 @@ def lonlat2xyz(lon, lat): def xyz2lonlat(x, y, z, asin=False): - """Convert cartesian to lon lat.""" + """Convert cartesian to lon lat. + + For a sphere with unit radius, convert cartesian coordinates to spherical + coordinates longitude and latitude. + + Args: + x (number or array of numbers): x-coordinate, unitless + y (number or array of numbers): y-coordinate, unitless + z (number or array of numbers): z-coordinate, unitless + asin (optional, bool): If true, use arcsin for calculations. + If false, use arctan2 for calculations. + + Returns: + (lon, lat): Longitude and latitude in °. + """ lon = np.rad2deg(np.arctan2(y, x)) if asin: lat = np.rad2deg(np.arcsin(z)) @@ -291,6 +319,7 @@ def atmospheric_path_length_correction(data, cos_zen, limit=88., max_sza=95.): def get_satpos( data_arr: xr.DataArray, preference: Optional[str] = None, + use_tle: bool = False ) -> tuple[float, float, float]: """Get satellite position from dataset attributes. @@ -309,9 +338,14 @@ def get_satpos( preference is not available then the original preference list is used. A warning is issued when projection values have to be used because nothing else is available and it wasn't provided as the ``preference``. + use_tle: If true, try to obtain position via satellite name + and TLE if it can't be determined otherwise. This requires pyorbital, skyfield, + and astropy to be installed and may need network access to obtain the TLE. + Note that even if ``use_tle`` is true, the TLE will not be used if + the dataset metadata contain the satellite position directly. Returns: - Geodetic longitude, latitude, altitude + Geodetic longitude, latitude, altitude [km] """ if preference is not None and preference not in ("nadir", "actual", "nominal", "projection"): @@ -323,6 +357,11 @@ def get_satpos( lon, lat = _get_sat_lonlat(data_arr, lonlat_prefixes) alt = _get_sat_altitude(data_arr, alt_prefixes) except KeyError: + if use_tle: + logger.warning( + "Orbital parameters missing from metadata. " + "Calculating from TLE using skyfield and astropy.") + return _get_satpos_from_platform_name(data_arr) raise KeyError("Unable to determine satellite position. Either the " "reader doesn't provide that information or " "geolocation datasets were not available.") @@ -363,6 +402,30 @@ def _get_sat_lonlat(data_arr, key_prefixes): return lon, lat +def _get_satpos_from_platform_name(cth_dataset): + """Get satellite position if no orbital parameters in metadata. + + Some cloud top height datasets lack orbital parameter information in + metadata. Here, orbital parameters are calculated based on the platform + name and start time, via Two Line Element (TLE) information. + + Needs pyorbital, skyfield, and astropy to be installed. + """ + from pyorbital.orbital import tlefile + from skyfield.api import EarthSatellite, load + from skyfield.toposlib import wgs84 + + name = cth_dataset.attrs["platform_name"] + tle = tlefile.read(name) + es = EarthSatellite(tle.line1, tle.line2, name) + ts = load.timescale() + gc = es.at(ts.from_datetime( + cth_dataset.attrs["start_time"].replace(tzinfo=datetime.timezone.utc))) + (lat, lon) = wgs84.latlon_of(gc) + height = wgs84.height_of(gc).to("km") + return (lon.degrees, lat.degrees, height.value) + + def _get_first_available_item(data_dict, possible_keys): for possible_key in possible_keys: try: diff --git a/setup.py b/setup.py index 236b66e346..ab8a55c315 100644 --- a/setup.py +++ b/setup.py @@ -30,7 +30,7 @@ except ImportError: pass -requires = ['numpy >=1.13', 'pillow', 'pyresample >=1.11.0', 'trollsift', +requires = ['numpy >=1.13', 'pillow', 'pyresample >=1.24.0', 'trollsift', 'trollimage >1.10.1', 'pykdtree', 'pyyaml', 'xarray >=0.10.1, !=0.13.0', 'dask[array] >=0.17.1', 'pyproj>=2.2', 'zarr', 'donfig', 'appdirs', 'pooch'] @@ -81,6 +81,7 @@ # Other 'geoviews': ['geoviews'], 'overlays': ['pycoast', 'pydecorate'], + 'satpos_from_tle': ['skyfield', 'astropy'], 'tests': test_requires, } all_extras = []