Wrong coordinates returned by `AreaDefintion.get_lonlats` for some projections in out-of-Earth locations

Abierto
  • ¿Algo no está bien?
#558 6 comentarios 0 reacciones 0 asignados Ver en GitHub

Los mantenedores suelen responder en 1 día

Nadie ha tomado este issue todavía.

Evaluación

Dificultad
4/5
Tiempo estimado
3-5 días
Aptitud para principiantes
35/100
Tipo de issue
Error
Claridad
Bastante claro
Estado de actividad
Estancado
Stack tecnológico
numpy, python
Área
data

Línea de trabajo

Comienza ejecutando el ejemplo de proyección de cartopy proporcionado e inspecciona AreaDefinition.get_lonlats(), get_proj_coords() y la llamada inversa a pyproj Transformer mostrada en el issue. Compara los resultados fuera de la Tierra con los valores inf esperados para las proyecciones indicadas; después, añade cobertura de regresión que demuestre el manejo correcto de coordenadas no válidas y verifica que el remuestreo ya no utilice los valores erróneos.

Escrito por el modelo de indexación a partir del texto del issue.

Descripción

In the context of https://github.com/pytroll/pyresample/pull/546, I was testing that the new boundary extraction methodology works for all projections available in cartopy, and I stumbled in wrong values returned by AreaDefinition.get_lonlats().
After a quick investigation, the error seems to arise from wrong results by pyproj transformations !

Here below I provide a code that enable to visualize the wrong coordinates as well as snapshot of what is going on for some of the problematic projections.

Problem description

The wrong results occurs only for some CRS in out-of-Earth locations.
Usually pyproj returns np.inf values in out-of-Earth locations, while with the projections listed below, AreaDefinition.get_lonlats() returns wrong/bad longitude/latitude values. In a couple of cases, even returns latitudes outside the (-90, 90) bounds.

NOTE1: These errors causes wrong resampling with pyresample !

NOTE2: Note that for all cartopy projections not listed below, pyproj returns the correct results !

Expected Output

Usually pyproj returns np.inf values in out-of-Earth locations.

Code Sample, a minimal, complete, and verifiable piece of code
import numpy as np 
import pyproj 
import cartopy
from pyresample.future.geometry import AreaDefinition 
import matplotlib.pyplot as plt 


list_bugs = [ 
    cartopy.crs.EquidistantConic(central_longitude=0.0, central_latitude=0.0, false_easting=0.0, false_northing=0.0, standard_parallels=(20.0, 50.0), globe=None),
    cartopy.crs.AzimuthalEquidistant(central_longitude=0.0, central_latitude=0.0, false_easting=0.0, false_northing=0.0, globe=None),
    cartopy.crs.AlbersEqualArea(central_longitude=0.0, central_latitude=0.0, false_easting=0.0, false_northing=0.0, standard_parallels=(20.0, 50.0), globe=None),
    cartopy.crs.LambertConformal(central_longitude=-96.0, central_latitude=39.0, false_easting=0.0, false_northing=0.0, standard_parallels=(33, 45), globe=None, cutoff=-30),
    cartopy.crs.EqualEarth(central_longitude=0, false_easting=None, false_northing=None, globe=None),
    cartopy.crs.Hammer(central_longitude=0, false_easting=None, false_northing=None, globe=None),
    cartopy.crs.Sinusoidal(central_longitude=0.0, false_easting=0.0, false_northing=0.0, globe=None),
    cartopy.crs.EckertI(central_longitude=0, false_easting=None, false_northing=None, globe=None),
    cartopy.crs.EckertII(central_longitude=0, false_easting=None, false_northing=None, globe=None),
    cartopy.crs.EckertIII(central_longitude=0, false_easting=None, false_northing=None, globe=None),
    cartopy.crs.EckertIV(central_longitude=0, false_easting=None, false_northing=None, globe=None),
    cartopy.crs.EckertV(central_longitude=0, false_easting=None, false_northing=None, globe=None),
    cartopy.crs.EckertVI(central_longitude=0, false_easting=None, false_northing=None, globe=None),
]

cartopy_crs = list_bugs[1]
for cartopy_crs in list_bugs:
    name = type(cartopy_crs).__name__
    x_min, x_max = cartopy_crs.x_limits
    y_min, y_max = cartopy_crs.y_limits
    area_extent = [x_min, y_min, x_max, y_max]
    proj_dict = cartopy_crs.to_dict()
    shape = (200, 300)
    area_def = AreaDefinition(
        crs=proj_dict,
        shape=shape,
        area_extent=area_extent
    )
    lons, lats = area_def.get_lonlats()
    
    plt.imshow(lons)
    plt.title(name+"Longitude"); plt.colorbar()
    plt.show()
    
    plt.imshow(lats)
    plt.title(name+" Latitude"); plt.colorbar()
    plt.show()

To further investigate the erroneous results, just take values in the lower left corner (which are out of Earth in the specified projections) and look at the returned values:

from pyresample.geometry import get_geodetic_crs_with_no_datum_shift, TransformDirection, CRS
self = area_def
proj_wkt = self.crs_wkt

x, y = area_def.get_proj_coords(data_slice=(-1, 0))
crs = CRS.from_wkt(proj_wkt)
gcrs = get_geodetic_crs_with_no_datum_shift(crs)
transformer = pyproj.Transformer.from_crs(gcrs, crs, always_xy=True)
lon, lat = transformer.transform(x, y, direction=TransformDirection.INVERSE)
print(lon, lat)   # should be inf, inf
Examples

Equidistant Conic
image
image

Azimuthal Equidistant
image
image

Albers Equal Area
image
image

Equal Earth
image
image

Lamber Conformal
image
image

Sinusoidal
image
image

Eckert
image
image

Lenguaje dominante
Python
Estrellas
388
Forks
104
Merge medio
2 d 16 h
PR fusionados (30 d)
3

Preparar el entorno

  • Sin Dockerfile ni archivo de Docker Compose
  • Tiene una plantilla de pull request
  • Sin guía de contribución

Primeros pasos

  1. Lee el issue completo y luego la guía de contribución del proyecto.
  2. Comenta en el issue que vas a ocuparte — evita que dos personas hagan lo mismo.
  3. Haz un fork del repositorio y trabaja en una rama.
  4. Abre un pull request que haga referencia al número del issue.

Más de pytroll/pyresample

Todos los issues de pytroll/pyresample

Issues similares

Más issues de Python

Recibe los nuevos issues en tu correo

Un resumen breve de issues de GitHub para principiantes.