#!/usr/bin/env python3 # Use the proper idiom in the main module ... # NOTE: See https://docs.python.org/3.12/library/multiprocessing.html#the-spawn-and-forkserver-start-methods if __name__ == "__main__": # Import standard modules ... import pathlib # Import special modules ... try: import cartopy cartopy.config.update( { "cache_dir" : pathlib.PosixPath("~/.local/share/cartopy").expanduser(), } ) import cartopy.geodesic except: raise Exception("\"cartopy\" is not installed; run \"pip install --user Cartopy\"") from None try: import matplotlib matplotlib.rcParams.update( { "axes.xmargin" : 0.01, "axes.ymargin" : 0.01, "backend" : "Agg", # NOTE: See https://matplotlib.org/stable/gallery/user_interfaces/canvasagg.html "figure.dpi" : 300, "figure.figsize" : (9.6, 7.2), # NOTE: See https://github.com/Guymer/misc/blob/main/README.md#matplotlib-figure-sizes "font.size" : 8, "image.interpolation" : "none", # NOTE: See https://matplotlib.org/stable/gallery/images_contours_and_fields/interpolation_methods.html "image.resample" : False, } ) import matplotlib.pyplot except: raise Exception("\"matplotlib\" is not installed; run \"pip install --user matplotlib\"") from None try: import numpy except: raise Exception("\"numpy\" is not installed; run \"pip install --user numpy\"") from None try: import shapely import shapely.geometry except: raise Exception("\"shapely\" is not installed; run \"pip install --user Shapely\"") from None # Import my modules ... try: import pyguymer3 import pyguymer3.geo import pyguymer3.image except: raise Exception("\"pyguymer3\" is not installed; run \"pip install --user PyGuymer3\"") from None # ************************************************************************** # Define point ... lonA = -0.4619 # [°] latA = 51.4706 # [°] # Define point ... lonB = -122.3750 # [°] latB = 37.6180 # [°] # Define interpolation parameters ... npoint = 99 # ************************************************************************** # Create the Euclidean line ... lons = numpy.linspace(lonA, lonB, num = npoint) # [°] lats = numpy.linspace(latA, latB, num = npoint) # [°] points = [] for ipoint in range(npoint): points.append((lons[ipoint], lats[ipoint])) euclideanLine = shapely.geometry.linestring.LineString(points) euclideanCoords = numpy.array(euclideanLine.coords) # [°] # Create Geodesic line ... geodesicLine = pyguymer3.geo.great_circle( lonA, latA, lonB, latB, npoint = npoint, ) geodesicCoords = numpy.array(geodesicLine.coords) # [°] # ************************************************************************** # Find the Euclidean and Geodetic distances of the Euclidean line ... euclideanDist = 0.0 # [°] geodeticDist = 0.0 # [m] for ipoint in range(npoint - 1): euclideanDist += numpy.hypot( euclideanCoords[ipoint + 1, 0] - euclideanCoords[ipoint, 0], euclideanCoords[ipoint + 1, 1] - euclideanCoords[ipoint, 1], ) # [°] geodeticDist += pyguymer3.geo.calc_dist_between_two_locs( euclideanCoords[ipoint, 0], euclideanCoords[ipoint, 1], euclideanCoords[ipoint + 1, 0], euclideanCoords[ipoint + 1, 1], )[0] # [m] print(f"The shortest line in Euclidean space is {euclideanDist:,.1f}° long and {geodeticDist:,.1f}m long.") # Find the Euclidean and Geodetic distances of the Geodesic line ... euclideanDist = 0.0 # [°] geodeticDist = 0.0 # [m] for ipoint in range(npoint - 1): euclideanDist += numpy.hypot( geodesicCoords[ipoint + 1, 0] - geodesicCoords[ipoint, 0], geodesicCoords[ipoint + 1, 1] - geodesicCoords[ipoint, 1], ) # [°] geodeticDist += pyguymer3.geo.calc_dist_between_two_locs( geodesicCoords[ipoint, 0], geodesicCoords[ipoint, 1], geodesicCoords[ipoint + 1, 0], geodesicCoords[ipoint + 1, 1], )[0] # [m] print(f"The shortest line in Geodesic space is {euclideanDist:,.1f}° long and {geodeticDist:,.1f}m long.") # ************************************************************************** # Create figure ... fg = matplotlib.pyplot.figure() # Create axis ... ax1 = pyguymer3.geo.add_axis( fg, index = 1, ncols = 2, nrows = 2, ) # Create axis ... ax2 = pyguymer3.geo.add_axis( fg, index = 2, lat = geodesicCoords[(npoint - 1) // 2, 1], lon = geodesicCoords[(npoint - 1) // 2, 0], ncols = 2, nrows = 2, ) # Create axis ... ax3 = fg.add_subplot( 2, 2, (3, 4), ) # Configure axis ... pyguymer3.geo.add_map_background(ax1, subName = "large8192px") # Configure axis ... pyguymer3.geo.add_map_background(ax2, subName = "large8192px") # Configure axis ... ax3.grid() ax3.set_aspect("equal") ax3.set_xlabel("Longitude [°]") ax3.set_xlim(-180.0, +180.0) ax3.set_xticks(range(-180, 225, 45)) ax3.set_ylabel("Latitude [°]") ax3.set_ylim(-90.0, +90.0) ax3.set_yticks(range(-90, 135, 45)) # Plot Euclidean line thrice ... ax1.add_geometries( [euclideanLine], cartopy.crs.PlateCarree(), edgecolor = (1.0, 0.0, 0.0, 1.0), facecolor = "none", linewidth = 1.0, ) ax2.add_geometries( [euclideanLine], cartopy.crs.PlateCarree(), edgecolor = (1.0, 0.0, 0.0, 1.0), facecolor = "none", linewidth = 1.0, ) ax3.plot( euclideanCoords[:, 0], euclideanCoords[:, 1], color = (1.0, 0.0, 0.0, 1.0), linewidth = 1.0, ) # Plot Geodesic line thrice ... ax1.add_geometries( [geodesicLine], cartopy.crs.PlateCarree(), edgecolor = (0.0, 0.0, 1.0, 1.0), facecolor = "none", linewidth = 1.0, ) ax2.add_geometries( [geodesicLine], cartopy.crs.PlateCarree(), edgecolor = (0.0, 0.0, 1.0, 1.0), facecolor = "none", linewidth = 1.0, ) ax3.plot( geodesicCoords[:, 0], geodesicCoords[:, 1], color = (0.0, 0.0, 1.0, 1.0), linewidth = 1.0, ) # Configure figure ... fg.suptitle(f"Flight from ({lonA:.1f}°,{latA:.1f}°) to ({lonB:.1f}°,{latB:.1f}°)\nshortest Euclidean distance in red; shortest Geodetic distance in blue") fg.tight_layout() # Save figure ... fg.savefig("flight.png") matplotlib.pyplot.close(fg) # Optimise PNG ... pyguymer3.image.optimise_image("flight.png", strip = True)