Reliably Buffering A Shape By A Geodetic Distance

metadata

This post has been a long time coming and it is broken down into sections:

I have also included some discussion on the background principles associated with this topic, such as great circles and Tissot’s indicatrix.

§1 Introduction

On the 25th of February 2017 I committed my first implementation of Vincenty’s formulae into my Python 2 module. Three and a half years later, on the 21st of August 2020 I committed my second implementation of Vincenty’s formulae into my FORTRAN module. What is the significance of these events? Well …

Some of the Python scripts that I write would be categorised as “geospatial themed”. These scripts all share a common methodology: use Shapely to store the points/lines/polygons; use Cartopy to handle map projections; and use MatPlotLib to plot the output. Shapely is an excellent and powerful Python module to handle geometries, it is able to perform tasks such as buffering and unions, which I find very handy. Shapely is very generic: it just describes points/lines/polygons as a series of numbers without units. Here lies the problem: if you wanted to define a point using longitude and latitude, and then find the circle on the surface of the Earth that had a 1 kilometre radius around that point, then Shapely cannot do that for you (as that task would involve mixing units). If you wanted to find the circle on the surface of the Earth that had a 1 degree radius around that point, then Shapely can do that for you but I am not sure why anyone would want that.

§2 Euclidean Space -vs- Geodesic Space

One of the tasks that Shapely is able to perform is buffering, which is the process of enlarging an object by a constant distance. Unfortunately, Shapely has no knowledge of degrees or metres, and so the point in the earlier example would be enlarged in all directions by X [units that the point is defined by], in this case: X degrees, rather than X metres. The relationship between degrees and metres on the surface of the Earth changes depending on where you are and what direction you are going. In mathematics, we call these two different spaces Euclidean space and Geodesic space: in the earlier example Shapely buffers the point in Euclidean space by a constant distance (in degrees) in all directions, yet I want to buffer the point in Geodesic space by a constant distance (in metres).

§3 Great Circles

This difference is best highlighted by aeroplanes: aeroplanes fly in a straight line from point A to point B over the surface of the Earth. Depending on the map projection that is used, this line appears curved. These lines are called great circles. For example, consider the flight from London to San Francisco:

Download:
  1. 512 px × 384 px (0.2 Mpx; 120.4 KiB)
  2. 1,024 px × 768 px (0.8 Mpx; 424.1 KiB)
  3. 2,048 px × 1,536 px (3.1 Mpx; 1.5 MiB)
  4. 2,880 px × 2,160 px (6.2 Mpx; 2.6 MiB)

Here is the Python script to create the above map:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
 

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

    # Import special modules ...
    try:
        import cartopy
        cartopy.config.update(
            {
                "cache_dir" : os.path.expanduser("~/.local/share/cartopy_cache"),
            }
        )
        import cartopy.geodesic
    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

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

    # 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, resolution = "large8192px")

    # Configure axis ...
    pyguymer3.geo.add_map_background(ax2, resolution = "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.optimize_image("flight.png", strip = True)
You may also download “flight.py” directly or view “flight.py” on GitHub Gist (you may need to manually checkout the “main” branch).

In the above map, the red line is the line connecting (-0.4619°,51.4706°) and (-122.3750°,37.6180°) in Euclidean space, i.e., treating all degrees as being of equal size in all directions. The blue line is the line connecting (-0.4619°,51.4706°) and (-122.3750°,37.6180°) in Geodesic space, i.e., having knowledge that the degrees actually describe a location on the surface of a spheroid and that degrees have different physical size depending on where you are and which direction you are pointing. The above script also reports that:

This demonstrates why aeroplanes fly the large arcs that they do: the red line is 1,148,588.3 metres shorter than the blue line, but it is also 9.4° longer. As a passenger on an aeroplane, which one of those two differences is more important to you?

§4 Tissot’s Indicatrix

In the above map, the red and blue lines appear differently depending on the map projection that is used. No map projection is perfect, just see xkcd for some hilarious examples. As a way of displaying how distorted, or not, a map projection is Tissot’s indicatrix are sometimes displayed on a map. These are simply circles which have a constant radius in metres and which, therefore, appear distorted on the surface of the Earth. Cartopy has an example script which displays Tissot’s indicatrix and produces the below map.

Download:
  1. 512 px × 256 px (0.1 Mpx; 126.8 KiB)
  2. 1,000 px × 500 px (0.5 Mpx; 375.2 KiB)

Here is the Python script to create the above map:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
 

#!/usr/bin/env python3

"""
Tissot's Indicatrix
-------------------

Visualize Tissot's indicatrix on a map.

"""
import matplotlib.pyplot as plt

import cartopy.crs as ccrs


def main():
    fg = plt.figure(figsize = (10, 5))
    ax = fg.add_subplot(1, 1, 1, projection = ccrs.PlateCarree())

    # make the map global rather than have it zoom in to
    # the extents of any plotted data
    ax.set_global()

    ax.stock_img()
    ax.coastlines()

    ax.tissot(facecolor = "orange", alpha = 0.4)

    plt.show()


if __name__ == "__main__":
    main()
You may also download “tissot.py” directly or view “tissot.py” on GitHub Gist (you may need to manually checkout the “main” branch).

§5 Vincenty’s Formulae

… and now all the way back to me mentioning Vincenty’s formulae in the introduction. Vincenty’s formulae are iterating functions which perform two tasks:

If you wanted to create your own Tissot’s indicatrix manually then all you would need to do is pick a starting location and call “the direct problem” of Vincenty’s formulae for a list of bearings (with a constant distance) to obtain a list of end locations. These end locations are the circumference of a circle used in Tissot’s indicatrix. You can then repeat the process for different starting locations to populate the whole map projection.

If you look at the source code of the tissot() method you will see that it calls the circle() method in the Geodesic class. If you look at the source code of the Geodesic class then you will see that Cartopy just calls its own implementation of “the direct problem”.

§6 Buffering

The below map shows Lisbon buffered by 1,000.0 kilometres by both Cartopy and my Python 3 module. Reassuringly, they both give the same answer.

Download:
  1. 512 px × 384 px (0.2 Mpx; 117.2 KiB)
  2. 1,024 px × 768 px (0.8 Mpx; 418.1 KiB)
  3. 2,048 px × 1,536 px (3.1 Mpx; 1.5 MiB)
  4. 2,880 px × 2,160 px (6.2 Mpx; 2.6 MiB)

Here is the Python script to create the above map:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
 

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

    # Import special modules ...
    try:
        import cartopy
        cartopy.config.update(
            {
                "cache_dir" : os.path.expanduser("~/.local/share/cartopy_cache"),
            }
        )
        import cartopy.geodesic
    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

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

    # Define point ...
    lon = -9.144020                                                             # [°]
    lat = 38.719144                                                             # [°]

    # Define buffering parameters ...
    dist = 1000.0e3                                                             # [m]
    nAng = 9

    # Create an array based off the scalars ...
    pnt = numpy.zeros((1, 2), dtype = numpy.float64)                            # [°]
    pnt[0, 0] = lon                                                             # [°]
    pnt[0, 1] = lat                                                             # [°]

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

    # Buffer point using Cartopy ...
    # NOTE: Cartopy does not close the loop, so I must use one fewer points
    #       during the buffer and then manually add the 12 o'clock position on
    #       the end.
    # NOTE: Cartopy just passes the list of point directly to Shapely when
    #       making a Tissot circle, see:
    #         * https://scitools.org.uk/cartopy/docs/latest/_modules/cartopy/mpl/geoaxes.html#GeoAxes.tissot
    ring1 = cartopy.geodesic.Geodesic().circle(
        lon,
        lat,
        dist,
        n_samples = nAng - 1,
    )                                                                           # [°]
    ring1 = numpy.concatenate([ring1, [ring1[0, :]]])                           # [°]
    poly1 = shapely.geometry.Polygon(ring1)

    # Buffer point using PyGuymer3 ...
    ring2 = pyguymer3.geo._buffer_points_crudely(                               # pylint: disable=W0212
        pnt,
        dist,
        nAng,
    )[0, :, :]                                                                  # [°]
    poly2 = shapely.geometry.Polygon(ring2)

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

    # 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 = lat,
          lon = lon,
        ncols = 2,
        nrows = 2,
    )

    # Create axis ...
    ax3 = fg.add_subplot(
        2,
        2,
        (3, 4),
    )

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

    # Configure axis ...
    pyguymer3.geo.add_map_background(ax2, resolution = "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 Cartopy buffer thrice ...
    ax1.plot(
        ring1[:, 0],
        ring1[:, 1],
            color = (1.0, 0.0, 0.0, 1.0),
        linestyle = "--",
           marker = "o",
        transform = cartopy.crs.PlateCarree(),
    )
    ax1.add_geometries(
        [poly1],
        cartopy.crs.PlateCarree(),
        edgecolor = "none",
        facecolor = (1.0, 0.0, 0.0, 0.5),
        linewidth = 1.0,
    )
    ax2.plot(
        ring1[:, 0],
        ring1[:, 1],
            color = (1.0, 0.0, 0.0, 1.0),
        linestyle = "--",
           marker = "o",
        transform = cartopy.crs.PlateCarree(),
    )
    ax2.add_geometries(
        [poly1],
        cartopy.crs.PlateCarree(),
        edgecolor = "none",
        facecolor = (1.0, 0.0, 0.0, 0.5),
        linewidth = 1.0,
    )
    ax3.plot(
        ring1[:, 0],
        ring1[:, 1],
            color = (1.0, 0.0, 0.0, 1.0),
        linestyle = "--",
           marker = "o",
    )

    # Plot PyGuymer3 buffer thrice ...
    ax1.plot(
        ring2[:, 0],
        ring2[:, 1],
            color = (0.0, 0.0, 1.0, 1.0),
        transform = cartopy.crs.PlateCarree(),
    )
    ax1.add_geometries(
        [poly2],
        cartopy.crs.PlateCarree(),
        edgecolor = "none",
        facecolor = (0.0, 0.0, 1.0, 0.5),
        linewidth = 1.0,
    )
    ax2.plot(
        ring2[:, 0],
        ring2[:, 1],
            color = (0.0, 0.0, 1.0, 1.0),
        transform = cartopy.crs.PlateCarree(),
    )
    ax2.add_geometries(
        [poly2],
        cartopy.crs.PlateCarree(),
        edgecolor = "none",
        facecolor = (0.0, 0.0, 1.0, 0.5),
        linewidth = 1.0,
    )
    ax3.plot(
        ring2[:, 0],
        ring2[:, 1],
        color = (0.0, 0.0, 1.0, 1.0),
    )

    # Configure figure ...
    fg.suptitle(f"({lon:.1f}°,{lat:.1f}°) buffered by {0.001 * dist:,.1f} km\nCartopy in red; PyGuymer3 in blue")
    fg.tight_layout()

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

    # Optimise PNG ...
    pyguymer3.image.optimize_image("lisbon.png", strip = True)
You may also download “lisbon.py” directly or view “lisbon.py” on GitHub Gist (you may need to manually checkout the “main” branch).

Now, if I try the same script but using a random point in the middle of the Pacific Ocean then things go horribly wrong.

Download:
  1. 512 px × 384 px (0.2 Mpx; 110.6 KiB)
  2. 1,024 px × 768 px (0.8 Mpx; 374.3 KiB)
  3. 2,048 px × 1,536 px (3.1 Mpx; 1.3 MiB)
  4. 2,880 px × 2,160 px (6.2 Mpx; 2.1 MiB)

Here is the Python script to create the above map:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
 

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

    # Import special modules ...
    try:
        import cartopy
        cartopy.config.update(
            {
                "cache_dir" : os.path.expanduser("~/.local/share/cartopy_cache"),
            }
        )
        import cartopy.geodesic
    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

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

    # Define point ...
    lon = 170.0                                                                 # [°]
    lat = 10.0                                                                  # [°]

    # Define buffering parameters ...
    dist = 4000.0e3                                                             # [m]
    nAng = 9

    # Create an array based off the scalars ...
    pnt = numpy.zeros((1, 2), dtype = numpy.float64)                            # [°]
    pnt[0, 0] = lon                                                             # [°]
    pnt[0, 1] = lat                                                             # [°]

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

    # Buffer point using Cartopy ...
    # NOTE: Cartopy does not close the loop, so I must use one fewer points
    #       during the buffer and then manually add the 12 o'clock position on
    #       the end.
    # NOTE: Cartopy just passes the list of point directly to Shapely when
    #       making a Tissot circle, see:
    #         * https://scitools.org.uk/cartopy/docs/latest/_modules/cartopy/mpl/geoaxes.html#GeoAxes.tissot
    ring1 = cartopy.geodesic.Geodesic().circle(
        lon,
        lat,
        dist,
        n_samples = nAng - 1,
    )                                                                           # [°]
    ring1 = numpy.concatenate([ring1, [ring1[0, :]]])                           # [°]
    poly1 = shapely.geometry.Polygon(ring1)

    # Buffer point using PyGuymer3 ...
    ring2 = pyguymer3.geo._buffer_points_crudely(                               # pylint: disable=W0212
        pnt,
        dist,
        nAng,
    )[0, :, :]                                                                  # [°]
    poly2 = shapely.geometry.Polygon(ring2)

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

    # 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 = lat,
          lon = lon,
        ncols = 2,
        nrows = 2,
    )

    # Create axis ...
    ax3 = fg.add_subplot(
        2,
        2,
        (3, 4),
    )

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

    # Configure axis ...
    pyguymer3.geo.add_map_background(ax2, resolution = "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 Cartopy buffer thrice ...
    ax1.plot(
        ring1[:, 0],
        ring1[:, 1],
            color = (1.0, 0.0, 0.0, 1.0),
        linestyle = "--",
           marker = "o",
        transform = cartopy.crs.PlateCarree(),
    )
    ax1.add_geometries(
        [poly1],
        cartopy.crs.PlateCarree(),
        edgecolor = "none",
        facecolor = (1.0, 0.0, 0.0, 0.5),
        linewidth = 1.0,
    )
    ax2.plot(
        ring1[:, 0],
        ring1[:, 1],
            color = (1.0, 0.0, 0.0, 1.0),
        linestyle = "--",
           marker = "o",
        transform = cartopy.crs.PlateCarree(),
    )
    ax2.add_geometries(
        [poly1],
        cartopy.crs.PlateCarree(),
        edgecolor = "none",
        facecolor = (1.0, 0.0, 0.0, 0.5),
        linewidth = 1.0,
    )
    ax3.plot(
        ring1[:, 0],
        ring1[:, 1],
            color = (1.0, 0.0, 0.0, 1.0),
        linestyle = "--",
           marker = "o",
    )

    # Plot PyGuymer3 buffer thrice ...
    ax1.plot(
        ring2[:, 0],
        ring2[:, 1],
            color = (0.0, 0.0, 1.0, 1.0),
        transform = cartopy.crs.PlateCarree(),
    )
    ax1.add_geometries(
        [poly2],
        cartopy.crs.PlateCarree(),
        edgecolor = "none",
        facecolor = (0.0, 0.0, 1.0, 0.5),
        linewidth = 1.0,
    )
    ax2.plot(
        ring2[:, 0],
        ring2[:, 1],
            color = (0.0, 0.0, 1.0, 1.0),
        transform = cartopy.crs.PlateCarree(),
    )
    ax2.add_geometries(
        [poly2],
        cartopy.crs.PlateCarree(),
        edgecolor = "none",
        facecolor = (0.0, 0.0, 1.0, 0.5),
        linewidth = 1.0,
    )
    ax3.plot(
        ring2[:, 0],
        ring2[:, 1],
        color = (0.0, 0.0, 1.0, 1.0),
    )

    # Configure figure ...
    fg.suptitle(f"({lon:.1f}°,{lat:.1f}°) buffered by {0.001 * dist:,.1f} km\nCartopy in red; PyGuymer3 in blue")
    fg.tight_layout()

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

    # Optimise PNG ...
    pyguymer3.image.optimize_image("pacific.png", strip = True)
You may also download “pacific.py” directly or view “pacific.py” on GitHub Gist (you may need to manually checkout the “main” branch).

Shapely tries to make a polygon from the list of points along the ring, but it has no idea about the anti-meridian and so whilst it does make a mathematically correct polygon it is certainly not the one that the user desired. In summary, Vincenty’s formulae correctly return the longitudes and latitudes of all of the points along a ring which is a fixed distance (in metres) around a point but assembling these points into a polygon in Shapely is not trivial when the ring crosses either the anti-meridian or one of the poles.

Therefore, a function needs writing which takes the raw list of longitude and latitude pairs and uses them to create a polygon for the user. This function needs to be clever enough to return two individual polygons when the ring around a point crosses either the anti-meridian or one of the poles. I wrote a function to perform this task and it served me well for numerous years, such as:

However, it was not completely robust and as I worked on ever more complex projects I found sets of input values which would break it. I believe that my new function pyguymer3.geo._points2polys() works more robustly than any of the other ones that I have written over the years. For example, see how it fixes the buffered point in the Pacific Ocean (from above) by returning two unconnected polygons:

Download:
  1. 512 px × 384 px (0.2 Mpx; 113.1 KiB)
  2. 1,024 px × 768 px (0.8 Mpx; 399.9 KiB)
  3. 2,048 px × 1,536 px (3.1 Mpx; 1.4 MiB)
  4. 2,880 px × 2,160 px (6.2 Mpx; 2.5 MiB)

Here is the Python script to create the above map:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
 

#!/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 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

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

    # Define point ...
    lon = 170.0                                                                 # [°]
    lat = 10.0                                                                  # [°]

    # Define buffering parameters ...
    dist = 4000.0e3                                                             # [m]
    nAng = 9

    # Create a point based off the scalars ...
    pnt = shapely.geometry.point.Point(lon, lat)

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

    # Buffer point using PyGuymer3 ...
    multipoly = pyguymer3.geo.buffer(
        pnt,
        dist,
        fill = 1.0,
        nAng = nAng,
        simp = -1.0,
    )

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

    # 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 = lat,
          lon = lon,
        ncols = 2,
        nrows = 2,
    )

    # Create axis ...
    ax3 = fg.add_subplot(
        2,
        2,
        (3, 4),
    )

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

    # Configure axis ...
    pyguymer3.geo.add_map_background(
        ax2,
        resolution = "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 PyGuymer3 buffer thrice ...
    ax1.add_geometries(
        pyguymer3.geo.extract_polys(multipoly),
        cartopy.crs.PlateCarree(),
        edgecolor = (0.0, 1.0, 0.0, 1.0),
        facecolor = (0.0, 1.0, 0.0, 0.5),
        linewidth = 1.0,
    )
    ax2.add_geometries(
        pyguymer3.geo.extract_polys(multipoly),
        cartopy.crs.PlateCarree(),
        edgecolor = (0.0, 1.0, 0.0, 1.0),
        facecolor = (0.0, 1.0, 0.0, 0.5),
        linewidth = 1.0,
    )
    for poly in pyguymer3.geo.extract_polys(multipoly):
        coords = numpy.array(poly.exterior.coords)                              # [°]
        ax3.plot(
            coords[:, 0],
            coords[:, 1],
            color = (0.0, 1.0, 0.0, 1.0),
        )

    # Configure figure ...
    fg.suptitle(f"({lon:.1f}°,{lat:.1f}°) buffered by {0.001 * dist:,.1f} km\nPyGuymer3 in green")
    fg.tight_layout()

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

    # Optimise PNG ...
    pyguymer3.image.optimize_image("working.png", strip = True)
You may also download “working.py” directly or view “working.py” on GitHub Gist (you may need to manually checkout the “main” branch).

§7 Testing

my Python 3 module contains a couple of test scripts (“bufferPoint” and “buffer”) which have slowly ballooned to contain new inputs which demonstrate previous shortcomings and which hope to guard against future regressions in robustness. After growing tired of adding yet more examples to these scripts, I decided to write a third script which just loops over all combinations of longitude and latitude on Earth in an effort to pre-emptively catch further errors: “animateBufferPoint”. I am pleased to include the video below, which has been generated by this testing script of my new function:

Download:
  1. 512 px × 512 px (0.3 Mpx; 29.3 MiB)
  2. 1,024 px × 1,024 px (1.0 Mpx; 59.1 MiB)
  3. 1,800 px × 1,800 px (3.2 Mpx; 130.3 MiB)

As you can see, across all locations on the surface of the Earth, the function is able to take the longitude and latitude pairs and correctly generate one or two Shapely polygons which correctly describe what the user is expecting: a perfect circle around a point with a fixed radius, defined in metres. Whether this point is near to one of the Poles or near to the anti-meridian is irrelevant: when viewed from an orthographic projection above the point then the circle appears as a circle. There are no glitches or weird bumps or excursions in the polygon(s).

Buffering a polygon is similarly straightforward: just split the polygon up into its individual points (on its exterior and its interiors) then return the union of these buffers with the original polygon.

§8 Summary

In summary, after over five years of writing “geospatial themed” Python scripts, I believe that I now have a robust Python function which can return one or two Shapely polygons of the buffer around a point on the surface of the Earth. With the tests scripts that are part of my my Python 3 module module, I believe that this function (and its wrappers) works correctly for all locations on the surface of the Earth.