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

[antoni-upc]

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:

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