Coverage for pyguymer3/geo/create_image_of_points.py: 1%
106 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_image_of_points(
5 pntLons,
6 pntLats,
7 zoom,
8 sess,
9 pngOut,
10 /,
11 *,
12 background = (255, 255, 255),
13 chunksize = 1048576,
14 cookies = None,
15 debug = __debug__,
16 drawGreatCircles = True,
17 drawPointBuffers = False,
18 drawPoints = True,
19 eps = 1.0e-12,
20 exiftoolPath = None,
21 fillColor = (255, 0, 0),
22 gifsiclePath = None,
23 globalExtent = False,
24 globalRatio = 16.0 / 9.0,
25 headers = None,
26 jpegtranPath = None,
27 nAng = 9,
28 nIter = 100,
29 onlyValid = False,
30 optipngPath = None,
31 padDist = 12.0 * 1852.0,
32 prefix = ".",
33 ramLimit = 1073741824,
34 repair = False,
35 route = None,
36 routeFillColor = ( 0, 128, 0),
37 scale = 1,
38 skipFillColor = (255, 165, 0),
39 skips = None,
40 thunderforestKey = None,
41 thunderforestMap = "atlas",
42 timeout = 60.0,
43 tol = 1.0e-10,
44 verify = True,
45):
46 """Save a PNG map of a sequence of points
48 This function accepts a sequence of longitudes and latitudes then saves a
49 PNG map containing all of them drawn together in a big line.
51 Parameters
52 ----------
53 pntLons : numpy.ndarray
54 the sequence of longitudes
55 pntLats : numpy.ndarray
56 the sequence of latitudes
57 zoom : int
58 the OpenStreetMap zoom level
59 sess : requests.Session
60 the session for any requests calls
61 pngOut : str
62 the name of the output PNG
63 background : tuple of int, optional
64 the background colour of the merged tile
65 chunksize : int, optional
66 the size of the chunks of any files which are read in (in bytes)
67 cookies : dict, optional
68 extra cookies for any requests calls
69 debug : bool, optional
70 print debug messages and draw the circle on the axis
71 drawGreatCircles : bool, optional
72 whether to draw the great circles between the points
73 drawPointBuffers : bool, optional
74 whether to draw the buffers around the points
75 drawPoints : bool, optional
76 whether to draw the points
77 eps : float, optional
78 the tolerance of the Vincenty formula iterations
79 exiftoolPath : str, optional
80 the path to the "exiftool" binary (if not provided then Python will
81 attempt to find the binary itself)
82 fillColor : tuple of int, optional
83 the fill colour of the points
84 gifsiclePath : str, optional
85 the path to the "gifsicle" binary (if not provided then Python will
86 attempt to find the binary itself)
87 globalExtent : bool, optional
88 whether to override the calculation of the extent of the points and just
89 make the image of global extent
90 globalRatio : float, optional
91 the ratio to make the image when it is global extent (because the
92 Mercator projection looks silly at the poles)
93 headers : dict, optional
94 extra headers for any requests calls
95 jpegtranPath : str, optional
96 the path to the "jpegtran" binary (if not provided then Python will
97 attempt to find the binary itself)
98 nAng : int, optional
99 the number of angles around the middle location to search over
100 nIter : int, optional
101 the maximum number of iterations (particularly the Vincenty formula)
102 onlyValid : bool, optional
103 only return valid Polygons (checks for validity can take a while, if
104 being called often)
105 optipngPath : str, optional
106 the path to the "optipng" binary (if not provided then Python will
107 attempt to find the binary itself)
108 padDist : float, optional
109 the padding to draw around the points (in metres)
110 prefix : str, optional
111 change the name of the output debugging CSVs
112 ramLimit : int, optional
113 the maximum RAM usage of each "large" array (in bytes)
114 repair : bool, optional
115 attempt to repair invalid Polygons
116 route : shapely.geometry.linestring.LineString, optional
117 an extra line to draw on the map
118 routeFillColor : tuple of int, optional
119 the fill colour of the extra route
120 scale : int, optional
121 the scale of the tiles
122 skipFillColor : tuple of int, optional
123 the fill colour of the skipped points
124 skips : numpy.ndarray, optional
125 an array of booleans as to whether to include/exclude each individual
126 point from calculating the image's field-of-view (this allows the great
127 circles from flights to be drawn but for them to not expand the image to
128 fit in the departing airport); if not provided then all points are used
129 thunderforestKey : string, optional
130 your personal API key for the Thunderforest service (if provided then it
131 is assumed that you want to use the Thunderforest service)
132 thunderforestMap : string, optional
133 the Thunderforest map style (see https://www.thunderforest.com/maps/)
134 timeout : float, optional
135 the timeout for any requests/subprocess calls (in seconds)
136 tol : float, optional
137 the Euclidean distance that defines two points as being the same (in
138 degrees)
139 verify : bool, optional
140 verify the server's certificates for any requests calls
142 Notes
143 -----
144 Copyright 2017 Thomas Guymer [1]_
146 References
147 ----------
148 .. [1] PyGuymer3, https://github.com/Guymer/PyGuymer3
149 """
151 # Import special modules ...
152 try:
153 import numpy
154 except:
155 raise Exception("\"numpy\" is not installed; run \"pip install --user numpy\"") from None
156 try:
157 import PIL
158 import PIL.Image
159 PIL.Image.MAX_IMAGE_PIXELS = 1024 * 1024 * 1024 # [px]
160 import PIL.ImageDraw
161 except:
162 raise Exception("\"PIL\" is not installed; run \"pip install --user Pillow\"") from None
163 try:
164 import shapely
165 import shapely.geometry
166 except:
167 raise Exception("\"shapely\" is not installed; run \"pip install --user Shapely\"") from None
169 # Import sub-functions ...
170 from .buffer import buffer
171 from .extract_lines import extract_lines
172 from .extract_polys import extract_polys
173 from .great_circle import great_circle
174 from .ll2mer import ll2mer
175 from .mer2ll import mer2ll
176 from ..image import image2png
177 from ..openstreetmap import tiles
179 # **************************************************************************
181 # Check inputs ...
182 if skips is None:
183 skips = numpy.zeros(pntLons.size, dtype = bool)
185 # **************************************************************************
187 # Create short-hands ...
188 n = pow(2, zoom)
190 # Create a [Multi]Point from the lists of longitudes and latitudes ...
191 pntsLonLat = []
192 for pntLon, pntLat, skip in zip(pntLons, pntLats, skips, strict = True):
193 if skip:
194 continue
195 pntsLonLat.append(shapely.geometry.point.Point(pntLon, pntLat))
196 pntsLonLat = shapely.geometry.multipoint.MultiPoint(pntsLonLat)
197 if debug:
198 print(f"DEBUG: The points extend from {pntsLonLat.bounds[0]:+.6f}° to {pntsLonLat.bounds[2]:+.6f}° longitude.")
199 print(f"DEBUG: The points extend from {pntsLonLat.bounds[1]:+.6f}° to {pntsLonLat.bounds[3]:+.6f}° latitude.")
201 # Buffer the [Multi]Point ...
202 polysLonLat = buffer(
203 pntsLonLat,
204 padDist,
205 debug = debug,
206 eps = eps,
207 fill = -1.0,
208 fillSpace = "EuclideanSpace",
209 keepInteriors = False,
210 nAng = nAng,
211 nIter = nIter,
212 prefix = prefix,
213 ramLimit = ramLimit,
214 simp = -1.0,
215 tol = tol,
216 )
217 if debug:
218 print(f"DEBUG: The {0.001 * padDist:,.1f} km buffer of the points extends from {polysLonLat.bounds[0]:+.6f}° to {polysLonLat.bounds[2]:+.6f}° longitude.")
219 print(f"DEBUG: The {0.001 * padDist:,.1f} km buffer of the points extends from {polysLonLat.bounds[1]:+.6f}° to {polysLonLat.bounds[3]:+.6f}° latitude.")
221 # Convert [Multi]Polygon to the Mercator projection and create short-hands ...
222 polysMer = ll2mer(
223 polysLonLat,
224 debug = debug,
225 prefix = prefix,
226 tol = tol,
227 )
228 if debug:
229 print(f"DEBUG: The Mercator projection of the {0.001 * padDist:,.1f} km buffer of the points extends from {polysMer.bounds[0]:.6f} to {polysMer.bounds[2]:.6f} in the x-axis of the Mercator projection.")
230 print(f"DEBUG: The Mercator projection of the {0.001 * padDist:,.1f} km buffer of the points extends from {polysMer.bounds[1]:.6f} to {polysMer.bounds[3]:.6f} in the y-axis of the Mercator projection.")
231 if globalExtent:
232 # NOTE: I want the final image to have an aspect ratio of 16:9, which
233 # means that if the width is 1.0 then the height should be 0.5625,
234 # which means that the top should start at 0.21875 and the bottom
235 # should finish at 0.78125.
236 globalBand = 1.0 - (1.0 / globalRatio) # [#]
237 minMerX = 0.0 # [#]
238 minMerY = 0.5 * globalBand # [#]
239 maxMerX = 1.0 # [#]
240 maxMerY = 1.0 - 0.5 * globalBand # [#]
241 midMerX = 0.5 # [#]
242 midMerY = 0.5 # [#]
243 else:
244 minMerX, minMerY, maxMerX, maxMerY = polysMer.bounds # [#], [#], [#], [#]
245 midMerX = 0.5 * (minMerX + maxMerX) # [#]
246 midMerY = 0.5 * (minMerY + maxMerY) # [#]
247 if debug:
248 print(f"DEBUG: The middle of the Mercator projection of the {0.001 * padDist:,.1f} km buffer of the points is at ({midMerX:.6f}, {midMerY:.6f}) in the Mercator projection.")
250 # Convert the middle from Mercator projection back in to longitude and
251 # latitude ...
252 midLon, midLat = mer2ll(
253 shapely.geometry.point.Point(midMerX, midMerY),
254 debug = debug,
255 prefix = prefix,
256 tol = tol,
257 ).coords[0] # [°], [°]
258 if debug:
259 print(f"DEBUG: The middle of the Mercator projection of the {0.001 * padDist:,.1f} km buffer of the points is at ({midLon:+.6f}°, {midLat:+.6f}°).")
261 # Calculate the size of the image ...
262 imgWidth = (maxMerX - minMerX) * float(n * scale * 256) # [px]
263 imgHeight = (maxMerY - minMerY) * float(n * scale * 256) # [px]
264 imgWidth = 2 * round(imgWidth / 2.0) # [px]
265 imgHeight = 2 * round(imgHeight / 2.0) # [px]
266 if debug:
267 print(f"DEBUG: The image is {imgWidth:,d} px × {imgHeight:,d} px.")
269 # Make the image ...
270 img = tiles(
271 midLon,
272 midLat,
273 zoom,
274 imgWidth,
275 imgHeight,
276 sess,
277 background = background,
278 chunksize = chunksize,
279 cookies = cookies,
280 debug = debug,
281 exiftoolPath = exiftoolPath,
282 fill = fillColor,
283 gifsiclePath = gifsiclePath,
284 headers = headers,
285 jpegtranPath = jpegtranPath,
286 optipngPath = optipngPath,
287 radius = None,
288 scale = scale,
289 thunderforestKey = thunderforestKey,
290 thunderforestMap = thunderforestMap,
291 timeout = timeout,
292 verify = verify,
293 )
295 # Create short-hands ...
296 midImgX = imgWidth // 2 # [px]
297 midImgY = imgHeight // 2 # [px]
298 if debug:
299 print(f"DEBUG: The middle of the image image is at ({midImgX:,d} px, {midImgY:,d} px).")
301 # Create short-hand ...
302 draw = PIL.ImageDraw.Draw(img, "RGBA")
304 # Check if the user wants to draw the buffers of the points ...
305 if drawPointBuffers:
306 # Loop over Polygons in the buffer of the points ...
307 for polyMer in extract_polys(polysMer, onlyValid = onlyValid, repair = repair):
308 # Convert LineString to the image projection ...
309 coordsMer = numpy.array(polyMer.exterior.coords) # [#]
310 coordsImgX = float(midImgX) + (coordsMer[:, 0] - midMerX) * float(n * scale * 256) # [px]
311 coordsImgY = float(midImgY) + (coordsMer[:, 1] - midMerY) * float(n * scale * 256) # [px]
313 # Draw the Polygon ...
314 draw.polygon(
315 list(zip(coordsImgX, coordsImgY, strict = True)),
316 fill = fillColor,
317 )
319 # Check if the user wants to draw the points ...
320 if drawPoints:
321 # Loop over points ...
322 for pntLon, pntLat, skip in zip(pntLons, pntLats, skips, strict = True):
323 # Draw the point ...
324 pntMerX, pntMerY = ll2mer(
325 shapely.geometry.point.Point(pntLon, pntLat),
326 debug = debug,
327 prefix = prefix,
328 tol = tol,
329 ).coords[0] # [#], [#]
330 difMerX = pntMerX - midMerX # [#]
331 difMerY = pntMerY - midMerY # [#]
332 difImgX = difMerX * float(n * scale * 256) # [px]
333 difImgY = difMerY * float(n * scale * 256) # [px]
334 pntImgX = float(midImgX) + difImgX # [px]
335 pntImgY = float(midImgY) + difImgY # [px]
336 draw.ellipse(
337 [
338 pntImgX - 10.0,
339 pntImgY - 10.0,
340 pntImgX + 10.0,
341 pntImgY + 10.0,
342 ],
343 fill = skipFillColor if skip else fillColor,
344 )
346 # Check if the user wants to draw the great circles between points ...
347 if drawGreatCircles:
348 # Loop over points ...
349 for iPnt in range(pntLons.size - 1):
350 # Find the great circle from this point to the next ...
351 circleLonLat = great_circle(
352 pntLons[iPnt],
353 pntLats[iPnt],
354 pntLons[iPnt + 1],
355 pntLats[iPnt + 1],
356 debug = debug,
357 eps = eps,
358 maxdist = 12.0 * 1852.0,
359 nIter = nIter,
360 npoint = None,
361 prefix = prefix,
362 ramLimit = ramLimit,
363 )
365 # Loop over LineStrings in the great circle ...
366 for lineLonLat in extract_lines(circleLonLat, onlyValid = onlyValid):
367 # Convert LineString to the Mercator projection ...
368 lineMer = ll2mer(
369 lineLonLat,
370 debug = debug,
371 prefix = prefix,
372 tol = tol,
373 )
375 # Convert LineString to the image projection ...
376 coordsMer = numpy.array(lineMer.coords) # [#]
377 coordsImgX = float(midImgX) + (coordsMer[:, 0] - midMerX) * float(n * scale * 256) # [px]
378 coordsImgY = float(midImgY) + (coordsMer[:, 1] - midMerY) * float(n * scale * 256) # [px]
380 # Draw the line ...
381 draw.line(
382 list(zip(coordsImgX, coordsImgY, strict = True)),
383 fill = skipFillColor if skips[iPnt] or skips[iPnt + 1] else fillColor,
384 width = 4,
385 )
387 # Check that an extra route was passed ...
388 if route is not None:
389 # Loop over LineStrings in the extra route ...
390 for lineLonLat in extract_lines(route, onlyValid = onlyValid):
391 # Convert LineString to the Mercator projection ...
392 lineMer = ll2mer(
393 lineLonLat,
394 debug = debug,
395 prefix = prefix,
396 tol = tol,
397 )
399 # Convert LineString to the image projection ...
400 coordsMer = numpy.array(lineMer.coords) # [#]
401 coordsImgX = float(midImgX) + (coordsMer[:, 0] - midMerX) * float(n * scale * 256) # [px]
402 coordsImgY = float(midImgY) + (coordsMer[:, 1] - midMerY) * float(n * scale * 256) # [px]
404 # Draw the line ...
405 draw.line(
406 list(zip(coordsImgX, coordsImgY, strict = True)),
407 fill = routeFillColor,
408 width = 4,
409 )
411 # Save map ...
412 image2png(
413 img,
414 pngOut,
415 chunksize = chunksize,
416 debug = debug,
417 exif = {
418 "Artist" : "OpenStreetMap contributors",
419 "Copyright" : "All Rights Reserved",
420 "ImageDescription" : "https://www.openstreetmap.org",
421 },
422 exiftoolPath = exiftoolPath,
423 gifsiclePath = gifsiclePath,
424 jpegtranPath = jpegtranPath,
425 mode = "RGB",
426 optimise = True,
427 optipngPath = optipngPath,
428 screenHeight = -1,
429 screenWidth = -1,
430 strip = False,
431 timeout = timeout,
432 )