My Blog

#!/usr/bin/env python3

# Use the proper idiom in the main module ...
# NOTE: See https://docs.python.org/3.11/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; you need to have the Python module from https://github.com/Guymer/PyGuymer3 located somewhere in your $PYTHONPATH") from None

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

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

    # Create axis ...
    ax = fg.add_subplot(projection = cartopy.crs.Robinson())

    # Configure axis ...
    ax.set_global()
    pyguymer3.geo.add_coastlines(ax)
    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 ...
        if i == 0:
            lands[i] = 0.0                                                      # [°]
            continue
        if i == lats.size - 1:
            lands[i] = 360.0                                                    # [°]
            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 ...
            lands[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(),
            )

    # Configure figure ...
    fg.tight_layout()

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

    # Optimize PNG ...
    pyguymer3.image.optimize_image(
        "land_vs_latitude.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 [°]\n")
        for i in range(lats.size):
            fObj.write(f"{lats[i]:e},{lands[i]:e}\n")

              

# Specify output settings ...
set output "land_vs_latitude.svg"
set terminal svg size 750, 550 enhanced mouse jsdir "/js" rounded dashed

# Make functions ...
radius = 6371.009
radians(x) = pi * x / 180.0
f(x) = 2.0 * pi * radius * sin(radians(90.0 - x))
g(x,y) = (y / 360.0) * f(x)

# Specify plot settings ...
unset key
set grid
set title "Land -vs- Water (On Planet Earth)"
set xlabel "Latitude [deg]"
set xrange [-90:90]
set ylabel "Distance [km]"
set yrange [0:2.0 * pi * radius]

# Draw labels ...
set arrow from -40, 0 to -40, f(-40) nohead front lc rgb "red"
set label "Roaring Forties"   at -45, 0.5 * f(-45) center rotate by 90 front textcolor rgb "red"
set arrow from -50, 0 to -50, f(-50) nohead front lc rgb "red"
set label "Furious Fifties"   at -55, 0.5 * f(-55) center rotate by 90 front textcolor rgb "red"
set arrow from -60, 0 to -60, f(-60) nohead front lc rgb "red"
set label "Screaming Sixties" at -65, 0.5 * f(-65) center rotate by 90 front textcolor rgb "red"
set arrow from -70, 0 to -70, f(-70) nohead front lc rgb "red"

# Plot the data ...
set datafile separator ","
plot "land_vs_latitude.csv" every ::1 using 1:(0):(g($1, $2))     with filledcurves lc rgb "green", \
     "land_vs_latitude.csv" every ::1 using 1:(g($1, $2)):(f($1)) with filledcurves lc rgb "blue"