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()