Source code for satpy.modifiers.parallax

# 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__)

[docs] class MissingHeightError(ValueError): """Raised when heights do not overlap with area to be corrected."""
[docs] class IncompleteHeightWarning(UserWarning): """Raised when heights only partially overlap with area to be corrected."""
[docs] 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])
[docs] 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
[docs] 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)
[docs] 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
[docs] 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))
[docs] 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): # noqa: D417 """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)
[docs] 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',!s}, " f"with base area {!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, 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)
[docs] @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 # and return SwathDefinition( xr.DataArray(lon, dims=("y", "x")), xr.DataArray(lat, dims=("y", "x")))
[docs] 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
[docs] 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 implemented. Waiting for " "fix for", stacklevel=3 )
[docs] 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)
[docs] 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
[docs] 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
[docs] 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)