matplotlib / matplotlib/basemap

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

Offen
#646 0 Kommentare 0 Reaktionen 0 zugewiesene Personen Auf GitHub ansehen

Dieses Issue hat noch niemand übernommen.

Vorherrschende Sprache
Python
Sterne
817
Forks
395
PR-Merge-Kennzahlen
Keine gemergten PRs in 30 T.

Beschreibung

### 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`

Beitragsleitfaden

Für dieses Repository ist kein Beitragsleitfaden indexiert

Erste Schritte

  1. Lies das ganze Issue und danach den Beitragsleitfaden des Projekts.
  2. Schreib ins Issue, dass du es übernimmst — das erspart doppelte Arbeit.
  3. Forke das Repository und arbeite in einem Branch.
  4. Öffne einen Pull Request, der die Issue-Nummer nennt.

Rechercherichtung

Reproduziere das bereitgestellte MVE mit Basemap's LCC-Projektion und seinen Eckargumenten und untersuche anschließend, wie die Kartengrenze durch die Initialisierung von Basemap und das Verhalten von drawmapboundary festgelegt wird. Vergleiche das projizierte Rechteck mit den interpolierten Breitengrad-/Längengradkanten; das Issue ist erst abgeschlossen, wenn ein vereinbartes Verhalten für die Grenze oder deren Erweiterung verhindert, dass gültige Daten geografischer Regionen abgeschnitten werden.

Vom Indexierungsmodell aus dem Issue-Text verfasst.

Bewertung

Tech-Stack
matplotlib, numpy, python
Bereich
data-visualization
Issue-Typ
Bug
Schwierigkeit
5/5
Geschätzter Aufwand
Über eine Woche
Aktivitätsstatus
Veraltet
Klarheit
Größtenteils klar
Anfängerfreundlichkeit
35/100

Neue Issues direkt in Ihr Postfach

Eine kurze Übersicht über anfängerfreundliche GitHub-Issues.