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