Coverage for pyguymer3/geo/create_map_of_points.py: 1%
87 statements
« prev ^ index » next coverage.py v7.9.2, created at 2025-07-08 18:47 +0000
1#!/usr/bin/env python3
3# Define function ...
4def create_map_of_points(
5 pntLons,
6 pntLats,
7 pngOut,
8 /,
9 *,
10 angConv = 0.1,
11 background = "NE",
12 chunksize = 1048576,
13 conv = 1.0e3,
14 debug = __debug__,
15 eps = 1.0e-12,
16 exiftoolPath = None,
17 extent = None,
18 fillColor = (255.0 / 255.0, 0.0 / 255.0, 0.0 / 255.0),
19 gifsiclePath = None,
20 jpegtranPath = None,
21 method = "GeodesicBox",
22 name = "natural-earth-1",
23 nAng = 9,
24 nIter = 100,
25 nRefine = 1,
26 onlyValid = False,
27 optipngPath = None,
28 padDist = 12.0 * 1852.0,
29 prefix = ".",
30 ramLimit = 1073741824,
31 repair = False,
32 resolution = "10m",
33 route = None,
34 routeFillColor = ( 0.0 / 255.0, 128.0 / 255.0, 0.0 / 255.0),
35 satellite_height = False,
36 scale = 1,
37 skipFillColor = (255.0 / 255.0, 165.0 / 255.0, 0.0 / 255.0),
38 skips = None,
39 timeout = 60.0,
40 title = None,
41 tol = 1.0e-10,
42 useSciPy = False,
43):
44 """Save a PNG map of a sequence of points
46 This function accepts a sequence of longitudes and latitudes then saves a
47 PNG map containing all of them drawn together in a big line.
49 Parameters
50 ----------
51 pntLons : numpy.ndarray
52 the sequence of longitudes
53 pntLats : numpy.ndarray
54 the sequence of latitudes
55 pngOut : str
56 the name of the output PNG
57 angConv : float, optional
58 the angle change which classifies as converged (in degrees)
59 background : str, optional
60 the type of background to add (recognised values are: "GSHHG"; "image";
61 "NE"; "none"; and "OSM")
62 chunksize : int, optional
63 the size of the chunks of any files which are read in (in bytes)
64 conv : float, optional
65 the Geodesic distance that defines the middle as being converged (in
66 metres)
67 debug : bool, optional
68 print debug messages and draw the circle on the axis
69 eps : float, optional
70 the tolerance of the Vincenty formula iterations
71 exiftoolPath : str, optional
72 the path to the "exiftool" binary (if not provided then Python will attempt to
73 find the binary itself)
74 extent : list of floats
75 for high-resolution images, save time by specifying the extent that is
76 to be added
77 fillColor : tuple of int, optional
78 the fill colour of the points
79 gifsiclePath : str, optional
80 the path to the "gifsicle" binary (if not provided then Python will attempt to
81 find the binary itself)
82 jpegtranPath : str, optional
83 the path to the "jpegtran" binary (if not provided then Python will attempt to
84 find the binary itself)
85 method : str, optional
86 the method for finding the middle of the points
87 name : str, optional
88 the name of the image in the database
89 nAng : int, optional
90 the number of angles around the middle location to search over
91 nIter : int, optional
92 the maximum number of iterations (particularly the Vincenty formula)
93 nRefine : int, optional
94 the number of refinements to make (each refinement halves the "conv"
95 distance)
96 onlyValid : bool, optional
97 only return valid Polygons (checks for validity can take a while, if
98 being called often)
99 optipngPath : str, optional
100 the path to the "optipng" binary (if not provided then Python will attempt to
101 find the binary itself)
102 padDist : float, optional
103 the padding to draw around the points (in metres)
104 prefix : str, optional
105 change the name of the output debugging CSVs
106 ramLimit : int, optional
107 the maximum RAM usage of each "large" array (in bytes)
108 repair : bool, optional
109 attempt to repair invalid Polygons
110 resolution : str, optional
111 the resolution of the image or NE dataset or GSHHG dataset
112 route : shapely.geometry.linestring.LineString, optional
113 an extra line to draw on the map
114 routeFillColor : tuple of int, optional
115 the fill colour of the extra route
116 satellite_height : float, optional
117 if a distance is provided then use a "NearsidePerspective" projection at
118 an altitude which has the same field-of-view as the distance
119 scale : int, optional
120 the scale of the tiles
121 skipFillColor : tuple of int, optional
122 the fill colour of the skipped points
123 skips : numpy.ndarray, optional
124 an array of booleans as to whether to include/exclude each individual
125 point from calculating the image's field-of-view (this allows the great
126 circles from flights to be drawn but for them to not expand the image to
127 fit in the departing airport); if not provided then all points are used
128 timeout : float, optional
129 the timeout for any requests/subprocess calls (in seconds)
130 title : str, optional
131 the title
132 tol : float, optional
133 the Euclidean distance that defines two points as being the same (in
134 degrees)
135 useSciPy : bool, optional
136 use "scipy.optimize.minimize" or my own minimizer
138 Notes
139 -----
140 Copyright 2017 Thomas Guymer [1]_
142 References
143 ----------
144 .. [1] PyGuymer3, https://github.com/Guymer/PyGuymer3
145 """
147 # Import standard modules ...
148 import os
150 # Import special modules ...
151 try:
152 import cartopy
153 cartopy.config.update(
154 {
155 "cache_dir" : os.path.expanduser("~/.local/share/cartopy_cache"),
156 }
157 )
158 except:
159 raise Exception("\"cartopy\" is not installed; run \"pip install --user Cartopy\"") from None
160 try:
161 import matplotlib
162 matplotlib.rcParams.update(
163 {
164 "backend" : "Agg", # NOTE: See https://matplotlib.org/stable/gallery/user_interfaces/canvasagg.html
165 "figure.dpi" : 300,
166 "figure.figsize" : (9.6, 7.2), # NOTE: See https://github.com/Guymer/misc/blob/main/README.md#matplotlib-figure-sizes
167 "font.size" : 8,
168 }
169 )
170 import matplotlib.pyplot
171 except:
172 raise Exception("\"matplotlib\" is not installed; run \"pip install --user matplotlib\"") from None
173 try:
174 import numpy
175 except:
176 raise Exception("\"numpy\" is not installed; run \"pip install --user numpy\"") from None
178 # Import sub-functions ...
179 from .add_axis import add_axis
180 from .add_GSHHG_map_underlay import add_GSHHG_map_underlay
181 from .add_map_background import add_map_background
182 from .add_NE_map_underlay import add_NE_map_underlay
183 from .add_OSM_map_background import add_OSM_map_background
184 from .extract_lines import extract_lines
185 from .find_middle_of_locs import find_middle_of_locs
186 from .great_circle import great_circle
187 from ..consts import RESOLUTION_OF_EARTH
188 from ..image import optimise_image
190 # **************************************************************************
192 # Check inputs ...
193 if skips is None:
194 skips = numpy.zeros(pntLons.size, dtype = bool)
196 # **************************************************************************
198 # Create figure ...
199 fg = matplotlib.pyplot.figure(figsize = (7.2, 7.2))
201 # **************************************************************************
203 # Find the centre of the points very quickly ...
204 midLonQuick, midLatQuick, maxDistQuick = find_middle_of_locs(
205 pntLons[numpy.logical_not(skips)],
206 pntLats[numpy.logical_not(skips)],
207 angConv = None,
208 conv = None,
209 debug = debug,
210 eps = eps,
211 method = "EuclideanBox",
212 midLat = None,
213 midLon = None,
214 nAng = None,
215 nIter = nIter,
216 nRefine = nRefine,
217 pad = -1.0,
218 useSciPy = None,
219 ) # [°]
221 # Check if the points are so widely spread that the map has to have global
222 # extent to show them all ...
223 if maxDistQuick > 90.0:
224 # Create axis ...
225 ax = add_axis(
226 fg,
227 add_coastlines = False,
228 add_gridlines = True,
229 configureAgain = bool(background == "OSM"),
230 debug = debug,
231 eps = eps,
232 gs = None,
233 index = None,
234 ncols = None,
235 nIter = nIter,
236 nrows = None,
237 onlyValid = onlyValid,
238 prefix = prefix,
239 ramLimit = ramLimit,
240 repair = repair,
241 tol = tol,
242 )
243 else:
244 # If the user asked for a Euclidean method then the padding distance
245 # needs converting from metres in to degrees ...
246 match method:
247 case "EuclideanBox" | "EuclideanCircle":
248 padDist /= RESOLUTION_OF_EARTH # [°]
249 case "GeodesicBox" | "GeodesicCircle":
250 pass
251 case _:
252 # Crash ...
253 raise ValueError(f"\"method\" is an unexpected value ({repr(method)})") from None
255 # Find the centre of the points ...
256 midLon, midLat, maxDist = find_middle_of_locs(
257 pntLons[numpy.logical_not(skips)],
258 pntLats[numpy.logical_not(skips)],
259 angConv = angConv,
260 conv = conv,
261 debug = debug,
262 eps = eps,
263 method = method,
264 midLat = midLatQuick,
265 midLon = midLonQuick,
266 nAng = nAng,
267 nIter = nIter,
268 nRefine = nRefine,
269 pad = padDist,
270 useSciPy = useSciPy,
271 ) # [°], [°], [°] or [m]
273 # Check what method the user wants ...
274 match method:
275 case "EuclideanBox" | "EuclideanCircle":
276 if debug:
277 print(f"INFO: Centre at (lon={midLon:+.6f}°, lat={midLat:+.6f}°) with a {maxDist:.6f}° radius.")
278 case "GeodesicBox" | "GeodesicCircle":
279 if debug:
280 print(f"INFO: Centre at (lon={midLon:+.6f}°, lat={midLat:+.6f}°) with a {0.001 * maxDist:,.1f} km radius.")
281 case _:
282 # Crash ...
283 raise ValueError(f"\"method\" is an unexpected value ({repr(method)})") from None
285 # If the user asked for a Euclidean method then the maximum distance
286 # needs converting from degrees in to metres ...
287 match method:
288 case "EuclideanBox" | "EuclideanCircle":
289 maxDist *= RESOLUTION_OF_EARTH # [m]
290 if debug:
291 print(f"INFO: Centre at (lon={midLon:+.6f}°, lat={midLat:+.6f}°) with a {0.001 * maxDist:,.1f} km radius.")
292 case "GeodesicBox" | "GeodesicCircle":
293 pass
294 case _:
295 # Crash ...
296 raise ValueError(f"\"method\" is an unexpected value ({repr(method)})") from None
298 # Create axis ...
299 ax = add_axis(
300 fg,
301 add_coastlines = False,
302 add_gridlines = True,
303 configureAgain = bool(background == "OSM"),
304 debug = debug,
305 dist = maxDist,
306 eps = eps,
307 gs = None,
308 index = None,
309 lat = midLat,
310 lon = midLon,
311 ncols = None,
312 nIter = nIter,
313 nrows = None,
314 onlyValid = onlyValid,
315 prefix = prefix,
316 ramLimit = ramLimit,
317 repair = repair,
318 satellite_height = satellite_height,
319 tol = tol,
320 )
322 # **************************************************************************
324 # Check which background the user wants ...
325 match background:
326 case "GSHHG":
327 # Add GSHHG background ...
328 add_GSHHG_map_underlay(
329 ax,
330 background = True,
331 debug = debug,
332 iceOcean = True,
333 islandLake = True,
334 lakeLand = True,
335 landOcean = True,
336 linewidth = 0.5,
337 onlyValid = onlyValid,
338 pondIsland = True,
339 repair = repair,
340 resolution = resolution,
341 )
342 case "image":
343 # Add image background ...
344 add_map_background(
345 ax,
346 debug = debug,
347 extent = extent,
348 name = name,
349 resolution = resolution,
350 )
351 case "NE":
352 # Add NE background ...
353 add_NE_map_underlay(
354 ax,
355 background = True,
356 cultural = True,
357 debug = debug,
358 linestyle = "solid",
359 linewidth = 0.5,
360 maxElev = 8850.0,
361 onlyValid = onlyValid,
362 physical = True,
363 repair = repair,
364 resolution = resolution,
365 )
366 case "none":
367 # Don't add any background ...
368 pass
369 case "OSM":
370 # Calculate the resolution depending on the half-height of the
371 # figure and the resolution of the figure ...
372 res = maxDist / (0.5 * fg.get_size_inches()[1] * fg.dpi) # [m/px]
374 # Add OpenStreetMap background ...
375 add_OSM_map_background(
376 ax,
377 midLat,
378 res,
379 debug = debug,
380 scale = scale,
381 )
382 case _:
383 # Crash ...
384 raise ValueError(f"\"background\" is an unexpected value ({repr(background)})") from None
386 # Plot locations ...
387 # NOTE: As of 5/Dec/2023, the default "zorder" of the coastlines is 1.5, the
388 # default "zorder" of the gridlines is 2.0 and the default "zorder" of
389 # the scattered points is 1.0.
390 ax.scatter(
391 pntLons[numpy.logical_not(skips)],
392 pntLats[numpy.logical_not(skips)],
393 edgecolor = "none",
394 facecolor = fillColor,
395 linewidth = 0.1,
396 s = 64.0,
397 transform = cartopy.crs.Geodetic(),
398 zorder = 5.0,
399 )
401 # Plot locations ...
402 # NOTE: As of 5/Dec/2023, the default "zorder" of the coastlines is 1.5, the
403 # default "zorder" of the gridlines is 2.0 and the default "zorder" of
404 # the scattered points is 1.0.
405 ax.scatter(
406 pntLons[skips],
407 pntLats[skips],
408 edgecolor = "none",
409 facecolor = skipFillColor,
410 linewidth = 0.1,
411 s = 64.0,
412 transform = cartopy.crs.Geodetic(),
413 zorder = 5.0,
414 )
416 # Loop over locations ...
417 for iPnt in range(pntLons.size - 1):
418 # Find the great circle ...
419 circle = great_circle(
420 pntLons[iPnt],
421 pntLats[iPnt],
422 pntLons[iPnt + 1],
423 pntLats[iPnt + 1],
424 debug = debug,
425 eps = eps,
426 maxdist = 12.0 * 1852.0,
427 nIter = nIter,
428 npoint = None,
429 prefix = prefix,
430 ramLimit = ramLimit,
431 )
433 # Draw the great circle ...
434 ax.add_geometries(
435 extract_lines(circle, onlyValid = onlyValid),
436 cartopy.crs.PlateCarree(),
437 edgecolor = skipFillColor if skips[iPnt] or skips[iPnt + 1] else fillColor,
438 facecolor = "none",
439 linewidth = 1.0,
440 zorder = 5.0,
441 )
443 # Check that an extra route was passed ...
444 if route is not None:
445 # Draw the extra route ...
446 ax.add_geometries(
447 extract_lines(route, onlyValid = onlyValid),
448 cartopy.crs.PlateCarree(),
449 edgecolor = routeFillColor,
450 facecolor = "none",
451 linewidth = 1.0,
452 zorder = 5.0,
453 )
455 # Configure axis ...
456 if title is not None:
457 ax.set_title(title)
459 # Configure figure ...
460 fg.tight_layout()
462 # Save figure ...
463 fg.savefig(pngOut)
464 matplotlib.pyplot.close(fg)
466 # Optimise PNG ...
467 optimise_image(
468 pngOut,
469 chunksize = chunksize,
470 debug = debug,
471 exiftoolPath = exiftoolPath,
472 gifsiclePath = gifsiclePath,
473 jpegtranPath = jpegtranPath,
474 optipngPath = optipngPath,
475 pool = None,
476 strip = True,
477 timeout = timeout,
478 )