-
Notifications
You must be signed in to change notification settings - Fork 361
/
geoaxes.py
1902 lines (1564 loc) · 76.2 KB
/
geoaxes.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# (C) British Crown Copyright 2011 - 2017, Met Office
#
# This file is part of cartopy.
#
# cartopy is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by the
# Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# cartopy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with cartopy. If not, see <https://www.gnu.org/licenses/>.
"""
This module defines the :class:`GeoAxes` class, for use with matplotlib.
When a matplotlib figure contains a GeoAxes the plotting commands can transform
plot results from source coordinates to the GeoAxes' target projection.
"""
from __future__ import (absolute_import, division, print_function)
import collections
import contextlib
import warnings
import weakref
import matplotlib as mpl
import matplotlib.artist
import matplotlib.axes
from matplotlib.image import imread
import matplotlib.transforms as mtransforms
import matplotlib.patches as mpatches
import matplotlib.path as mpath
import matplotlib.ticker as mticker
import numpy as np
import numpy.ma as ma
import shapely.geometry as sgeom
from cartopy import config
import cartopy.crs as ccrs
import cartopy.feature
import cartopy.img_transform
from cartopy.mpl.clip_path import clip_path
import cartopy.mpl.feature_artist as feature_artist
import cartopy.mpl.patch as cpatch
from cartopy.mpl.slippy_image_artist import SlippyImageArtist
from cartopy.vector_transform import vector_scalar_to_grid
assert mpl.__version__ >= '1.3', ('Cartopy is only supported with '
'Matplotlib 1.3 or greater.')
_PATH_TRANSFORM_CACHE = weakref.WeakKeyDictionary()
"""
A nested mapping from path, source CRS, and target projection to the
resulting transformed paths::
{path: {(source_crs, target_projection): list_of_paths}}
Provides a significant performance boost for contours which, at
matplotlib 1.2.0 called transform_path_non_affine twice unnecessarily.
"""
_BACKG_IMG_CACHE = {}
"""
A dictionary of pre-loaded images for large background images, kept as a
dictionary so that large images are loaded only once.
"""
_USER_BG_IMGS = {}
"""
A dictionary of background images in the directory specified by the
CARTOPY_USER_BACKGROUNDS environment variable.
"""
# XXX call this InterCRSTransform
class InterProjectionTransform(mtransforms.Transform):
"""
Transforms coordinates from the source_projection to
the ``target_projection``.
"""
input_dims = 2
output_dims = 2
is_separable = False
has_inverse = True
def __init__(self, source_projection, target_projection):
"""
Create the transform object from the given projections.
Args:
* source_projection - A :class:`~cartopy.crs.CRS`.
* target_projection - A :class:`~cartopy.crs.CRS`.
"""
# assert target_projection is cartopy.crs.Projection
# assert source_projection is cartopy.crs.CRS
self.source_projection = source_projection
self.target_projection = target_projection
mtransforms.Transform.__init__(self)
def __repr__(self):
return ('< {!s} {!s} -> {!s} >'.format(self.__class__.__name__,
self.source_projection,
self.target_projection))
def transform_non_affine(self, xy):
"""
Transforms from source to target coordinates.
Args:
* xy - An (n,2) array of points in source coordinates.
Returns:
* An (n,2) array of transformed points in target coordinates.
"""
prj = self.target_projection
if isinstance(xy, np.ndarray):
return prj.transform_points(self.source_projection,
xy[:, 0], xy[:, 1])[:, 0:2]
else:
x, y = xy
x, y = prj.transform_point(x, y, self.source_projection)
return x, y
def transform_path_non_affine(self, src_path):
"""
Transforms from source to target coordinates.
Caches results, so subsequent calls with the same *src_path* argument
(and the same source and target projections) are faster.
Args:
* src_path - A matplotlib :class:`~matplotlib.path.Path` object
with vertices in source coordinates.
Returns
* A matplotlib :class:`~matplotlib.path.Path` with vertices
in target coordinates.
"""
mapping = _PATH_TRANSFORM_CACHE.get(src_path)
if mapping is not None:
key = (self.source_projection, self.target_projection)
result = mapping.get(key)
if result is not None:
return result
# Allow the vertices to be quickly transformed, if
# quick_vertices_transform allows it.
new_vertices = self.target_projection.quick_vertices_transform(
src_path.vertices, self.source_projection)
if new_vertices is not None:
if new_vertices is src_path.vertices:
return src_path
else:
return mpath.Path(new_vertices, src_path.codes)
if src_path.vertices.shape == (1, 2):
return mpath.Path(self.transform(src_path.vertices))
transformed_geoms = []
# Check whether this transform has the "force_path_ccw" attribute set.
# This is a cartopy extension to the Transform API to allow finer
# control of Path orientation handling (Path ordering is not important
# in matplotlib, but is in Cartopy).
geoms = cpatch.path_to_geos(src_path,
getattr(self, 'force_path_ccw', False))
for geom in geoms:
proj_geom = self.target_projection.project_geometry(
geom, self.source_projection)
transformed_geoms.append(proj_geom)
if not transformed_geoms:
result = mpath.Path(np.empty([0, 2]))
else:
paths = cpatch.geos_to_path(transformed_geoms)
if not paths:
return mpath.Path(np.empty([0, 2]))
points, codes = list(zip(*[cpatch.path_segments(path,
curves=False,
simplify=False)
for path in paths]))
result = mpath.Path(np.concatenate(points, 0),
np.concatenate(codes))
# store the result in the cache for future performance boosts
key = (self.source_projection, self.target_projection)
if mapping is None:
_PATH_TRANSFORM_CACHE[src_path] = {key: result}
else:
mapping[key] = result
return result
def inverted(self):
"""
Return a matplotlib :class:`~matplotlib.transforms.Transform`
from target to source coordinates.
"""
return InterProjectionTransform(self.target_projection,
self.source_projection)
class GeoAxes(matplotlib.axes.Axes):
"""
A subclass of :class:`matplotlib.axes.Axes` which represents a
map :class:`~cartopy.crs.Projection`.
This class replaces the matplotlib :class:`~matplotlib.axes.Axes` class
when created with the *projection* keyword. For example::
# Set up a standard map for latlon data.
geo_axes = pyplot.axes(projection=cartopy.crs.PlateCarree())
# Set up an OSGB map.
geo_axes = pyplot.subplot(2, 2, 1, projection=cartopy.crs.OSGB())
When a source projection is provided to one of it's plotting methods,
using the *transform* keyword, the standard matplotlib plot result is
transformed from source coordinates to the target projection. For example::
# Plot latlon data on an OSGB map.
pyplot.axes(projection=cartopy.crs.OSGB())
pyplot.contourf(x, y, data, transform=cartopy.crs.PlateCarree())
"""
def __init__(self, *args, **kwargs):
"""
Create a GeoAxes object using standard matplotlib
:class:`~matplotlib.axes.Axes` args and kwargs.
Kwargs:
* map_projection - The target :class:`~cartopy.crs.Projection` of
this Axes object.
All other args and keywords are passed through to
:class:`matplotlib.axes.Axes`.
"""
self.projection = kwargs.pop('map_projection')
"""The :class:`cartopy.crs.Projection` of this GeoAxes."""
self.outline_patch = None
"""The patch that provides the line bordering the projection."""
self.background_patch = None
"""The patch that provides the filled background of the projection."""
super(GeoAxes, self).__init__(*args, **kwargs)
self._gridliners = []
self.img_factories = []
self._done_img_factory = False
def add_image(self, factory, *args, **kwargs):
"""
Adds an image "factory" to the Axes.
Any image "factory" added, will be asked to retrieve an image
with associated metadata for a given bounding box at draw time.
The advantage of this approach is that the limits of the map
do not need to be known when adding the image factory, but can
be deferred until everything which can effect the limits has been
added.
Currently an image "factory" is just an object with
a ``image_for_domain`` method. Examples of image factories
are :class:`cartopy.io.img_nest.NestedImageCollection` and
:class:`cartopy.io.image_tiles.GoogleTiles`.
"""
if hasattr(factory, 'image_for_domain'):
# XXX TODO: Needs deprecating.
self.img_factories.append([factory, args, kwargs])
else:
# Args and kwargs not allowed.
assert not bool(args) and not bool(kwargs)
image = factory
try:
super(GeoAxes, self).add_image(image)
except AttributeError:
# If add_image method doesn't exist (only available from
# v1.4 onwards) we implement it ourselves.
self._set_artist_props(image)
self.images.append(image)
image._remove_method = lambda h: self.images.remove(h)
return image
@contextlib.contextmanager
def hold_limits(self, hold=True):
"""
Keep track of the original view and data limits for the life of this
context manager, optionally reverting any changes back to the original
values after the manager exits.
Parameters
----------
hold : bool (default True)
Whether to revert the data and view limits after the context
manager exits.
"""
data_lim = self.dataLim.frozen().get_points()
view_lim = self.viewLim.frozen().get_points()
other = (self.ignore_existing_data_limits,
self._autoscaleXon, self._autoscaleYon)
try:
yield
finally:
if hold:
self.dataLim.set_points(data_lim)
self.viewLim.set_points(view_lim)
(self.ignore_existing_data_limits,
self._autoscaleXon, self._autoscaleYon) = other
@matplotlib.artist.allow_rasterization
def draw(self, renderer=None, inframe=False):
"""
Extends the standard behaviour of :func:`matplotlib.axes.Axes.draw`.
Draws grid lines and image factory results before invoking standard
matplotlib drawing. A global range is used if no limits have yet
been set.
"""
# If data has been added (i.e. autoscale hasn't been turned off)
# then we should autoscale the view.
if self.get_autoscale_on() and self.ignore_existing_data_limits:
self.autoscale_view()
if self.outline_patch.reclip or self.background_patch.reclip:
clipped_path = clip_path(self.outline_patch.orig_path,
self.viewLim)
self.outline_patch._path = clipped_path
self.background_patch._path = clipped_path
for gl in self._gridliners:
gl._draw_gridliner(background_patch=self.background_patch)
self._gridliners = []
# XXX This interface needs a tidy up:
# image drawing on pan/zoom;
# caching the resulting image;
# buffering the result by 10%...;
if not self._done_img_factory:
for factory, args, kwargs in self.img_factories:
img, extent, origin = factory.image_for_domain(
self._get_extent_geom(factory.crs), args[0])
self.imshow(img, extent=extent, origin=origin,
transform=factory.crs, *args[1:], **kwargs)
self._done_img_factory = True
return matplotlib.axes.Axes.draw(self, renderer=renderer,
inframe=inframe)
def __str__(self):
return '< GeoAxes: %s >' % self.projection
def cla(self):
"""Clears the current axes and adds boundary lines."""
result = matplotlib.axes.Axes.cla(self)
self.xaxis.set_visible(False)
self.yaxis.set_visible(False)
# Enable tight autoscaling.
self._tight = True
self.set_aspect('equal')
with self.hold_limits():
self._boundary()
# XXX consider a margin - but only when the map is not global...
# self._xmargin = 0.15
# self._ymargin = 0.15
self.dataLim.intervalx = self.projection.x_limits
self.dataLim.intervaly = self.projection.y_limits
return result
def format_coord(self, x, y):
"""Return a string formatted for the matplotlib GUI status bar."""
lon, lat = ccrs.Geodetic().transform_point(x, y, self.projection)
ns = 'N' if lat >= 0.0 else 'S'
ew = 'E' if lon >= 0.0 else 'W'
return u'%.4g, %.4g (%f\u00b0%s, %f\u00b0%s)' % (x, y, abs(lat),
ns, abs(lon), ew)
def coastlines(self, resolution='110m', color='black', **kwargs):
"""
Adds coastal **outlines** to the current axes from the Natural Earth
"coastline" shapefile collection.
Kwargs:
* resolution - a named resolution to use from the Natural Earth
dataset. Currently can be one of "110m", "50m", and
"10m".
.. note::
Currently no clipping is done on the coastlines before adding
them to the axes. This means, if very high resolution coastlines
are being used, performance is likely to be severely effected.
This should be resolved transparently by v0.5.
"""
kwargs['edgecolor'] = color
kwargs['facecolor'] = 'none'
feature = cartopy.feature.NaturalEarthFeature('physical', 'coastline',
resolution, **kwargs)
return self.add_feature(feature)
def tissot(self, rad_km=500, lons=None, lats=None, n_samples=80, **kwargs):
"""
Adds Tissot's indicatrices to the axes.
Kwargs:
* rad_km - The radius in km of the the circles to be drawn.
* lons - A numpy.ndarray, list or tuple of longitude values that
locate the centre of each circle. Specifying more than one
dimension allows individual points to be drawn whereas a
1D array produces a grid of points.
* lats - A numpy.ndarray, list or tuple of latitude values that
that locate the centre of each circle. See lons.
* n_samples - Integer number of points sampled around the
circumference of each circle.
**kwargs are passed through to `class:ShapelyFeature`.
"""
from cartopy import geodesic
geod = geodesic.Geodesic()
geoms = []
if lons is None:
lons = np.linspace(-180, 180, 6, endpoint=False)
else:
lons = np.asarray(lons)
if lats is None:
lats = np.linspace(-80, 80, 6)
else:
lats = np.asarray(lats)
if lons.ndim == 1 or lats.ndim == 1:
lons, lats = np.meshgrid(lons, lats)
lons, lats = lons.flatten(), lats.flatten()
if lons.shape != lats.shape:
raise ValueError('lons and lats must have the same shape.')
for i in range(len(lons)):
circle = geod.circle(lons[i], lats[i], rad_km*1e3,
n_samples=n_samples)
geoms.append(sgeom.Polygon(circle))
feature = cartopy.feature.ShapelyFeature(geoms, ccrs.Geodetic(),
**kwargs)
return self.add_feature(feature)
def natural_earth_shp(self, name='land', resolution='110m',
category='physical', **kwargs):
"""
Adds the geometries from the specified Natural Earth shapefile to the
Axes as a :class:`~matplotlib.collections.PathCollection`.
``**kwargs`` are passed through to the
:class:`~matplotlib.collections.PathCollection` constructor.
Returns the created :class:`~matplotlib.collections.PathCollection`.
.. note::
Currently no clipping is done on the geometries before adding them
to the axes. This means, if very high resolution geometries are
being used, performance is likely to be severely effected. This
should be resolved transparently by v0.5.
"""
warnings.warn('This method has been deprecated.'
' Please use `add_feature` instead.')
kwargs.setdefault('edgecolor', 'face')
kwargs.setdefault('facecolor', cartopy.feature.COLORS['land'])
feature = cartopy.feature.NaturalEarthFeature(category, name,
resolution, **kwargs)
return self.add_feature(feature)
def add_feature(self, feature, **kwargs):
"""
Adds the given :class:`~cartopy.feature.Feature` instance to the axes.
Args:
* feature:
An instance of :class:`~cartopy.feature.Feature`.
Kwargs:
Keyword arguments to be used when drawing the feature. This allows
standard matplotlib control over aspects such as 'facecolor',
'alpha', etc.
Returns:
* A :class:`cartopy.mpl.feature_artist.FeatureArtist`
instance responsible for drawing the feature.
"""
# Instantiate an artist to draw the feature and add it to the axes.
artist = feature_artist.FeatureArtist(feature, **kwargs)
return self.add_artist(artist)
def add_geometries(self, geoms, crs, **kwargs):
"""
Add the given shapely geometries (in the given crs) to the axes.
Args:
* geoms:
A collection of shapely geometries.
* crs:
The cartopy CRS in which the provided geometries are defined.
Kwargs:
Keyword arguments to be used when drawing this feature.
Returns:
A :class:`cartopy.mpl.feature_artist.FeatureArtist`
instance responsible for drawing the geometries.
"""
feature = cartopy.feature.ShapelyFeature(geoms, crs, **kwargs)
return self.add_feature(feature)
def get_extent(self, crs=None):
"""
Get the extent (x0, x1, y0, y1) of the map in the given coordinate
system.
If no crs is given, the returned extents' coordinate system will be
the CRS of this Axes.
"""
p = self._get_extent_geom(crs)
r = p.bounds
x1, y1, x2, y2 = r
return x1, x2, y1, y2
def _get_extent_geom(self, crs=None):
# Perform the calculations for get_extent(), which just repackages it.
with self.hold_limits():
if self.get_autoscale_on():
self.autoscale_view()
[x1, y1], [x2, y2] = self.viewLim.get_points()
domain_in_src_proj = sgeom.Polygon([[x1, y1], [x2, y1],
[x2, y2], [x1, y2],
[x1, y1]])
# Determine target projection based on requested CRS.
if crs is None:
proj = self.projection
elif isinstance(crs, ccrs.Projection):
proj = crs
else:
# Attempt to select suitable projection for
# non-projection CRS.
if isinstance(crs, ccrs.RotatedGeodetic):
proj = ccrs.RotatedPole(crs.proj4_params['lon_0'] - 180,
crs.proj4_params['o_lat_p'])
warnings.warn('Approximating coordinate system {!r} with a '
'RotatedPole projection.'.format(crs))
elif hasattr(crs, 'is_geodetic') and crs.is_geodetic():
proj = ccrs.PlateCarree(crs.globe)
warnings.warn('Approximating coordinate system {!r} with the '
'PlateCarree projection.'.format(crs))
else:
raise ValueError('Cannot determine extent in'
' coordinate system {!r}'.format(crs))
# Calculate intersection with boundary and project if necesary.
boundary_poly = sgeom.Polygon(self.projection.boundary)
if proj != self.projection:
# Erode boundary by threshold to avoid transform issues.
# This is a workaround for numerical issues at the boundary.
eroded_boundary = boundary_poly.buffer(-self.projection.threshold)
geom_in_src_proj = eroded_boundary.intersection(
domain_in_src_proj)
geom_in_crs = proj.project_geometry(geom_in_src_proj,
self.projection)
else:
geom_in_crs = boundary_poly.intersection(domain_in_src_proj)
return geom_in_crs
def set_extent(self, extents, crs=None):
"""
Set the extent (x0, x1, y0, y1) of the map in the given
coordinate system.
If no crs is given, the extents' coordinate system will be assumed
to be the Geodetic version of this axes' projection.
"""
# TODO: Implement the same semantics as plt.xlim and
# plt.ylim - allowing users to set None for a minimum and/or
# maximum value
x1, x2, y1, y2 = extents
domain_in_crs = sgeom.polygon.LineString([[x1, y1], [x2, y1],
[x2, y2], [x1, y2],
[x1, y1]])
projected = None
# Sometimes numerical issues cause the projected vertices of the
# requested extents to appear outside the projection domain.
# This results in an empty geometry, which has an empty `bounds`
# tuple, which causes an unpack error.
# This workaround avoids using the projection when the requested
# extents are obviously the same as the projection domain.
try_workaround = ((crs is None and
isinstance(self.projection, ccrs.PlateCarree)) or
crs == self.projection)
if try_workaround:
boundary = self.projection.boundary
if boundary.equals(domain_in_crs):
projected = boundary
if projected is None:
projected = self.projection.project_geometry(domain_in_crs, crs)
try:
# This might fail with an unhelpful error message ('need more
# than 0 values to unpack') if the specified extents fall outside
# the projection extents, so try and give a better error message.
x1, y1, x2, y2 = projected.bounds
except ValueError:
msg = ('Failed to determine the required bounds in projection '
'coordinates. Check that the values provided are within '
'the valid range (x_limits=[{xlim[0]}, {xlim[1]}], '
'y_limits=[{ylim[0]}, {ylim[1]}]).')
raise ValueError(msg.format(xlim=self.projection.x_limits,
ylim=self.projection.y_limits))
self.set_xlim([x1, x2])
self.set_ylim([y1, y2])
def set_global(self):
"""
Set the extent of the Axes to the limits of the projection.
.. note::
In some cases where the projection has a limited sensible range
the ``set_global`` method does not actually make the whole globe
visible. Instead, the most appropriate extents will be used (e.g.
Ordnance Survey UK will set the extents to be around the British
Isles.
"""
self.set_xlim(self.projection.x_limits)
self.set_ylim(self.projection.y_limits)
def autoscale_view(self, tight=None, scalex=True, scaley=True):
matplotlib.axes.Axes.autoscale_view(self, tight=tight,
scalex=scalex, scaley=scaley)
# Limit the resulting bounds to valid area.
if scalex and self._autoscaleXon:
bounds = self.get_xbound()
self.set_xbound(max(bounds[0], self.projection.x_limits[0]),
min(bounds[1], self.projection.x_limits[1]))
if scaley and self._autoscaleYon:
bounds = self.get_ybound()
self.set_ybound(max(bounds[0], self.projection.y_limits[0]),
min(bounds[1], self.projection.y_limits[1]))
autoscale_view.__doc__ = matplotlib.axes.Axes.autoscale_view.__doc__
def set_xticks(self, ticks, minor=False, crs=None):
"""
Set the x ticks.
Args:
* ticks - list of floats denoting the desired position of x ticks.
Kwargs:
* minor - boolean flag indicating whether the ticks should be minor
ticks i.e. small and unlabelled (default is False).
* crs - An instance of :class:`~cartopy.crs.CRS` indicating the
coordinate system of the provided tick values. If no
coordinate system is specified then the values are assumed
to be in the coordinate system of the projection.
Only transformations from one rectangular coordinate system
to another rectangular coordinate system are supported.
.. note::
This interface is subject to change whilst functionality is added
to support other map projections.
"""
# Project ticks if crs differs from axes' projection
if crs is not None and crs != self.projection:
if not isinstance(crs, (ccrs._RectangularProjection,
ccrs.Mercator)) or \
not isinstance(self.projection,
(ccrs._RectangularProjection,
ccrs.Mercator)):
raise RuntimeError('Cannot handle non-rectangular coordinate '
'systems.')
proj_xyz = self.projection.transform_points(crs,
np.asarray(ticks),
np.zeros(len(ticks)))
xticks = proj_xyz[..., 0]
else:
xticks = ticks
# Switch on drawing of x axis
self.xaxis.set_visible(True)
return super(GeoAxes, self).set_xticks(xticks, minor)
def set_yticks(self, ticks, minor=False, crs=None):
"""
Set the y ticks.
Args:
* ticks - list of floats denoting the desired position of y ticks.
Kwargs:
* minor - boolean flag indicating whether the ticks should be minor
ticks i.e. small and unlabelled (default is False).
* crs - An instance of :class:`~cartopy.crs.CRS` indicating the
coordinate system of the provided tick values. If no
coordinate system is specified then the values are assumed
to be in the coordinate system of the projection.
Only transformations from one rectangular coordinate system
to another rectangular coordinate system are supported.
.. note::
This interface is subject to change whilst functionality is added
to support other map projections.
"""
# Project ticks if crs differs from axes' projection
if crs is not None and crs != self.projection:
if not isinstance(crs, (ccrs._RectangularProjection,
ccrs.Mercator)) or \
not isinstance(self.projection,
(ccrs._RectangularProjection,
ccrs.Mercator)):
raise RuntimeError('Cannot handle non-rectangular coordinate '
'systems.')
proj_xyz = self.projection.transform_points(crs,
np.zeros(len(ticks)),
np.asarray(ticks))
yticks = proj_xyz[..., 1]
else:
yticks = ticks
# Switch on drawing of y axis
self.yaxis.set_visible(True)
return super(GeoAxes, self).set_yticks(yticks, minor)
def stock_img(self, name='ne_shaded'):
"""
Add a standard image to the map.
Currently, the only (and default) option is a downsampled version of
the Natural Earth shaded relief raster.
"""
if name == 'ne_shaded':
import os
source_proj = ccrs.PlateCarree()
fname = os.path.join(config["repo_data_dir"],
'raster', 'natural_earth',
'50-natural-earth-1-downsampled.png')
return self.imshow(imread(fname), origin='upper',
transform=source_proj,
extent=[-180, 180, -90, 90])
else:
raise ValueError('Unknown stock image %r.' % name)
def background_img(self, name='ne_shaded', resolution='low', extent=None,
cache=False):
"""
Adds a background image to the map, from a selection of pre-prepared
images held in a directory specified by the CARTOPY_USER_BACKGROUNDS
environment variable. That directory is checked with
func:`self.read_user_background_images` and needs to contain a JSON
file which defines for the image metadata.
Kwargs:
* name - the name of the image to read according to the contents
of the JSON file. A typical file might have, for instance:
'ne_shaded' : Natural Earth Shaded Relief
'ne_grey' : Natural Earth Grey Earth
* resolution - the resolution of the image to read, according to
the contents of the JSON file. A typical file might
have the following for each name of the image:
'low', 'med', 'high', 'vhigh', 'full'.
* extent - using a high resolution background image, zoomed into
a small area, will take a very long time to render as
the image is prepared globally, even though only a small
area is used. Adding the extent will only render a
particular geographic region. Specified as
[longitude start, longitude end,
latitude start, latitude end].
e.g. [-11, 3, 48, 60] for the UK
or [167.0, 193.0, 47.0, 68.0] to cross the date line.
* cache - logical flag as to whether or not to cache the loaded
images into memory. The images are stored before the
extent is used.
"""
# read in the user's background image directory:
if len(_USER_BG_IMGS) == 0:
self.read_user_background_images()
import os
bgdir = os.getenv('CARTOPY_USER_BACKGROUNDS')
if bgdir is None:
bgdir = os.path.join(config["repo_data_dir"],
'raster', 'natural_earth')
# now get the filename we want to use:
try:
fname = _USER_BG_IMGS[name][resolution]
except KeyError:
msg = ('Image "{}" and resolution "{}" are not present in '
'the user background image metadata in directory "{}"')
raise ValueError(msg.format(name, resolution, bgdir))
# Now obtain the image data from file or cache:
fpath = os.path.join(bgdir, fname)
if cache:
if fname in _BACKG_IMG_CACHE:
img = _BACKG_IMG_CACHE[fname]
else:
img = imread(fpath)
_BACKG_IMG_CACHE[fname] = img
else:
img = imread(fpath)
if len(img.shape) == 2:
# greyscale images are only 2-dimensional, so need replicating
# to 3 colour channels:
img = np.repeat(img[:, :, np.newaxis], 3, axis=2)
# now get the projection from the metadata:
if _USER_BG_IMGS[name]['__projection__'] == 'PlateCarree':
# currently only PlateCarree is defined:
source_proj = ccrs.PlateCarree()
else:
raise NotImplementedError('Background image projection undefined')
if extent is None:
# not specifying an extent, so return all of it:
return self.imshow(img, origin='upper',
transform=source_proj,
extent=[-180, 180, -90, 90])
else:
# return only a subset of the image:
# set up coordinate arrays:
d_lat = 180.0 / img.shape[0]
d_lon = 360.0 / img.shape[1]
# latitude starts at 90N for this image:
lat_pts = (np.arange(img.shape[0]) * -d_lat - (d_lat / 2.0)) + 90.0
lon_pts = (np.arange(img.shape[1]) * d_lon + (d_lon / 2.0)) - 180.0
# which points are in range:
lat_in_range = np.logical_and(lat_pts >= extent[2],
lat_pts <= extent[3])
if extent[0] < 180 and extent[1] > 180:
# we have a region crossing the dateline
# this is the westerly side of the input image:
lon_in_range1 = np.logical_and(lon_pts >= extent[0],
lon_pts <= 180.0)
img_subset1 = img[lat_in_range, :, :][:, lon_in_range1, :]
# and the eastward half:
lon_in_range2 = lon_pts + 360. <= extent[1]
img_subset2 = img[lat_in_range, :, :][:, lon_in_range2, :]
# now join them up:
img_subset = np.concatenate((img_subset1, img_subset2), axis=1)
# now define the extent for output that matches those points:
ret_extent = [lon_pts[lon_in_range1][0] - d_lon / 2.0,
lon_pts[lon_in_range2][-1] + d_lon / 2.0 + 360,
lat_pts[lat_in_range][-1] - d_lat / 2.0,
lat_pts[lat_in_range][0] + d_lat / 2.0]
else:
# not crossing the dateline, so just find the region:
lon_in_range = np.logical_and(lon_pts >= extent[0],
lon_pts <= extent[1])
img_subset = img[lat_in_range, :, :][:, lon_in_range, :]
# now define the extent for output that matches those points:
ret_extent = [lon_pts[lon_in_range][0] - d_lon / 2.0,
lon_pts[lon_in_range][-1] + d_lon / 2.0,
lat_pts[lat_in_range][-1] - d_lat / 2.0,
lat_pts[lat_in_range][0] + d_lat / 2.0]
return self.imshow(img_subset, origin='upper',
transform=source_proj,
extent=ret_extent)
def read_user_background_images(self, verify=True):
"""
Reads the metadata in the specified CARTOPY_USER_BACKGROUNDS
environment variable to populate the dictionaries for background_img.
If CARTOPY_USER_BACKGROUNDS is not set then by default the image in
lib/cartopy/data/raster/natural_earth/ will be made available.
The metadata should be a standard JSON file which specifies a two
level dictionary. The first level is the image type.
For each image type there must be the fields:
__comment__, __source__ and __projection__
and then an element giving the filename for each resolution.
An example JSON file can be found at:
lib/cartopy/data/raster/natural_earth/images.json
"""
import os
import json
bgdir = os.getenv('CARTOPY_USER_BACKGROUNDS')
if bgdir is None:
bgdir = os.path.join(config["repo_data_dir"],
'raster', 'natural_earth')
json_file = os.path.join(bgdir, 'images.json')
with open(json_file, 'r') as js_obj:
dict_in = json.load(js_obj)
for img_type in dict_in:
_USER_BG_IMGS[img_type] = dict_in[img_type]
if verify:
required_info = ['__comment__', '__source__', '__projection__']
for img_type in _USER_BG_IMGS:
if img_type == '__comment__':
# the top level comment doesn't need verifying:
pass
else:
# check that this image type has the required info:
for required in required_info:
if required not in _USER_BG_IMGS[img_type]:
msg = ('User background metadata file "{}", '
'image type "{}", does not specify '
'metadata item "{}"')
raise ValueError(msg.format(json_file, img_type,
required))
for resln in _USER_BG_IMGS[img_type]:
# the required_info items are not resolutions:
if resln not in required_info:
img_it_r = _USER_BG_IMGS[img_type][resln]
test_file = os.path.join(bgdir, img_it_r)
if not os.path.isfile(test_file):
msg = 'File "{}" not found'
raise ValueError(msg.format(test_file))
def add_raster(self, raster_source, **slippy_image_kwargs):
"""
Add the given raster source to the GeoAxes.
Parameters
----------
raster_source : :class:`cartopy.io.RasterSource` like instance
``raster_source`` may be any object which implements the
RasterSource interface, including instances of objects such as
:class:`~cartopy.io.ogc_clients.WMSRasterSource` and
:class:`~cartopy.io.ogc_clients.WMTSRasterSource`. Note that image
retrievals are done at draw time, not at creation time.