Coverage for pyguymer3/geo/_add_topDown_axis.py: 1%
83 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 _add_topDown_axis(
5 fg,
6 lon,
7 lat,
8 /,
9 *,
10 add_coastlines = True,
11 add_gridlines = True,
12 coastlines_edgecolor = "black",
13 coastlines_facecolor = "none",
14 coastlines_levels = None,
15 coastlines_linestyle = "solid",
16 coastlines_linewidth = 0.5,
17 coastlines_resolution = "i",
18 coastlines_zorder = 1.5,
19 configureAgain = False,
20 debug = __debug__,
21 dist = 1.0e99,
22 eps = 1.0e-12,
23 gridlines_int = None,
24 gridlines_linecolor = "black",
25 gridlines_linestyle = ":",
26 gridlines_linewidth = 0.5,
27 gridlines_zorder = 2.0,
28 gs = None,
29 index = None,
30 ncols = None,
31 nIter = 100,
32 nrows = None,
33 onlyValid = False,
34 prefix = ".",
35 ramLimit = 1073741824,
36 repair = False,
37 satellite_height = False,
38 tol = 1.0e-10,
39):
40 """Add an Orthographic axis centred above a point with optionally a
41 field-of-view based on a circle around the point on the surface of the Earth
43 Parameters
44 ----------
45 fg : matplotlib.figure.Figure
46 the figure to add the axis to
47 lon : float
48 the longitude of the point (in degrees)
49 lat : float
50 the latitude of the point (in degrees)
51 add_coastlines : bool, optional
52 add coastline boundaries
53 add_gridlines : bool, optional
54 add gridlines of longitude and latitude
55 coastlines_edgecolor : str, optional
56 the colour of the edges of the coastline Polygons
57 coastlines_facecolor : str, optional
58 the colour of the faces of the coastline Polygons
59 coastlines_levels : list of int, optional
60 the levels of the coastline boundaries (if None then default to
61 ``[1, 6]``)
62 coastlines_linestyle : str, optional
63 the linestyle to draw the coastline boundaries with
64 coastlines_linewidth : float, optional
65 the linewidth to draw the coastline boundaries with
66 coastlines_resolution : str, optional
67 the resolution of the coastline boundaries
68 coastlines_zorder : float, optional
69 the zorder to draw the coastline boundaries with (the default value has
70 been chosen to match the value that it ends up being if the coastline
71 boundaries are not drawn with the zorder keyword specified -- obtained
72 by manual inspection on 5/Dec/2023)
73 configureAgain : bool, optional
74 configure the axis a second time (this is a hack to make narrow
75 field-of-view top-down axes work correctly with OpenStreetMap tiles)
76 debug : bool, optional
77 print debug messages and draw the circle on the axis
78 dist : float, optional
79 the radius of the circle around the point, if larger than the Earth then
80 make the axis of global extent (in metres)
81 eps : float, optional
82 the tolerance of the Vincenty formula iterations
83 gridlines_int : int, optional
84 the interval between gridlines, best results if ``90 % gridlines_int == 0``;
85 if the axis is of global extent then the default will be 45° else it
86 will be 1° (in degrees)
87 gridlines_linecolor : str, optional
88 the colour of the gridlines
89 gridlines_linestyle : str, optional
90 the style of the gridlines
91 gridlines_linewidth : float, optional
92 the width of the gridlines
93 gridlines_zorder : float, optional
94 the zorder to draw the gridlines with (the default value has been chosen
95 to match the value that it ends up being if the gridlines are not drawn
96 with the zorder keyword specified -- obtained by manual inspection on
97 5/Dec/2023)
98 gs : matplotlib.gridspec.SubplotSpec, optional
99 the subset of a gridspec to locate the axis
100 index : int or tuple of int, optional
101 the index of the axis in the array of axes
102 ncols : int, optional
103 the number of columns in the array of axes
104 nIter : int, optional
105 the maximum number of iterations (particularly the Vincenty formula)
106 nrows : int, optional
107 the number of rows in the array of axes
108 onlyValid : bool, optional
109 only return valid Polygons (checks for validity can take a while, if
110 being called often)
111 prefix : str, optional
112 change the name of the output debugging CSVs
113 ramLimit : int, optional
114 the maximum RAM usage of each "large" array (in bytes)
115 repair : bool, optional
116 attempt to repair invalid Polygons
117 satellite_height : bool, optional
118 if a distance is provided then use a "NearsidePerspective" projection at
119 an altitude which has the same field-of-view as the distance
120 tol : float, optional
121 the Euclidean distance that defines two points as being the same (in
122 degrees)
124 Returns
125 -------
126 ax : cartopy.mpl.geoaxes.GeoAxesSubplot
127 the axis
129 Notes
130 -----
131 There are two arguments relating to the `Global Self-Consistent Hierarchical
132 High-Resolution Geography dataset <https://www.ngdc.noaa.gov/mgg/shorelines/>`_ :
134 * *coastlines_levels*; and
135 * *coastlines_resolution*.
137 There are six levels to choose from:
139 * boundary between land and ocean (1);
140 * boundary between lake and land (2);
141 * boundary between island-in-lake and lake (3);
142 * boundary between pond-in-island and island-in-lake (4);
143 * boundary between Antarctica ice and ocean (5); and
144 * boundary between Antarctica grounding-line and ocean (6).
146 There are five resolutions to choose from:
148 * crude ("c");
149 * low ("l");
150 * intermediate ("i");
151 * high ("h"); and
152 * full ("f").
154 Copyright 2017 Thomas Guymer [1]_
156 References
157 ----------
158 .. [1] PyGuymer3, https://github.com/Guymer/PyGuymer3
159 """
161 # Import standard modules ...
162 import math
163 import os
165 # Import special modules ...
166 try:
167 import cartopy
168 cartopy.config.update(
169 {
170 "cache_dir" : os.path.expanduser("~/.local/share/cartopy_cache"),
171 }
172 )
173 except:
174 raise Exception("\"cartopy\" is not installed; run \"pip install --user Cartopy\"") from None
175 try:
176 import matplotlib
177 matplotlib.rcParams.update(
178 {
179 "backend" : "Agg", # NOTE: See https://matplotlib.org/stable/gallery/user_interfaces/canvasagg.html
180 "figure.dpi" : 300,
181 "figure.figsize" : (9.6, 7.2), # NOTE: See https://github.com/Guymer/misc/blob/main/README.md#matplotlib-figure-sizes
182 "font.size" : 8,
183 }
184 )
185 import matplotlib.pyplot
186 except:
187 raise Exception("\"matplotlib\" is not installed; run \"pip install --user matplotlib\"") from None
188 try:
189 import numpy
190 except:
191 raise Exception("\"numpy\" is not installed; run \"pip install --user numpy\"") from None
192 try:
193 import shapely
194 import shapely.geometry
195 except:
196 raise Exception("\"shapely\" is not installed; run \"pip install --user Shapely\"") from None
198 # Import sub-functions ...
199 from ._add_coastlines import _add_coastlines
200 from ._add_horizontal_gridlines import _add_horizontal_gridlines
201 from ._add_vertical_gridlines import _add_vertical_gridlines
202 from .buffer import buffer
203 from .calc_loc_from_loc_and_bearing_and_dist import calc_loc_from_loc_and_bearing_and_dist
204 from .clean import clean
205 from ..consts import CIRCUMFERENCE_OF_EARTH, RADIUS_OF_EARTH
207 # **************************************************************************
209 # Create short-hand ...
210 globalFieldOfView = bool(dist > 0.25 * CIRCUMFERENCE_OF_EARTH)
212 # Check inputs ...
213 if gridlines_int is None:
214 if globalFieldOfView:
215 gridlines_int = 45 # [°]
216 else:
217 gridlines_int = 1 # [°]
218 if not globalFieldOfView and satellite_height:
219 alt = RADIUS_OF_EARTH / math.cos(dist / RADIUS_OF_EARTH) - RADIUS_OF_EARTH # [m]
221 # Create a Point ...
222 point1 = shapely.geometry.point.Point(lon, lat)
224 # Check where the axis should be created ...
225 # NOTE: See https://scitools.org.uk/cartopy/docs/latest/reference/projections.html
226 if gs is not None:
227 # Check if a NearsidePerspective axis can be used ...
228 if not globalFieldOfView and satellite_height:
229 # Create NearsidePerspective axis ...
230 ax = fg.add_subplot(
231 gs,
232 projection = cartopy.crs.NearsidePerspective(
233 central_longitude = point1.x,
234 central_latitude = point1.y,
235 satellite_height = alt,
236 ),
237 )
238 else:
239 # Create Orthographic axis ...
240 ax = fg.add_subplot(
241 gs,
242 projection = cartopy.crs.Orthographic(
243 central_longitude = point1.x,
244 central_latitude = point1.y,
245 ),
246 )
247 elif nrows is not None and ncols is not None and index is not None:
248 # Check if a NearsidePerspective axis can be used ...
249 if not globalFieldOfView and satellite_height:
250 # Create NearsidePerspective axis ...
251 ax = fg.add_subplot(
252 nrows,
253 ncols,
254 index,
255 projection = cartopy.crs.NearsidePerspective(
256 central_longitude = point1.x,
257 central_latitude = point1.y,
258 satellite_height = alt,
259 ),
260 )
261 else:
262 # Create Orthographic axis ...
263 ax = fg.add_subplot(
264 nrows,
265 ncols,
266 index,
267 projection = cartopy.crs.Orthographic(
268 central_longitude = point1.x,
269 central_latitude = point1.y,
270 ),
271 )
272 else:
273 # Check if a NearsidePerspective axis can be used ...
274 if not globalFieldOfView and satellite_height:
275 # Create NearsidePerspective axis ...
276 ax = fg.add_subplot(
277 projection = cartopy.crs.NearsidePerspective(
278 central_longitude = point1.x,
279 central_latitude = point1.y,
280 satellite_height = alt,
281 ),
282 )
283 else:
284 # Create Orthographic axis ...
285 ax = fg.add_subplot(
286 projection = cartopy.crs.Orthographic(
287 central_longitude = point1.x,
288 central_latitude = point1.y,
289 ),
290 )
292 # Check if the field-of-view is too large ...
293 if globalFieldOfView:
294 # Configure axis ...
295 ax.set_global()
296 elif not satellite_height:
297 # Buffer the Point ...
298 polygon1 = buffer(
299 point1,
300 dist,
301 debug = debug,
302 eps = eps,
303 fill = +1.0,
304 fillSpace = "EuclideanSpace",
305 keepInteriors = False,
306 nAng = 361,
307 nIter = nIter,
308 prefix = prefix,
309 ramLimit = ramLimit,
310 simp = -1.0,
311 tol = tol,
312 )
314 # Calculate Northern extent ...
315 tmpLon, latMax, _ = calc_loc_from_loc_and_bearing_and_dist(
316 lon,
317 lat,
318 0.0,
319 dist,
320 eps = 1.0e-12,
321 nIter = 100,
322 ) # [°], [°]
323 tmpPnt = shapely.geometry.point.Point(tmpLon, latMax)
324 yMax = ax.projection.project_geometry(tmpPnt).y # [?]
326 # Calculate Eastern extent ...
327 lonIter, tmpLat, _ = calc_loc_from_loc_and_bearing_and_dist(
328 lon,
329 lat,
330 90.0,
331 dist,
332 eps = 1.0e-12,
333 nIter = 100,
334 ) # [°], [°]
335 tmpPnt = shapely.geometry.point.Point(lonIter, tmpLat)
336 xMax = ax.projection.project_geometry(tmpPnt).x # [?]
338 # Calculate Southern extent ...
339 tmpLon, latMin, _ = calc_loc_from_loc_and_bearing_and_dist(
340 lon,
341 lat,
342 180.0,
343 dist,
344 eps = 1.0e-12,
345 nIter = 100,
346 ) # [°], [°]
347 tmpPnt = shapely.geometry.point.Point(tmpLon, latMin)
348 yMin = ax.projection.project_geometry(tmpPnt).y # [?]
350 # Calculate Western extent ...
351 lonMin, tmpLat, _ = calc_loc_from_loc_and_bearing_and_dist(
352 lon,
353 lat,
354 270.0,
355 dist,
356 eps = 1.0e-12,
357 nIter = 100,
358 ) # [°], [°]
359 tmpPnt = shapely.geometry.point.Point(lonMin, tmpLat)
360 xMin = ax.projection.project_geometry(tmpPnt).x # [?]
362 # Project the Point ...
363 point2 = ax.projection.project_geometry(point1)
365 # Create a correctly oriented Polygon from scratch that is the Point
366 # buffered in MatPlotLib space with the same fuzziness as Cartopy does
367 # internally ...
368 radius2 = numpy.array(
369 [
370 point2.x - xMin,
371 xMax - point2.x,
372 point2.y - yMin,
373 yMax - point2.y,
374 ],
375 dtype = numpy.float64,
376 ).mean() # [?]
377 polygon2 = point2.buffer(radius2 * 0.99999)
378 polygon2 = clean(
379 polygon2,
380 debug = debug,
381 prefix = prefix,
382 tol = tol,
383 )
385 # Convert the exterior ring of the Polygon to a Path ...
386 path = matplotlib.path.Path(polygon2.exterior.coords)
388 # Configure axis ...
389 # NOTE: The orthographic projection does not have the ability to set
390 # either the altitude or the field-of-view. I manually do this,
391 # which involves setting the boundary and the limits for the
392 # MatPlotLib axis.
393 ax.set_boundary(path)
394 ax.set_xlim(xMin, xMax)
395 ax.set_ylim(yMin, yMax)
397 # Check if the user wants to configure the axis a second time ...
398 if configureAgain:
399 # Configure axis again ...
400 # NOTE: For some reason, "cartopy.io.img_tiles.OSM()" doesn't work
401 # unless the first of the following protected members is also
402 # set. All other interactions with the axis appear to be fine
403 # without it being set though. Annoyingly, once the first
404 # protected member is set, all other interactions with the
405 # axis fail unless the second and third protected members are
406 # also set. If the second and third protected members are set
407 # then the resulting gridlines do not extend all the way to
408 # the edge of the map. Therefore, I have chosen to not set the
409 # second and third protected members and instead I protect
410 # setting the first protect member by checking if the
411 # background is going to be "OSM".
412 ax.projection._boundary = polygon2.exterior # pylint: disable=W0212
413 # ax.projection._cw_boundary = polygon2.exterior.reverse()
414 # ax.projection._ccw_boundary = polygon2.exterior
416 # Check if the user wants to draw the circle ...
417 if debug:
418 # Draw the circle ...
419 ax.add_geometries(
420 [polygon1],
421 cartopy.crs.PlateCarree(),
422 edgecolor = (0.0, 0.0, 1.0, 1.0),
423 facecolor = (0.0, 0.0, 1.0, 0.5),
424 linewidth = 1.0,
425 )
427 # Check if the user wants to add coastline boundaries ...
428 if add_coastlines:
429 # Add coastline boundaries ...
430 _add_coastlines(
431 ax,
432 debug = debug,
433 edgecolor = coastlines_edgecolor,
434 facecolor = coastlines_facecolor,
435 levels = coastlines_levels,
436 linestyle = coastlines_linestyle,
437 linewidth = coastlines_linewidth,
438 onlyValid = onlyValid,
439 repair = repair,
440 resolution = coastlines_resolution,
441 zorder = coastlines_zorder,
442 )
444 # Check if the user wants to add gridlines ...
445 if add_gridlines:
446 # Add gridlines ...
447 _add_horizontal_gridlines(
448 ax,
449 color = gridlines_linecolor,
450 linestyle = gridlines_linestyle,
451 linewidth = gridlines_linewidth,
452 locs = range( -90, +90 + gridlines_int, gridlines_int),
453 ngrid = -1,
454 npoint = 3601,
455 zorder = gridlines_zorder,
456 )
457 _add_vertical_gridlines(
458 ax,
459 color = gridlines_linecolor,
460 linestyle = gridlines_linestyle,
461 linewidth = gridlines_linewidth,
462 locs = range(-180, +180 + gridlines_int, gridlines_int),
463 ngrid = -1,
464 npoint = 1801,
465 zorder = gridlines_zorder,
466 )
468 # Return answer ...
469 return ax