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