matplotlib / matplotlib/basemap

Discrepancy in LCC (and maybe others) between Lat/Lon corners and Projected Map Boundary

Abierto
#646 0 comentarios 0 reacciones 0 asignados Ver en GitHub

Nadie ha tomado este issue todavía.

Lenguaje dominante
Python
Estrellas
817
Forks
395
Métricas de merge de PR
Sin PR fusionados en 30 d

Descripción

### Description:

The reason for this issue is to point out the discrepancy in the behaviour of Basemap when using the Lambert Conformal Conic (LCC) projection with the corner arguments (`llcrnrlon`, `llcrnrlat`, `urcrnrlon`, `urcrnrlat`).

### Issue:

When using an LCC map, I expect the resulting figure to respect the geographic boundaries defined by the Lower-Left (BL) and Upper-Right (TR) corners, e.g. _lines of constant Latitude and of constant Longitude_. However, it appears that Basemap uses the projected X/Y coordinates of these two points and then draws a rectangle in the projected plane connecting them.

Because LCC parallels are curved and meridians converge, a rectangle drawn in the projected X/Y plane does **not** correspond to the geographic box defined by the input coordinates.

The consequence is that data plotted near the boundaries may not be displayed. Parts of the region that are technically within the Latitude/Longitude limits provided are "cut off" because they fall outside the calculated rectangular X/Y box.

### Expected Behaviour:

When a user provides specific lower-left and upper-right coordinates, the expectation is that the resulting map respects those geographic boundaries. I propose two potential solutions:

- **Non-Rectangular Clipping (Cone Section)**: The map boundary (axis spines) should follow the actual meridians and parallels connecting the corners, rather than being forced into a rectangle. This would result in a "cone section" or wedge-shaped plot that accurately reflects the geometry of the projection and the user's requested region.

- **Expand the Projection Box**: Alternatively, if the plot must remain rectangular, the calculated clipping box in the projected plane should be automatically expanded to encompass the maximum extent of the curved geographic region, ensuring that no data strictly within the provided Lat/Lon limits is cut off.

### Example of the issue:

Image

From the image above:

- Red Box: The map boundary created by Basemap (Projected Rectangle).
- Blue Box: The geographic region requested (following Lat/Lon lines).

### MVE:

```python
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np

# --- 1. Map Initialisation ---
fig = plt.figure(figsize=(12, 10))

m = Basemap(
projection='lcc',
resolution='l', # Low resolution
llcrnrlon=-75.0, # Map Lower-Left Longitude
llcrnrlat=25.0, # Map Lower-Left Latitude
urcrnrlon=30.0, # Map Upper-Right Longitude
urcrnrlat=55.0, # Map Upper-Right Latitude
lat_0=55, # Central Latitude (approx mid-point)
lon_0=-30, # Central Longitude (approx mid-point)
lat_1=30, lat_2=60 # Standard Parallels
)

# --- 2. Visual Styling ---
# Pastel colors for clear contrast against the red/blue boxes
water_color = '#AEE1FA'
land_color = '#B2D8B2'

m.drawmapboundary(fill_color=water_color)
m.fillcontinents(color=land_color, lake_color=water_color)
m.drawcoastlines(linewidth=0.6)
m.drawcountries(linewidth=0.5, color='gray')

# Draw grid to visualize curvature
m.drawparallels(np.arange(0, 90, 10), labels=[1, 0, 0, 0], color='gray', linewidth=0.5)
m.drawmeridians(np.arange(-100, 40, 10), labels=[0, 0, 0, 1], color='gray', linewidth=0.5)

# Shared Corner Coordinates
# BL: 40N, 70W | TR: 60N, 0E
bl_lat, bl_lon = 40.0, -70.0
tr_lat, tr_lon = 60.0, 0.0

# Convert Geographic corners to Projected (Meters) corners
bl_x, bl_y = m(bl_lon, bl_lat)
tr_x, tr_y = m(tr_lon, tr_lat)

# --- A. RED BOX: Projected Rectangle (90-degree corners) ---
# This box is defined in the X/Y projection plane.
# It draws straight lines between the projected corners.
# Result: A perfect rectangle on screen, but does not follow lat/lon lines.
proj_x = [bl_x, tr_x, tr_x, bl_x, bl_x]
proj_y = [bl_y, bl_y, tr_y, tr_y, bl_y]

m.plot(proj_x, proj_y, color='red', linewidth=3, label='Projected Box (X/Y Rectangle)')

# --- B. BLUE BOX: Geographic Box (Follows Lat/Lon) ---
# This box is defined by constant Latitude and Longitude.
# Because LCC is conic, these lines curve in the projection.
# We interpolate points to show the true geographic boundary.
lons = np.linspace(bl_lon, tr_lon, 50)
lats = np.linspace(bl_lat, tr_lat, 50)

# Top & Bottom Edges (Constant Latitude)
m.plot(lons, [bl_lat] * 50, latlon=True, color='blue', linewidth=3)
m.plot(lons, [tr_lat] * 50, latlon=True, color='blue', linewidth=3)

# Left & Right Edges (Constant Longitude)
m.plot([bl_lon] * 50, lats, latlon=True, color='blue', linewidth=3)
m.plot([tr_lon] * 50, lats, latlon=True, color='blue', linewidth=3, label='Geographic Box (Lat/Lon)')

# 1. Plot Shared Corners
m.plot(bl_x, bl_y, 'ko', markersize=8, zorder=10)
plt.text(bl_x, bl_y - 150000, 'BL\n(40N, 70W)', ha='center', va='top', fontweight='bold')

m.plot(tr_x, tr_y, 'ko', markersize=8, zorder=10)
plt.text(tr_x, tr_y + 100000, 'TR\n(60N, 0E)', ha='center', va='bottom', fontweight='bold')

# 2. Projected Centre (Red Square)
# The mathematical midpoint in the X/Y plane
mid_x = (bl_x + tr_x) / 2.0
mid_y = (bl_y + tr_y) / 2.0
m.plot(mid_x, mid_y, 'rs', markersize=8, zorder=10)
plt.text(mid_x, mid_y - 200000, 'Proj\nCenter', color='red', ha='center', va='top', fontsize=9, fontweight='bold')

# 3. Geographic Centre (Blue Triangle)
# The mathematical midpoint in Degrees, then projected
mid_lon = (bl_lon + tr_lon) / 2.0
mid_lat = (bl_lat + tr_lat) / 2.0
geo_mid_x, geo_mid_y = m(mid_lon, mid_lat)
m.plot(geo_mid_x, geo_mid_y, 'b^', markersize=8, zorder=10)
plt.text(geo_mid_x, geo_mid_y + 150000, 'Geo\nCenter', color='blue', ha='center', va='bottom', fontsize=9,
fontweight='bold')

plt.legend(loc='lower right')
plt.show()
```

Environment:

- Python version: `3.12.8`
- Matplotlib version: `3.10.0`
- Basemap version: `1.4.1`

Guía de contribución

No hay ninguna guía de contribución indexada para este repositorio

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.

Línea de trabajo

Reproduce el MVE proporcionado con la proyección LCC de Basemap y sus argumentos de esquina; después, inspecciona cómo se establece el límite del mapa mediante la inicialización de Basemap y el comportamiento de drawmapboundary. Compara el rectángulo proyectado con los bordes interpolados de latitud/longitud; el issue solo estará terminado cuando un comportamiento acordado para el límite o su expansión evite que se recorten datos válidos de regiones geográficas.

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

Evaluación

Stack tecnológico
matplotlib, numpy, python
Área
data-visualization
Tipo de issue
Error
Dificultad
5/5
Tiempo estimado
Más de una semana
Estado de actividad
Estancado
Claridad
Bastante claro
Aptitud para principiantes
35/100

Recibe los nuevos issues en tu correo

Un resumen breve de issues de GitHub para principiantes.