My Blog

#!/usr/bin/env python3

# Define function ...
def circumferenceAtLatitude(usrLat, /):
    """Find the circumference of the Earth at a given latitude

    This function finds the circumference (in metres) of the Earth at a given
    latitude (in degrees).

    Parameters
    ----------
    usrLat : float, numpy.ndarray
        the latitude (in degrees)

    Returns
    -------
    usrCir : float, numpy.ndarray
        the circumference of the Earth (in metres)

    Notes
    -----

    Copyright 2017 Thomas Guymer [1]_

    References
    ----------
    .. [1] PyGuymer3, https://github.com/Guymer/PyGuymer3
    """

    # Import special modules ...
    try:
        import numpy
    except:
        raise Exception("\"numpy\" is not installed; run \"pip install --user numpy\"") from None

    # Import my modules ...
    try:
        import pyguymer3
    except:
        raise Exception("\"pyguymer3\" is not installed; run \"pip install --user PyGuymer3\"") from None

    # Return answer ...
    return 2.0 * numpy.pi * pyguymer3.RADIUS_OF_EARTH * numpy.sin(numpy.radians(90.0 - usrLat))

# Define function ...
def degreesOfLongitude2Metres(usrLon, usrLat, /):
    """Find the length of a line of longitude at a given latitude

    This function finds the length (in metres) of a line of longitude (in
    degrees) on the surface of the Earth at a given latitude (in degrees).

    Parameters
    ----------
    usrLon : float
        the length of the line of longitude (in degrees)
    usrLat : float
        the latitude of the line of longitude (in degrees)

    Returns
    -------
    usrLen : float
        the length of the line of longitude (in metres)

    Notes
    -----

    Copyright 2017 Thomas Guymer [1]_

    References
    ----------
    .. [1] PyGuymer3, https://github.com/Guymer/PyGuymer3
    """

    # Return answer ...
    return (usrLon / 360.0) * circumferenceAtLatitude(usrLat)

# 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 os

    # Import special modules ...
    try:
        import cartopy
        cartopy.config.update(
            {
                "cache_dir" : os.path.expanduser("~/.local/share/cartopy_cache"),
            }
        )
    except:
        raise Exception("\"cartopy\" is not installed; run \"pip install --user Cartopy\"") from None
    try:
        import matplotlib
        matplotlib.rcParams.update(
            {
                       "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,
            }
        )
        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

    # Create latitudes and arrays to hold results ...
    lats = numpy.linspace(90.0, -90.0, num = 181)                               # [°]
    landLensDeg = numpy.zeros(lats.size, dtype = numpy.float64)                 # [°]
    landLensMet = numpy.zeros(lats.size, dtype = numpy.float64)                 # [m]

    # Create figure ...
    fg = matplotlib.pyplot.figure(figsize = (12.8, 7.2))

    # Create axis ...
    ax = pyguymer3.geo.add_axis(fg)

    # Configure axis ...
    pyguymer3.geo.add_map_background(ax, resolution = "large8192px")

    # Find file containing all the country shapes ...
    sfile = cartopy.io.shapereader.natural_earth(
          category = "cultural",
              name = "admin_0_countries",
        resolution = "10m",
    )

    # Loop over latitudes ...
    for i in range(lats.size):
        # Catch the special cases of the poles ...
        # NOTE: The North Pole has no land; the South Pole has land (hence a
        #       full 360° in Geodetic space, which corresponds to only 0m in
        #       Euclidean space).
        if i == 0:
            landLensDeg[i] = 0.0                                                # [°]
            landLensMet[i] = 0.0                                                # [m]
            continue
        if i == lats.size - 1:
            landLensDeg[i] = 360.0                                              # [°]
            landLensMet[i] = 0.0                                                # [m]
            continue

        # Make grid line ...
        # HACK: Some of the regions go past 180 degrees (presumably due to
        #       floating-point rounding errors) so the grid line has to wrap
        #       around to ensure that it crosses at the anti-meridian.
        grid_line = shapely.geometry.LineString(
            [
                (-181.0, lats[i]),
                ( 181.0, lats[i]),
            ],
        )

        # Plot grid line ...
        ax.plot(
            [-180.0, 180.0],
            [lats[i], lats[i]],
                color = "blue",
            linewidth = 0.1,
            transform = cartopy.crs.PlateCarree(),
        )

        # Make stack of longitudes ...
        lons = numpy.empty(0, dtype = numpy.float64)                            # [°]

        # Loop over records ...
        for record in cartopy.io.shapereader.Reader(sfile).records():
            # Loop over Polygons ...
            for poly in pyguymer3.geo.extract_polys(record.geometry):
                # Find the Points which intersect the exterior ring and the grid
                # line ...
                points = pyguymer3.geo.extract_points(poly.exterior.intersection(grid_line))

                # Skip this intersection if there aren't any Points ...
                if len(points) == 0:
                    continue

                # Check that a sensible number of Points have been found ...
                if len(points) % 2 == 1:
                    print(f"WARNING: An odd number of Points was found at {lats[i]:.1f} (n = {len(points):d})")

                    # Loop over Points ...
                    for point in points:
                        # Extract coordinates ...
                        lon, lat = point.coords[0]                              # [°], [°]

                        # Add location to plot ...
                        ax.plot(
                            lon,
                            lat,
                                      color = "green",
                                     marker = "o",
                            markeredgewidth = 0.0,
                                 markersize = 1.0,
                                  transform = cartopy.crs.PlateCarree(),
                        )
                else:
                    # Loop over Points ...
                    for point in points:
                        # Extract coordinates and add to stack ...
                        lon, lat = point.coords[0]                              # [°], [°]
                        lons = numpy.append(lons, lon)                          # [°]

        # Skip if there aren't any intersections ...
        if lons.size == 0:
            continue

        # Sort array ...
        lons.sort()

        # Loop over longitude pairs ...
        for j in range(0, lons.size, 2):
            # Add length to total ...
            landLensDeg[i] += (lons[j + 1] - lons[j])                           # [°]

            # Plot length ...
            ax.plot(
                [lons[j], lons[j + 1]],
                [lats[i], lats[i]],
                    color = "red",
                linewidth = 1.0,
                transform = cartopy.crs.PlateCarree(),
            )

    # Populate array ...
    for i in range(1, lats.size - 1):
        landLensMet[i] = degreesOfLongitude2Metres(landLensDeg[i], lats[i])     # [m]

    # Configure figure ...
    fg.tight_layout()

    # Save figure ...
    fg.savefig("land_vs_latitude_map.png")
    matplotlib.pyplot.close(fg)

    # Optimize PNG ...
    pyguymer3.image.optimize_image(
        "land_vs_latitude_map.png",
        strip = True,
    )

    # **************************************************************************

    # Save data as CSV ...
    with open("land_vs_latitude.csv", "wt", encoding = "utf-8") as fObj:
        fObj.write("latitude [°],length of land [°],length of land [m]\n")
        for i in range(lats.size):
            fObj.write(f"{lats[i]:.15e},{landLensDeg[i]:.15e},{landLensMet[i]:.15e}\n")

    # **************************************************************************

    # Create figure ...
    fg = matplotlib.pyplot.figure()

    # Create axis ...
    ax = fg.add_subplot()

    # Plot data ...
    ax.fill_between(
        lats,
        numpy.zeros(lats.size, dtype = numpy.float64),
        landLensMet / 1.0e3,
        color = "green",
    )
    ax.fill_between(
        lats,
        landLensMet / 1.0e3,
        circumferenceAtLatitude(lats) / 1.0e3,
        color = "blue",
    )

    # Draw labels ...
    for lat in [-70.0, -60.0, -50.0, -40.0]:
        ax.plot(
            [lat, lat,],
            [0.0, circumferenceAtLatitude(lat) / 1.0e3,],
            color = "red",
        )
    ax.text(
        -65.0,
        0.5 * circumferenceAtLatitude(-65.0) / 1.0e3,
        "Screaming Sixties",
                      color = "red",
        horizontalalignment = "center",
                   rotation = "vertical",
          verticalalignment = "center",
    )
    ax.text(
        -55.0,
        0.5 * circumferenceAtLatitude(-55.0) / 1.0e3,
        "Furious Fifties",
                      color = "red",
        horizontalalignment = "center",
                   rotation = "vertical",
          verticalalignment = "center",
    )
    ax.text(
        -45.0,
        0.5 * circumferenceAtLatitude(-45.0) / 1.0e3,
        "Roaring Forties",
                      color = "red",
        horizontalalignment = "center",
                   rotation = "vertical",
          verticalalignment = "center",
    )

    # Configure axis ...
    ax.grid()
    ax.set_xlabel("Latitude [°]")
    ax.set_xlim(-90.0, +90.0)
    ax.set_ylabel("Distance [km]")
    ax.set_ylim(0.0, 2.0 * numpy.pi * pyguymer3.RADIUS_OF_EARTH / 1.0e3)

    # Configure figure ...
    fg.tight_layout()

    # Save figure ...
    fg.savefig("land_vs_latitude_plot.png")
    matplotlib.pyplot.close(fg)

    # Optimize PNG ...
    pyguymer3.image.optimize_image(
        "land_vs_latitude_plot.png",
        strip = True,
    )