# showMap.py
# Andrew Davison, ad@coe.psu.ac.th, July 2025
'''
  Draw GreatCircle and Rhumb line between two cities
  on an Orthographic, PlateCarree or Mercator map

  See 
    https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html

  Possible cities: = New York  Plymouth  Tokyo Berkeley
         Cape Town  Panama City  Seattle   Melbourne
         Reykjavik  Helsinki

Plot transform Notes
--------------------
If the plot involves lines, transform=crs.PlateCarree() transforms/draws straight lines along supplied data points

transform=crs.Geodetic() will add extra points along the supplied vertices to obtain the plots of long lines as greatcircle arcs that make the plots more accurate. – 

transform=crs.Geodetic() will do much more computations, so if you don't need greatcircle arcs in your plots, you should use transform=crs.PlateCarree(). 

Both provide similar coordinates transformation. 
'''
import cartopy.crs as ccrs
import cartopy.feature as cf
import matplotlib.pyplot as plt
import nav


def gnomonic(lat1, lon1, lat2, lon2):
  if lon1*lon2 < 0:  # a positive and negative coord
    dist1 = -lon1 if lon1 < 0 else lon1  # from prime meridian
    dist2 = -lon2 if lon2 < 0 else lon2
    cLon = (dist1 + dist2) / 2
  else:
    cLon = (lon1 + lon2) / 2
  cLat = (lat1 + lat2) / 2
  return ccrs.Gnomonic(central_longitude=cLon, 
                       central_latitude=cLat)



def plotCity(ax, city, lat, lon):
  # Add marker
  ax.plot(lon, lat, 'o', color='green', 
         transform=ccrs.PlateCarree())
  # label city
  ax.text(lon, lat, " "+city, color='green', 
         transform=ccrs.PlateCarree())


# get coordinates for the cities
startCity = "Perth"  
endCity =  "Nord"    # "Seattle"  # 
cities = nav.loadCityData()
lat1, lon1 = nav.readCoord(cities, startCity) 
lat2, lon2 = nav.readCoord(cities, endCity)
print(f"{startCity} (lat,lon): ({lat1:.2f}, {lon1:.2f})")
print(f"{endCity} (lat,lon): ({lat2:.2f}, {lon2:.2f})")
usesDL = nav.usesDateline(lon1, lon2)
if usesDL:
  print("Crosses the dateline")
print()

mapStr = input("Orthographic, PlateCarree, Gnomonic, or Mercator?  (o, p, g, or m) ")
if mapStr:
  mapType = mapStr[0].upper()
  if mapType not in ['O','P','G','M']:
    mapType = 'O'
else:
  mapType = 'O'
print("Map type:", mapType)


plt.figure(f"{startCity} to {endCity}", figsize=(8,6))

usesDL = nav.usesDateline(lon1, lon2)
if usesDL:
  print("Crosses the dateline")

if mapType == "M":
  proj = ccrs.Mercator()
elif mapType == "P":
  proj = ccrs.PlateCarree()
elif mapType == "G":
  proj = gnomonic(lat1, lon1, lat2, lon2)
else: 
  if usesDL:
    proj = ccrs.Orthographic(180, 30)  # on dateline
  else:
    proj = ccrs.Orthographic(70, 30)

proj.threshold = proj.threshold/20  # finer threshold
ax = plt.axes(projection=proj)

ax.set_global()  # show the entire globe
lineSep = 10 if mapType == "O" else 30
ax.gridlines(draw_labels=True, ls='--',
             xlocs=range(-180, 181, lineSep),
             ylocs=range(-90, 91, lineSep))
ax.coastlines()
ax.add_feature(cf.BORDERS, ls=':')


# use Cartopy to draw the Great Circle
plt.plot([lon1, lon2], [lat1, lat2],
     color='blue', lw=5, alpha=0.2,
     label="Great Circle", transform=ccrs.Geodetic())  
                 # spherical; great circle


# draw the Great Circle trip as a series of coordinates
GCStep = 20
if mapType == "M":  # motre points needed
  GCStep = 40
coords = nav.greatCircleCoords(lat1, lon1, lat2, lon2, GCStep)
print("\nGreat Circle (lat,lon)s and deg Bearing:")
for coord in coords:
  cLat, cLon, cBear = coord
  print(f"  ({cLat:.2f}, {cLon:.2f}); {cBear:.2f}")

lats, lons, _ = zip(*coords)
if usesDL:
  lons = [ lon+360 if lon < 0 else lon for lon in lons]
plt.plot(lons, lats, color='blue', linestyle='--',
         label="Nav Great Circle",
         transform=ccrs.PlateCarree())
                # straight lines link the coords

# show GC length on the map
dist = nav.distTo(lat1, lon1, lat2, lon2)
print(f"\nGreat Circle distance: {dist:.0f} km")
lat, lon, _ = coords[len(coords)//2]
plt.text(lon, lat, f" {dist:.0f}km",
         color='blue', transform=ccrs.PlateCarree())
            # no need for curves

# draw the Rhumb line trip as a series of coordinates
coords = nav.rhumbCoords(lat1, lon1, lat2, lon2, 20)
print("\nRhumb Line (lat,lon)s and deg Bearing:")
for coord in coords:
  cLat, cLon, cBear = coord
  print(f"  ({cLat:.2f}, {cLon:.2f}); {cBear:.2f}")

lats, lons, _ = zip(*coords)
if usesDL:
  lons = [ lon+360 if lon < 0 else lon for lon in lons]
plt.plot(lons, lats, color='red', ls='-.',
    label="Nav Rhumb Line", transform=ccrs.PlateCarree())
                # straight lines link the coords

# show Rhumb line length on the map
dist = nav.rhumbDistTo(lat1, lon1, lat2, lon2)
print(f"\nRhumb line distance: {dist:.0f} km")
lat, lon, _ = coords[len(coords)//2]
plt.text(lon, lat, f" {dist:.0f}km",
         color='red', transform=ccrs.PlateCarree())
                      # no need for curves

plotCity(ax, startCity, lat1, lon1)
plotCity(ax, endCity, lat2, lon2)
plt.title(f"{startCity} to {endCity}; projection: {mapType}\n")
plt.legend()
plt.show()