[Git][debian-gis-team/pyorbital][master] 6 commits: New upstream version 1.10.2
Antonio Valentino (@antonio.valentino)
gitlab at salsa.debian.org
Sun Aug 10 17:51:16 BST 2025
Antonio Valentino pushed to branch master at Debian GIS Project / pyorbital
Commits:
78f1e478 by Antonio Valentino at 2025-08-09T09:10:22+00:00
New upstream version 1.10.2
- - - - -
c5132147 by Antonio Valentino at 2025-08-09T09:10:22+00:00
Update upstream source from tag 'upstream/1.10.2'
Update to upstream version '1.10.2'
with Debian dir b50d6e8155d912c0f99b3798baf39d6d3ed3fe38
- - - - -
d756127b by Antonio Valentino at 2025-08-09T09:10:51+00:00
New upstream release
- - - - -
2812efc1 by Antonio Valentino at 2025-08-09T09:20:40+00:00
Update dependencies
- - - - -
4f8be6dc by Antonio Valentino at 2025-08-10T16:48:27+00:00
Set upstream metadata fields: Documentation.
Changes-By: lintian-brush
- - - - -
a94a6960 by Antonio Valentino at 2025-08-10T16:48:29+00:00
Set distribution to unstable
- - - - -
9 changed files:
- .pre-commit-config.yaml
- CHANGELOG.md
- debian/changelog
- debian/control
- debian/upstream/metadata
- pyorbital/geoloc.py
- + pyorbital/geoloc_avhrr.py
- pyorbital/tests/test_geoloc.py
- pyproject.toml
Changes:
=====================================
.pre-commit-config.yaml
=====================================
@@ -3,7 +3,7 @@ fail_fast: false
repos:
- repo: https://github.com/astral-sh/ruff-pre-commit
# Ruff version.
- rev: 'v0.11.4'
+ rev: 'v0.11.12'
hooks:
- id: ruff
- repo: https://github.com/pre-commit/pre-commit-hooks
@@ -20,7 +20,7 @@ repos:
- id: bandit
args: [--ini, .bandit]
- repo: https://github.com/pre-commit/mirrors-mypy
- rev: 'v1.15.0' # Use the sha / tag you want to point at
+ rev: 'v1.16.0' # Use the sha / tag you want to point at
hooks:
- id: mypy
additional_dependencies:
=====================================
CHANGELOG.md
=====================================
@@ -1,3 +1,27 @@
+## Version 1.10.2 (2025/06/09)
+
+
+### Pull Requests Merged
+
+#### Features added
+
+* [PR 192](https://github.com/pytroll/pyorbital/pull/192) - Return distances from GCP minimization
+
+In this release 1 pull request was closed.
+
+
+## Version 1.10.1 (2025/05/08)
+
+
+### Pull Requests Merged
+
+#### Features added
+
+* [PR 189](https://github.com/pytroll/pyorbital/pull/189) - Work with gcps for avhrr
+
+In this release 1 pull request was closed.
+
+
## Version 1.10.0 (2025/04/10)
### Issues Closed
=====================================
debian/changelog
=====================================
@@ -1,3 +1,12 @@
+pyorbital (1.10.2-1) unstable; urgency=medium
+
+ * New upstream release.
+ * debian/control:
+ - Add dependency on pyproj.
+ * Set upstream metadata fields: Documentation.
+
+ -- Antonio Valentino <antonio.valentino at tiscali.it> Sat, 09 Aug 2025 09:16:30 +0000
+
pyorbital (1.10.0-1) unstable; urgency=medium
[ Bas Couwenberg ]
=====================================
debian/control
=====================================
@@ -17,6 +17,8 @@ Build-Depends: debhelper-compat (= 13),
python3-hatchling,
python3-hatch-vcs,
python3-numpy,
+ python3-numpy-dev,
+ python3-pyproj,
python3-pytest <!nocheck>,
python3-requests,
python3-scipy,
=====================================
debian/upstream/metadata
=====================================
@@ -1,5 +1,6 @@
---
Bug-Database: https://github.com/pytroll/pyorbital/issues
Bug-Submit: https://github.com/pytroll/pyorbital/issues/new
+Documentation: https://pyorbital.readthedocs.io/en/latest/
Repository: https://github.com/pytroll/pyorbital.git
Repository-Browse: https://github.com/pytroll/pyorbital
=====================================
pyorbital/geoloc.py
=====================================
@@ -25,18 +25,18 @@
# TODO:
# - Attitude correction
-# - project on an ellipsoid instead of a sphere
# - optimize !!!
# - test !!!
from __future__ import print_function
import numpy as np
+from pyproj import Transformer
# DIRTY STUFF. Needed the get_lonlatalt function to work on pos directly if
# we want to print out lonlats in the end.
from pyorbital import astronomy
-from pyorbital.orbital import XKMPER, F, Orbital
+from pyorbital.orbital import Orbital
A = 6378.137 # WGS84 Equatorial radius (km)
B = 6356.75231414 # km, GRS80
@@ -183,29 +183,12 @@ def qrotate(vector, axis, angle):
def get_lonlatalt(pos, utc_time):
- """Calculate sublon, sublat and altitude of satellite, considering the earth an ellipsoid.
-
- http://celestrak.com/columns/v02n03/
-
- """
- (pos_x, pos_y, pos_z) = pos / XKMPER
- lon = ((np.arctan2(pos_y * XKMPER, pos_x * XKMPER) - astronomy.gmst(utc_time)) % (2 * np.pi))
- lon = np.where(lon > np.pi, lon - np.pi * 2, lon)
- lon = np.where(lon <= -np.pi, lon + np.pi * 2, lon)
-
- r = np.sqrt(pos_x ** 2 + pos_y ** 2)
- lat = np.arctan2(pos_z, r)
- e2 = F * (2 - F)
-
- while True:
- lat2 = lat
- c = 1 / (np.sqrt(1 - e2 * (np.sin(lat2) ** 2)))
- lat = np.arctan2(pos_z + c * e2 * np.sin(lat2), r)
- if np.all(abs(lat - lat2) < 1e-10):
- break
- alt = r / np.cos(lat) - c
- alt *= A
- return np.rad2deg(lon), np.rad2deg(lat), alt
+ """Calculate sublon, sublat and altitude of satellite, considering the earth an ellipsoid."""
+ trf = Transformer.from_crs(dict(proj="geocent"), dict(proj="latlong"))
+ lon, lat, alt = trf.transform(*(pos * 1000))
+ lon = lon - np.rad2deg(astronomy.gmst(utc_time))
+ lon = np.where(lon < -180, lon + 360, lon)
+ return lon, lat, alt
# END OF DIRTY STUFF
=====================================
pyorbital/geoloc_avhrr.py
=====================================
@@ -0,0 +1,157 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+# Copyright (c) 2025 PyTroll Community
+
+# Author(s):
+
+# Martin Raspaud <martin.raspaud at smhi.se>
+
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""This module provides geoloc operations specific to the avhrr instrument.
+
+In particular, it provides functions that allow matching gcp location in swath coordinates to reference positions, and
+then minimise the distance to these positions by adjusting the time offset and attitude.
+"""
+
+import logging
+
+import numpy as np
+from pyproj import Geod
+
+from pyorbital.geoloc import ScanGeometry, compute_pixels, get_lonlatalt
+
+logger = logging.getLogger(__name__)
+geod = Geod(ellps="WGS84")
+
+def compute_avhrr_gcps_lonlatalt(gcps, max_scan_angle, rpy, start_time, tle) -> None:
+ """Compute the longitute, latitude and altitude of given gcps (scanlines, columns of the swath).
+
+ The gcps are arbitrary location in swath coordinates, for example (10.3, 7.7) for a gcp at line 10.3 in the swath,
+ and column 7.7. This function returns the geographical coordinates of the gcps.
+
+ The scanlines are relative to the pass scanline numbers, zero-based.
+ """
+ time_line_interval = 1/6
+ time_row_interval = 25e-6
+
+ fov_x = gcps[:, 1]
+ fov_y = gcps[:, 0]
+
+ scan_angles_across = (fov_x / 1023.5 - 1) * np.deg2rad(-max_scan_angle)
+ scan_angles_along = np.zeros_like(scan_angles_across)
+ scan_angles = np.vstack((scan_angles_across, scan_angles_along))
+ time_offsets = np.array(fov_x * time_row_interval + fov_y * time_line_interval)
+ geom = ScanGeometry(scan_angles, time_offsets)
+ start_time = np.datetime64(start_time)
+ s_times = geom.times(start_time)
+
+ pixels_pos = compute_pixels(tle, geom, s_times, rpy)
+ return get_lonlatalt(pixels_pos, s_times)
+
+
+def estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, start_time, tle, max_scan_angle):
+ """Estimate time offset and attitude deviations from gcps.
+
+ Provided reference longitudes and latitudes for the gcps, this function minimises the attitude and time offset
+ needed to match the gcp coordinates to the reference coordinates.
+ """
+ from scipy.optimize import minimize
+
+ original_distances = compute_gcp_distances_to_reference_lonlats((0, 0, 0, 0), gcps, start_time, tle, max_scan_angle,
+ (ref_lons, ref_lats))
+ original_median_distance = np.median(original_distances)
+ logger.debug(f"GCP distances: median {original_median_distance}, std {np.std(original_distances)}")
+ # we need to work in seconds*1e3 to avoid the nanosecond precision issue
+ res = minimize(compute_gcp_accumulated_squared_distances_to_reference_lonlats,
+ x0=(0, 0, 0, 0),
+ args=(gcps, start_time, tle, max_scan_angle, (ref_lons, ref_lats)),
+ bounds=((-0.007, 0.007) , (-0.5, 0.5), (-0.5, 0.5), (-0.5, 0.5)))
+ if not res.success:
+ raise RuntimeError("Time and attitude estimation did not converge")
+ time_diff, roll, pitch, yaw = res.x * [1e3, 1, 1, 1]
+ logger.debug(f"Estimated time difference to {time_diff} seconds, "
+ f"attitude to {np.rad2deg(roll)}, {np.rad2deg(pitch)}, {np.rad2deg(yaw)} degrees")
+ distances = compute_gcp_distances_to_reference_lonlats(res.x, gcps, start_time, tle, max_scan_angle,
+ (ref_lons, ref_lats))
+
+ minimized_median_distance = np.median(distances)
+ logger.debug(f"Remaining GCP distances: median {minimized_median_distance}, std {np.std(distances)}")
+
+ return time_diff, (roll, pitch, yaw), (original_distances, distances)
+
+
+def estimate_time_offset(gcps, ref_lons, ref_lats, start_time, tle, max_scan_angle):
+ """Estimate time offset from gcps.
+
+ Provided reference longitudes and latitudes for the gcps, this function minimises the time offset
+ needed to match the gcp coordinates to the reference coordinates.
+ """
+ from scipy.optimize import minimize
+
+ original_distances = compute_gcp_distances_to_reference_lonlats((0, 0, 0, 0), gcps, start_time, tle, max_scan_angle,
+ (ref_lons, ref_lats))
+ original_median_distance = np.median(original_distances)
+ logger.debug(f"GCP distances: median {original_median_distance}, std {np.std(original_distances)}")
+
+ def gcp_distance_for_time(time):
+ dist = compute_gcp_accumulated_squared_distances_to_reference_lonlats((time[0], 0, 0, 0), gcps, start_time, tle,
+ max_scan_angle, (ref_lons, ref_lats))
+ return dist
+
+ # we need to work in seconds*1e3 to avoid the nanosecond precision issue
+ res = minimize(gcp_distance_for_time,
+ x0=(0,),
+ bounds=((-0.03, 0.03),),
+ options=dict(ftol=1e-1),
+ )
+ if not res.success:
+ raise RuntimeError("Time offset estimation did not converge")
+ time_diff, = res.x * [1e3,]
+ logger.debug(f"Estimated time difference to {time_diff} seconds")
+ distances = compute_gcp_distances_to_reference_lonlats((res.x[0], 0, 0, 0), gcps, start_time, tle, max_scan_angle,
+ (ref_lons, ref_lats))
+
+ minimized_median_distance = np.median(distances)
+ logger.debug(f"Remaining GCP distances: median {minimized_median_distance}, std {np.std(distances)}")
+
+ return time_diff, (original_distances, distances)
+
+
+def compute_gcp_accumulated_squared_distances_to_reference_lonlats(
+ variables, gcps, start_time, tle, max_scan_angle, refs):
+ """Compute the summed squared distance fot gcps to reference lonlats.
+
+ Given the gcps (in swath coordinates) along with attitude and time offset, compute the sum of squared distances to
+ the reference lons and lats of the gcps.
+ """
+ distances = compute_gcp_distances_to_reference_lonlats(variables, gcps, start_time, tle, max_scan_angle, refs)
+ return np.sum(distances**2)
+
+
+def compute_gcp_distances_to_reference_lonlats(variables, gcps, start_time, tle, max_scan_angle, refs):
+ """Compute the gcp distances to references lonlats."""
+ time_diff, roll, pitch, yaw = variables
+ # we need to work in seconds*1e3 to avoid the nanosecond precision issue
+ time = np.datetime64(start_time) + np.timedelta64(int(time_diff * 1e12), "ns")
+ lons, lats, _ = compute_avhrr_gcps_lonlatalt(gcps, max_scan_angle, (roll, pitch, yaw), time, tle)
+ valid = np.isfinite(lons)
+ lons = lons[valid]
+ lats = lats[valid]
+ ref_lons, ref_lats = refs
+ ref_lons = np.array(ref_lons)[valid]
+ ref_lats = np.array(ref_lats)[valid]
+ _, _, distances = geod.inv(ref_lons, ref_lats, lons, lats)
+ return distances
=====================================
pyorbital/tests/test_geoloc.py
=====================================
@@ -26,8 +26,14 @@
import datetime as dt
import numpy as np
+import pytest
from pyorbital.geoloc import ScanGeometry, geodetic_lat, qrotate, subpoint
+from pyorbital.geoloc_avhrr import (
+ compute_avhrr_gcps_lonlatalt,
+ estimate_time_and_attitude_deviations,
+ estimate_time_offset,
+)
from pyorbital.geoloc_instrument_definitions import (
amsua,
ascat,
@@ -173,6 +179,81 @@ class TestGeoloc:
4497.06396339]), rtol=1e-8, atol=1e-8)
+
+
+def test_arbitrary_point_geoloc():
+ """Test geolocating an arbitrary point in the swath."""
+ from pyorbital.geoloc_avhrr import compute_avhrr_gcps_lonlatalt
+
+ # Couple of example Two Line Elements
+ tle1 = "1 33591U 09005A 12345.45213434 .00000391 00000-0 24004-3 0 6113"
+ tle2 = "2 33591 098.8821 283.2036 0013384 242.4835 117.4960 14.11432063197875"
+
+ # Choosing a specific time, this should be relatively close to the issue date of the TLE
+ t = dt.datetime(2012, 12, 12, 4, 16, 1, 575000)
+ rpy = (0, 0, 0)
+
+ max_scan_angle = 55.37
+
+ gcps = np.array([[2, 500], [1500, 700], [20, 1000]])
+
+ lons, lats, alts = compute_avhrr_gcps_lonlatalt(gcps, max_scan_angle, rpy, t, (tle1, tle2))
+
+ assert lons[0] == pytest.approx(-34.69996894)
+ assert lats[0] == pytest.approx(56.69799502)
+
+ assert lons[2] == pytest.approx(-27.573052737698944)
+ assert lats[2] == pytest.approx(55.626740897592654)
+
+
+def test_minimize_geoloc_error():
+ """Test minimizing the distance to a set of gcps."""
+ # Couple of example Two Line Elements
+ tle1 = "1 33591U 09005A 12345.45213434 .00000391 00000-0 24004-3 0 6113"
+ tle2 = "2 33591 098.8821 283.2036 0013384 242.4835 117.4960 14.11432063197875"
+ tle = (tle1, tle2)
+
+ # Choosing a specific time, this should be relatively close to the issue date of the TLE
+ t = dt.datetime(2012, 12, 12, 4, 16, 1, 575000)
+
+ ref_time_displacement = 0.51
+ ref_time = t + dt.timedelta(seconds=ref_time_displacement)
+ ref_yaw = 0.1
+ rpy = (0, 0, ref_yaw)
+ max_scan_angle = 55.37
+ # gcps are line/col
+ gcps = np.array([[2, 500], [1500, 700], [20, 1000], [500, 1100], [100, 2000]])
+ ref_lons, ref_lats, _ = compute_avhrr_gcps_lonlatalt(gcps, max_scan_angle, rpy, ref_time, tle)
+ time_diff, (roll, pitch, yaw), (do, dm) = estimate_time_and_attitude_deviations(gcps, ref_lons, ref_lats, t,
+ tle, max_scan_angle)
+ assert time_diff == pytest.approx(ref_time_displacement, abs=1e-2)
+ assert yaw == pytest.approx(ref_yaw, abs=1e-2)
+ assert min(do) > max(dm)
+
+
+def test_minimize_time_error():
+ """Test minimizing the distance to a set of gcps using only time offset."""
+ # Couple of example Two Line Elements
+ tle1 = "1 33591U 09005A 12345.45213434 .00000391 00000-0 24004-3 0 6113"
+ tle2 = "2 33591 098.8821 283.2036 0013384 242.4835 117.4960 14.11432063197875"
+ tle = (tle1, tle2)
+
+ # Choosing a specific time, this should be relatively close to the issue date of the TLE
+ t = dt.datetime(2012, 12, 12, 4, 16, 1, 575000)
+
+ ref_time_displacement = 20
+ ref_time = t + dt.timedelta(seconds=ref_time_displacement)
+ rpy = (0, 0, 0)
+ max_scan_angle = 55.37
+ # gcps are line/col
+ gcps = np.array([[2, 500], [1500, 700], [20, 1000], [500, 1100], [100, 2000]])
+ ref_lons, ref_lats, _ = compute_avhrr_gcps_lonlatalt(gcps, max_scan_angle, rpy, ref_time, tle)
+ time_diff, (do, dm) = estimate_time_offset(gcps, ref_lons, ref_lats, t,
+ tle, max_scan_angle)
+ assert time_diff == pytest.approx(ref_time_displacement, abs=1e-2)
+ assert min(do) > max(dm)
+
+
class TestGeolocDefs:
"""Test the instrument definitions."""
=====================================
pyproject.toml
=====================================
@@ -9,6 +9,7 @@ dependencies = ["numpy>=1.19.0",
"scipy",
"requests",
"defusedxml",
+ "pyproj>=3.7.1",
]
readme = "README.md"
requires-python = ">=3.10"
@@ -32,6 +33,11 @@ fetch_tles = "pyorbital.fetch_tles:run"
[project.optional-dependencies]
doc = ["sphinx", "sphinx_rtd_theme", "sphinxcontrib-apidoc"]
+test = [
+ "dask>=2025.3.0",
+ "pytest>=8.3.5",
+ "xarray>=2025.3.1",
+]
[build-system]
requires = ["hatchling", "hatch-vcs"]
View it on GitLab: https://salsa.debian.org/debian-gis-team/pyorbital/-/compare/a756c6320eb84285503640e56d02054e8f11c75c...a94a69604d411e15673c3ac74ce4f132b2aa850f
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/pyorbital/-/compare/a756c6320eb84285503640e56d02054e8f11c75c...a94a69604d411e15673c3ac74ce4f132b2aa850f
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/pkg-grass-devel/attachments/20250810/7a4a2b60/attachment-0001.htm>
More information about the Pkg-grass-devel
mailing list