matplotlib / matplotlib/basemap

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

Ouverte
#646 0 commentaires 0 réactions 0 personnes assignées Voir sur GitHub

Personne n'a encore pris cette issue.

Langage dominant
Python
Étoiles
817
Forks
395
Métriques de merge des PR
Aucune PR mergée en 30 j

Description

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:
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

Guide de contribution

Aucun guide de contribution indexé pour ce dépôt

Par où commencer

  1. Lisez l'issue en entier, puis le guide de contribution du projet.
  2. Signalez en commentaire que vous la prenez — cela évite que deux personnes fassent le même travail.
  3. Forkez le dépôt et travaillez sur une branche.
  4. Ouvrez une pull request qui référence le numéro de l'issue.

Piste de recherche

Reproduisez le MVE fourni avec la projection LCC de Basemap et ses arguments de coin, puis examinez comment la limite de la carte est établie par l'initialisation de Basemap et le comportement de drawmapboundary. Comparez le rectangle projeté avec les bords interpolés de latitude/longitude ; l'issue ne sera terminée que lorsqu'un comportement convenu pour la limite ou son expansion empêchera que des données valides de régions géographiques soient tronquées.

Rédigé par le modèle d'indexation à partir du texte de l'issue.

Évaluation

Stack technique
matplotlib, numpy, python
Domaine
data-visualization
Type d'issue
Bug
Difficulté
5/5
Temps estimé
Plus d'une semaine
Activité
À l'abandon
Clarté
Plutôt claire
Accessibilité débutants
35/100

Recevez les nouvelles issues par e-mail

Un résumé court des issues GitHub adaptées aux débutants.