-
Notifications
You must be signed in to change notification settings - Fork 8
/
cfplot.py
10548 lines (8909 loc) · 372 KB
/
cfplot.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
"""
Climate contour/vector plots using cf-python, matplotlib and cartopy.
Andy Heaps NCAS-CMS November 2023
"""
import numpy as np
import subprocess
from scipy import interpolate
import matplotlib
from copy import deepcopy
import os
import sys
import matplotlib.pyplot as plot
from matplotlib.collections import PolyCollection
from distutils.version import StrictVersion
import cartopy
import cartopy.crs as ccrs
import cartopy.util as cartopy_util
import cartopy.feature as cfeature
from scipy.interpolate import griddata
import shapely.geometry as sgeom
import shapely
from matplotlib.collections import PatchCollection
import matplotlib.patches as mpatches
# Check for the minimum cf-python version
cf_version_min = '3.0.0b2'
errstr = '\n\n cf-python > ' + cf_version_min
errstr += '\n needs to be installed to use cf-plot \n\n'
try:
import cf
if StrictVersion(cf.__version__) < StrictVersion(cf_version_min):
raise Warning(errstr)
except ImportError:
raise Warning(errstr)
# Initiate the pvars class
# This is used for storing plotting variables in cfp.plotvars
class pvars(object):
def __init__(self, **kwargs):
'''Initialize a new Pvars instance'''
for attr, value in kwargs.items():
setattr(self, attr, value)
def __str__(self):
'''x.__str__() <==> str(x)'''
a = None
v = None
out = ['%s = %s' % (a, repr(v))]
for a, v in self.__dict__.items():
return '\n'.join(out)
# Check for a display and use the Agg backing store if none is present
# This is for batch mode processing
try:
disp = os.environ["DISPLAY"]
except Exception:
matplotlib.use('Agg')
# Check for user setting of pre_existing_data_dir pointing to central cartopy setup
# This is used in the cfview simple setup process
try:
pre_existing_data_dir = os.environ["pre_existing_data_dir"]
cartopy.config["pre_existing_data_dir"] = pre_existing_data_dir
except:
pass
# Code to check if the ImageMagick display command is available
def which(program):
def is_exe(fpath):
return os.path.exists(fpath) and os.access(fpath, os.X_OK)
def ext_candidates(fpath):
yield fpath
for ext in os.environ.get("PATHEXT", "").split(os.pathsep):
yield fpath + ext
for path in os.environ["PATH"].split(os.pathsep):
exe_file = os.path.join(path, program)
for candidate in ext_candidates(exe_file):
if is_exe(candidate):
return candidate
return None
# Default colour scales
# cscale1 is a differential data scale - blue to red
cscale1 = ['#0a3278', '#0f4ba5', '#1e6ec8', '#3ca0f0', '#50b4fa', '#82d2ff',
'#a0f0ff', '#c8faff', '#e6ffff', '#fffadc', '#ffe878', '#ffc03c',
'#ffa000', '#ff6000', '#ff3200', '#e11400', '#c00000', '#a50000']
# viridis is a continuous data scale - blue, green, yellow
viridis = ['#440154', '#440255', '#440357', '#450558', '#45065a', '#45085b',
'#46095c', '#460b5e', '#460c5f', '#460e61', '#470f62', '#471163',
'#471265', '#471466', '#471567', '#471669', '#47186a', '#48196b',
'#481a6c', '#481c6e', '#481d6f', '#481e70', '#482071', '#482172',
'#482273', '#482374', '#472575', '#472676', '#472777', '#472878',
'#472a79', '#472b7a', '#472c7b', '#462d7c', '#462f7c', '#46307d',
'#46317e', '#45327f', '#45347f', '#453580', '#453681', '#443781',
'#443982', '#433a83', '#433b83', '#433c84', '#423d84', '#423e85',
'#424085', '#414186', '#414286', '#404387', '#404487', '#3f4587',
'#3f4788', '#3e4888', '#3e4989', '#3d4a89', '#3d4b89', '#3d4c89',
'#3c4d8a', '#3c4e8a', '#3b508a', '#3b518a', '#3a528b', '#3a538b',
'#39548b', '#39558b', '#38568b', '#38578c', '#37588c', '#37598c',
'#365a8c', '#365b8c', '#355c8c', '#355d8c', '#345e8d', '#345f8d',
'#33608d', '#33618d', '#32628d', '#32638d', '#31648d', '#31658d',
'#31668d', '#30678d', '#30688d', '#2f698d', '#2f6a8d', '#2e6b8e',
'#2e6c8e', '#2e6d8e', '#2d6e8e', '#2d6f8e', '#2c708e', '#2c718e',
'#2c728e', '#2b738e', '#2b748e', '#2a758e', '#2a768e', '#2a778e',
'#29788e', '#29798e', '#287a8e', '#287a8e', '#287b8e', '#277c8e',
'#277d8e', '#277e8e', '#267f8e', '#26808e', '#26818e', '#25828e',
'#25838d', '#24848d', '#24858d', '#24868d', '#23878d', '#23888d',
'#23898d', '#22898d', '#228a8d', '#228b8d', '#218c8d', '#218d8c',
'#218e8c', '#208f8c', '#20908c', '#20918c', '#1f928c', '#1f938b',
'#1f948b', '#1f958b', '#1f968b', '#1e978a', '#1e988a', '#1e998a',
'#1e998a', '#1e9a89', '#1e9b89', '#1e9c89', '#1e9d88', '#1e9e88',
'#1e9f88', '#1ea087', '#1fa187', '#1fa286', '#1fa386', '#20a485',
'#20a585', '#21a685', '#21a784', '#22a784', '#23a883', '#23a982',
'#24aa82', '#25ab81', '#26ac81', '#27ad80', '#28ae7f', '#29af7f',
'#2ab07e', '#2bb17d', '#2cb17d', '#2eb27c', '#2fb37b', '#30b47a',
'#32b57a', '#33b679', '#35b778', '#36b877', '#38b976', '#39b976',
'#3bba75', '#3dbb74', '#3ebc73', '#40bd72', '#42be71', '#44be70',
'#45bf6f', '#47c06e', '#49c16d', '#4bc26c', '#4dc26b', '#4fc369',
'#51c468', '#53c567', '#55c666', '#57c665', '#59c764', '#5bc862',
'#5ec961', '#60c960', '#62ca5f', '#64cb5d', '#67cc5c', '#69cc5b',
'#6bcd59', '#6dce58', '#70ce56', '#72cf55', '#74d054', '#77d052',
'#79d151', '#7cd24f', '#7ed24e', '#81d34c', '#83d34b', '#86d449',
'#88d547', '#8bd546', '#8dd644', '#90d643', '#92d741', '#95d73f',
'#97d83e', '#9ad83c', '#9dd93a', '#9fd938', '#a2da37', '#a5da35',
'#a7db33', '#aadb32', '#addc30', '#afdc2e', '#b2dd2c', '#b5dd2b',
'#b7dd29', '#bade27', '#bdde26', '#bfdf24', '#c2df22', '#c5df21',
'#c7e01f', '#cae01e', '#cde01d', '#cfe11c', '#d2e11b', '#d4e11a',
'#d7e219', '#dae218', '#dce218', '#dfe318', '#e1e318', '#e4e318',
'#e7e419', '#e9e419', '#ece41a', '#eee51b', '#f1e51c', '#f3e51e',
'#f6e61f', '#f8e621', '#fae622', '#fde724']
# Read in defaults if they exist and overlay
# for contour options of fill, blockfill and lines
global_fill = True
global_lines = True
global_blockfill = False
global_degsym = False
global_viewer = 'display'
defaults_file = os.path.expanduser("~") + '/.cfplot_defaults'
if os.path.exists(defaults_file):
with open(defaults_file) as file:
for line in file:
vals = line.split(' ')
com, val = vals
v = False
if val.strip() == 'True':
v = True
if com == 'blockfill':
global_blockfill = v
if com == 'lines':
global_lines = v
if com == 'fill':
global_fill = v
if com == 'degsym':
global_degsym = v
if com == 'viewer':
global_viewer = val.strip()
# plotvars - global plotting variables
plotvars = pvars(lonmin=-180, lonmax=180, latmin=-90, latmax=90, proj='cyl',
resolution='110m', plot_type=1, boundinglat=0, lon_0=0,
levels=None,
levels_min=None, levels_max=None, levels_step=None,
norm=None, levels_extend='both', xmin=None,
xmax=None, ymin=None, ymax=None, xlog=None, ylog=None,
rows=1, columns=1, file=None, orientation='landscape',
user_mapset=0, user_gset=0, cscale_flag=0, user_levs=0,
user_plot=0, master_plot=None, plot=None, cs=cscale1,
cs_user='cscale1', mymap=None, xticks=None, yticks=None,
xticklabels=None, yticklabels=None, xstep=None, ystep=None,
xlabel=None, ylabel=None, title=None, title_fontsize=15,
axis_label_fontsize=11, text_fontsize=11,
text_fontweight='normal', axis_label_fontweight='normal',
colorbar_fontsize=11, colorbar_fontweight='normal',
title_fontweight='normal', continent_thickness=None,
continent_color=None, continent_linestyle=None,
pos=1, viewer=global_viewer, global_viewer=global_viewer,
tspace_year=None, tspace_month=None, tspace_day=None,
tspace_hour=None, xtick_label_rotation=0,
xtick_label_align='center', ytick_label_rotation=0,
ytick_label_align='right', legend_text_size=11,
legend_text_weight='normal',
cs_uniform=True, master_title=None,
master_title_location=[0.5, 0.95], master_title_fontsize=30,
master_title_fontweight='normal', dpi=None,
plot_xmin=None, plot_xmax=None, plot_ymin=None,
plot_ymax=None, land_color=None, ocean_color=None,
lake_color=None, feature_zorder=99, twinx=False, twiny=False,
rotated_grid_thickness=2.0, rotated_grid_spacing=10,
rotated_deg_spacing=0.75, rotated_continents=True,
rotated_grid=True, rotated_labels=True,
legend_frame=True, legend_frame_edge_color='k',
legend_frame_face_color=None, degsym=global_degsym,
axis_width=None, grid_x_spacing=60, grid_y_spacing=30,
grid_colour='k', grid_linestyle='--', grid_zorder=100,
grid_thickness=1.0, aspect='equal',
graph_xmin=None, graph_xmax=None,
graph_ymin=None, graph_ymax=None,
level_spacing=None, tight=False, gpos_called=False,
titles_con_called=False)
# Check for iPython notebook inline
# and set the viewer to None if found
is_inline = 'inline' in matplotlib.get_backend()
if is_inline:
plotvars.viewer = None
# Check for OSX and if so use matplotlib for for the viewer
# Not all users will have ImageMagick installed / XQuartz running
# Users can still select this with cfp.setvars(viewer='display')
if sys.platform == 'darwin':
plotvars.global_viewer = 'matplotlib'
plotvars.viewer = 'matplotlib'
def con(f=None, x=None, y=None, fill=global_fill, lines=global_lines, line_labels=True,
title=None, colorbar_title=None, colorbar=True,
colorbar_label_skip=None, ptype=0, negative_linestyle='solid',
blockfill=global_blockfill, zero_thick=False, colorbar_shrink=None,
colorbar_orientation=None, colorbar_position=None, xlog=False,
ylog=False, axes=True, xaxis=True, yaxis=True, xticks=None,
xticklabels=None, yticks=None, yticklabels=None, xlabel=None,
ylabel=None, colors='k', swap_axes=False, verbose=None,
linewidths=None, alpha=1.0, colorbar_text_up_down=False,
colorbar_fontsize=None, colorbar_fontweight=None,
colorbar_text_down_up=False, colorbar_drawedges=True,
colorbar_fraction=None, colorbar_thick=None,
colorbar_anchor=None, colorbar_labels=None,
linestyles=None, zorder=1, level_spacing=None,
irregular=None, face_lons=False, face_lats=False, face_connectivity=False,
titles=False, mytest=False, transform_first=None, blockfill_fast=None,
nlevs=False, orca=None, orca_skip=None, grid=False):
"""
| con is the interface to contouring in cf-plot. The minimum use is con(f)
| where f is a 2 dimensional array. If a cf field is passed then an
| appropriate plot will be produced i.e. a longitude-latitude or
| latitude-height plot for example. If a 2d numeric array is passed then
| the optional arrays x and y can be used to describe the x and y data
| point locations.
|
| f - array to contour
| x - x locations of data in f (optional)
| y - y locations of data in f (optional)
| fill=True - colour fill contours
| lines=True - draw contour lines and labels
| line_labels=True - label contour lines
| title=title - title for the plot
| ptype=0 - plot type - not needed for cf fields.
| 0 = no specific plot type,
| 1 = longitude-latitude,
| 2 = latitude - height,
| 3 = longitude - height,
| 4 = latitude - time,
| 5 = longitude - time
| 6 = rotated pole
| negative_linestyle='solid' - set to one of 'solid', 'dashed'
| zero_thick=False - add a thick zero contour line. Set to 3 for example.
| blockfill=False - set to True for a blockfill plot
| colorbar_title=colbar_title - title for the colour bar
| colorbar=True - add a colour bar if a filled contour plot
| colorbar_label_skip=None - skip colour bar labels. Set to 2 to skip
| every other label.
| colorbar_orientation=None - options are 'horizontal' and 'vertical'
| The default for most plots is horizontal but
| for polar stereographic plots this is vertical.
| colorbar_shrink=None - value to shrink the colorbar by. If the colorbar
| exceeds the plot area then values of 1.0, 0.55
| or 0.5m may help it better fit the plot area.
| colorbar_position=None - position of colorbar
| [xmin, ymin, x_extent,y_extent] in normalised
| coordinates. Use when a common colorbar
| is required for a set of plots. A typical set
| of values would be [0.1, 0.05, 0.8, 0.02]
| colorbar_fontsize=None - text size for colorbar labels and title
| colorbar_fontweight=None - font weight for colorbar labels and title
| colorbar_text_up_down=False - if True horizontal colour bar labels alternate
| above (start) and below the colour bar
| colorbar_text_down_up=False - if True horizontal colour bar labels alternate
| below (start) and above the colour bar
| colorbar_drawedges=True - draw internal divisions in the colorbar
| colorbar_fraction=None - space for the colorbar - default = 0.21, in normalised
| coordinates
| colorbar_thick=None - thickness of the colorbar - default = 0.015, in normalised
| coordinates
| colorbar_anchor=None - default=0.5 - anchor point of colorbar within the fraction space.
| 0.0 = close to plot, 1.0 = further away
| colorbar_labels=None - labels to use for colorbar. The default is to use the contour
| levels as labels
| colorbar_text_up_down=False - on a horizontal colorbar alternate the
| labels top and bottom starting in the up position
| colorbar_text_down_up=False - on a horizontal colorbar alternate the
| labels bottom and top starting in the bottom position
| colorbar_drawedges=True - draw internal delimeter lines in the colorbar
| colors='k' - contour line colors - takes one or many values.
| xlog=False - logarithmic x axis
| ylog=False - logarithmic y axis
| axes=True - plot x and y axes
| xaxis=True - plot xaxis
| yaxis=True - plot y axis
| xticks=None - xtick positions
| xticklabels=None - xtick labels
| yticks=None - y tick positions
| yticklabels=None - ytick labels
| xlabel=None - label for x axis
| ylabel=None - label for y axis
| swap_axes=False - swap plotted axes - only valid for X, Y, Z vs T plots
| verbose=None - change to 1 to get a verbose listing of what con
| is doing
| linewidths=None - contour linewidths. Either a single number for all
| lines or array of widths
| linestyles=None - takes 'solid', 'dashed', 'dashdot' or 'dotted'
| alpha=1.0 - transparency setting. The default is no transparency.
| zorder=1 - order of drawing
| level_spacing=None - Default of 'linear' level spacing. Also takes 'log', 'loglike',
| 'outlier' and 'inspect'
| irregular=None - flag for contouring irregular data
| face_lons=None - longitude points for face vertices
| face_lats=None - latitude points for face verticies
| face_connectivity=None - connectivity for face verticies
| titles=False - set to True to have a dimensions title
| transform_first=None - Cartopy should transform the points before calling the contouring algorithm,
| which can have a significant impact on speed (it is much faster to transform
| points than it is to transform patches) If this is unset and the number of points
| in the x direction is > 400 then it is set to True.
| blockfill_fast=None - Use pcolormesh blockfill. This is possibly less reliable that the usual code but is
| faster for higher resolution datasets
| nlevs=False - Let Matplotlib work out the levels for the contour plot
| orca=None - User specifies this is an orca tripolar grid. Internally cf-plot tries to detect this by looking
| for a single discontinuity in the logitude 2D array. If found a fix it make to the longitudes so
| that they are no longer discontinuous.
| orca_skip=None - Only plot every nth grid point in the 2D longitude and latitude arrays. This is useful for when
| plotting his resolution data over the whole globe which would otherwise be very slow to visualize.
| grid=False - Draw a grid on the map using the parameters set by cfp.setvars. Defaults are grid_x_spacing=60,
| grid_y_spacing=30, grid_colour='k', grid_linestyle = '--', grid_thickness=1.0
|
:Returns:
None
"""
# Turn off divide warning in contour routine which is a numpy issue
old_settings = np.seterr(all='ignore')
np.seterr(divide='ignore')
# Set potential user axis labels
user_xlabel = xlabel
user_ylabel = ylabel
# Set blockfill to True if blockfill_fast is not None
if blockfill_fast is not None:
blockfill=True
# Check if the field is a CF ugrid field with cell faces
blockfill_ugrid = False
if isinstance(f, cf.Field) and blockfill:
if f.domain_topologies():
if f.domain_topology('cell:face', default=None) is not None:
face_lons_array = f.aux('X').bounds.array
face_lats_array = f.aux('Y').bounds.array
face_connectivity_array = f.domain_topology('cell:face').array
blockfill_ugrid = True
fill = False
lines = False
irregular = True
else:
errstr = '\n\nError - field does not contain the UGRID face information to plot a blockfill plot\n\n\n'
raise TypeError(errstr)
# Set blockfill_2d if blockfill and x and y are 2D
blockfill_2d = False
if blockfill and not isinstance(f, cf.Field):
if np.ndim(x) == 2 and np.ndim(y) == 2:
blockfill_2d = True
# Call gpos(1) if not already called
if plotvars.rows > 1 or plotvars.columns > 1:
if plotvars.gpos_called is False:
gpos(1)
# Extract required data for contouring
# If a cf-python field
if isinstance(f, cf.Field):
ndims = np.squeeze(f.data).ndim
if ndims > 2:
errstr = "\n\ncfp.con error need a 1 or 2 dimensional field to contour\n"
errstr += "received " + str(np.squeeze(f.data).ndim) + " dimensions\n\n"
errstr += str(f)
raise TypeError(errstr)
# Extract data
if verbose:
print('con - calling cf_data_assign')
# Subset the data if a user map is set
# This is use to speed up the plotting
# myfield is used for calculating the contour levels
# myfield_extended is used to make the contour plot
if plotvars.user_mapset:
if plotvars.proj == 'npstere':
f = f.subspace(Y = cf.wi(plotvars.boundinglat, 90.0))
if plotvars.proj == 'spstere':
f = f.subspace(Y = cf.wi(-90.0, plotvars.boundinglat))
# Extract the data
field, x, y, ptype, colorbar_title, xlabel, ylabel, xpole, ypole =\
cf_data_assign(f, colorbar_title, verbose=verbose)
if user_xlabel is not None:
xlabel = user_xlabel
if user_ylabel is not None:
ylabel = user_ylabel
elif isinstance(f, cf.FieldList):
raise TypeError("\n\ncfp.con - cannot contour a field list\n\n")
else:
if verbose:
print('con - using user assigned data')
field = f # field data passed in as f
if x is None:
x = np.arange(np.shape(field)[1])
if y is None:
y = np.arange(np.shape(field)[0])
check_data(field, x, y)
xlabel = ''
ylabel = ''
# Assign irregular and orca keywords unless already set
if irregular is None:
if np.size(x) == np.size(np.unique(x)):
irregular = False
else:
irregular = True
if np.ndim(x) == 2 and np.ndim(y) == 2:
if orca is None:
orca = orca_check(x)
if orca:
# Apply orca_skip if set
if orca_skip is not None:
print('applying orca_skip value of ', orca_skip)
x = x[::orca_skip, ::orca_skip]
y = y[::orca_skip, ::orca_skip]
field = field[::orca_skip, ::orca_skip]
# orca grids have a discontinuity in the longitude grid
# use the method at https://gist.github.com/pelson/79cf31ef324774c97ae7
# to remove the discontinuity
fixed_x = x.copy()
for i, start in enumerate(np.argmax(np.abs(np.diff(x)) > 180, axis=1)):
fixed_x[i, start+1:] += 360
x = fixed_x
if np.ndim(x) == 2:
irregular = False
# Set contour line styles
matplotlib.rcParams['contour.negative_linestyle'] = negative_linestyle
# Set contour lines off on block plots
if blockfill:
fill = False
field_orig = deepcopy(field)
x_orig = deepcopy(x)
y_orig = deepcopy(y)
# Check number of colours and levels match if user has modified the
# number of colours
if plotvars.cscale_flag == 2:
ncols = np.size(plotvars.cs)
nintervals = np.size(plotvars.levels) - 1
if plotvars.levels_extend == 'min':
nintervals += 1
if plotvars.levels_extend == 'max':
nintervals += 1
if plotvars.levels_extend == 'both':
nintervals += 2
if ncols != nintervals:
errstr = "\n\ncfp.con - blockfill error \n"
errstr += "need to match number of colours and contour intervals\n"
errstr += "Don't forget to take account of the colorbar "
errstr += "extensions\n\n"
raise TypeError(errstr)
# Turn off colorbar if fill is turned off
if not fill and not blockfill and not blockfill_ugrid:
colorbar = False
# Revert to default colour scale if cscale_flag flag is set
if plotvars.cscale_flag == 0:
plotvars.cs = cscale1
# Set the orientation of the colorbar
if plotvars.plot_type == 1:
if plotvars.proj == 'npstere' or plotvars.proj == 'spstere':
if colorbar_orientation is None:
colorbar_orientation = 'vertical'
if colorbar_orientation is None:
colorbar_orientation = 'horizontal'
# Store original map resolution
resolution_orig = plotvars.resolution
# Set size of color bar if not specified
if colorbar_shrink is None:
colorbar_shrink = 1.0
if plotvars.proj == 'npstere' or plotvars.proj == 'spstere':
colorbar_shrink = 0.8
# Set plot type if user specified
if (ptype is not None):
plotvars.plot_type = ptype
# Get contour levels if none are defined
spacing = 'linear'
if plotvars.level_spacing is not None:
spacing = plotvars.level_spacing
if level_spacing is not None:
spacing = level_spacing
if plotvars.levels is None:
if isinstance(f, cf.Field):
field, x, y, ptype, colorbar_title, xlabel, ylabel, xpole, ypole =\
cf_data_assign(f, colorbar_title, verbose=verbose)
clevs, mult, fmult = calculate_levels(field=field,
level_spacing=spacing,
verbose=verbose)
else:
clevs = plotvars.levels
mult = 0
fmult = 1
# Set the colour scale if nothing is defined
includes_zero = False
if plotvars.cscale_flag == 0:
col_zero = 0
for cval in clevs:
if includes_zero is False:
col_zero = col_zero + 1
if cval == 0:
includes_zero = True
if includes_zero:
cs_below = col_zero
cs_above = np.size(clevs) - col_zero + 1
if plotvars.levels_extend == 'max' or plotvars.levels_extend == 'neither':
cs_below = cs_below - 1
if plotvars.levels_extend == 'min' or plotvars.levels_extend == 'neither':
cs_above = cs_above - 1
uniform = True
if plotvars.cs_uniform is False:
uniform = False
cscale('scale1', below=cs_below, above=cs_above, uniform=uniform)
else:
ncols = np.size(clevs)+1
if plotvars.levels_extend == 'min' or plotvars.levels_extend == 'max':
ncols = ncols-1
if plotvars.levels_extend == 'neither':
ncols = ncols-2
cscale('viridis', ncols=ncols)
plotvars.cscale_flag = 0
# User selected colour map but no mods so fit to levels
if plotvars.cscale_flag == 1:
ncols = np.size(clevs)+1
if plotvars.levels_extend == 'min' or plotvars.levels_extend == 'max':
ncols = ncols-1
if plotvars.levels_extend == 'neither':
ncols = ncols-2
cscale(plotvars.cs_user, ncols=ncols)
plotvars.cscale_flag = 1
# Set colorbar labels
# Set a sensible label spacing if the user hasn't already done so
if colorbar_label_skip is None:
if colorbar_orientation == 'horizontal':
nchars = 0
for lev in clevs:
nchars = nchars + len(str(lev))
colorbar_label_skip = int(nchars / 80 + 1)
if plotvars.columns > 1:
colorbar_label_skip = int(nchars * (plotvars.columns) / 80)
else:
colorbar_label_skip = 1
if colorbar_label_skip > 1:
if includes_zero:
# include zero in the colorbar labels
zero_pos = [i for i, mylev in enumerate(clevs) if mylev == 0][0]
cbar_labels = clevs[zero_pos]
i = zero_pos + colorbar_label_skip
while i <= len(clevs) - 1:
cbar_labels = np.append(cbar_labels, clevs[i])
i = i + colorbar_label_skip
i = zero_pos - colorbar_label_skip
if i >= 0:
while i >= 0:
cbar_labels = np.append(clevs[i], cbar_labels)
i = i - colorbar_label_skip
else:
cbar_labels = clevs[0]
i = int(colorbar_label_skip)
while i <= len(clevs) - 1:
cbar_labels = np.append(cbar_labels, clevs[i])
i = i + colorbar_label_skip
else:
cbar_labels = clevs
if colorbar_label_skip is None:
colorbar_label_skip = 1
# Make a list of strings of the colorbar levels for later labelling
clabels = []
for i in cbar_labels:
clabels.append(str(i))
if colorbar_label_skip > 1:
for skip in np.arange(colorbar_label_skip - 1):
clabels.append('')
if colorbar_labels is not None:
cbar_labels = colorbar_labels
else:
cbar_labels = clabels
# Turn off line_labels if the field is all the same
# Matplotlib 3.2.2 throws an error if there are no line labels
if np.nanmin(field) == np.nanmax(field):
line_labels = False
# Add mult to colorbar_title if used
if (colorbar_title is None):
colorbar_title = ''
if (mult != 0):
colorbar_title = colorbar_title + ' *10$^{' + str(mult) + '}$'
# Catch null title
if title is None:
title = ''
if plotvars.title is not None:
title = plotvars.title
# Set plot variables
title_fontsize = plotvars.title_fontsize
text_fontsize = plotvars.text_fontsize
if colorbar_fontsize is None:
colorbar_fontsize = plotvars.colorbar_fontsize
if colorbar_fontweight is None:
colorbar_fontweight = plotvars.colorbar_fontweight
continent_thickness = plotvars.continent_thickness
continent_color = plotvars.continent_color
continent_linestyle = plotvars.continent_linestyle
land_color = plotvars.land_color
ocean_color = plotvars.ocean_color
lake_color = plotvars.lake_color
title_fontweight = plotvars.title_fontweight
if continent_thickness is None:
continent_thickness = 1.5
if continent_color is None:
continent_color = 'k'
if continent_linestyle is None:
continent_linestyle = 'solid'
cb_orient = colorbar_orientation
# Retrieve any user defined axis labels
if xlabel == '' and plotvars.xlabel is not None:
xlabel = plotvars.xlabel
if ylabel == '' and plotvars.ylabel is not None:
ylabel = plotvars.ylabel
if xticks is None and plotvars.xticks is not None:
xticks = plotvars.xticks
if plotvars.xticklabels is not None:
xticklabels = plotvars.xticklabels
else:
xticklabels = list(map(str, xticks))
if yticks is None and plotvars.yticks is not None:
yticks = plotvars.yticks
if plotvars.yticklabels is not None:
yticklabels = plotvars.yticklabels
else:
yticklabels = list(map(str, yticks))
# Calculate a set of dimension titles if requested
if titles:
plotvars.titles_con_called = True
title_dims = generate_titles(f)
if not colorbar:
title_dims = colorbar_title + '\n' + title_dims
# Check if data is well formed
# i.e. dimensions have only recognizable X, Y, Z, T or a subset
well_formed = False
if isinstance(f, cf.Field):
well_formed = check_well_formed(f)
if nlevs is not False:
clevs = nlevs
plotvars.levels_extend = 'neither'
if plotvars.cscale_flag == 0:
if np.min(field) < 0 and np.max(field) > 0:
cscale('scale1', ncols=nlevs)
else:
cscale('viridis', ncols=nlevs)
plotvars.cscale_flag = 0
else:
cscale(plotvars.cs_user, ncols=nlevs)
##################
# Map contour plot
##################
if ptype == 1:
if verbose:
print('con - making a map plot')
# Open a new plot if necessary
if plotvars.user_plot == 0:
gopen(user_plot=0)
# Reset the stored mapping
if plotvars.user_mapset == 0:
plotvars.lonmin = -180
plotvars.lonmax = 180
plotvars.latmin = -90
plotvars.latmax = 90
# Set up mapping
mylonmin = np.nanmin(x)
mylonmax = np.nanmax(x)
mylatmin = np.nanmin(y)
mylatmax = np.nanmax(y)
lonrange = mylonmax - mylonmin
latrange = mylatmax - mylatmin
if lonrange > 360.0:
mylonmax = mylonmin + 360.0
lonrange = 360.0
if (lonrange > 350 and latrange > 160) or plotvars.user_mapset == 1:
set_map()
else:
mapset(lonmin=mylonmin, lonmax=mylonmax, latmin=mylatmin, latmax=mylatmax,
user_mapset=0, resolution=resolution_orig)
set_map()
mymap = plotvars.mymap
user_mapset = plotvars.user_mapset
lonrange = np.nanmax(x) - np.nanmin(x)
if not blockfill_ugrid and not blockfill_2d:
if not irregular:
if lonrange > 350 and np.ndim(y) == 1:
# Add cyclic information if missing.
if lonrange < 360:
# field, x = cartopy_util.add_cyclic_point(field, x)
# Call add_cyclic_point it spacing is regular
x_regular = True
xspacing = x[1] - x[0]
for ix in np.arange(len(x) - 1):
if x[ix+1] - x[ix] != xspacing:
x_regular = False
if x_regular:
field, x = add_cyclic(field, x)
lonrange = np.nanmax(x) - np.nanmin(x)
# cartopy line drawing fix
if x[-1] - x[0] == 360.0:
x[-1] = x[-1] + 0.001
# Shift grid if needed
if plotvars.lonmin < np.nanmin(x):
# Cartopy feature at version 0.20.0
# -360 to 0 creates strange contours
vers = cartopy.__version__.split('.')
val = int(vers[0] +vers[1])
if val < 20:
x = x - 360
if plotvars.lonmin > np.nanmax(x):
x = x + 360
elif not orca:
# Get the irregular data within the map coordinates
# Matplotlib tricontour cannot plot missing data so we need to split
# the missing data into a separate field to deal with this
field_modified = deepcopy(field)
pts_nan = np.where(np.isnan(field_modified))
field_modified[pts_nan] = -1e30
field_irregular, lons_irregular, lats_irregular = irregular_window(field_modified, x, y)
#pts_real = np.where(np.isfinite(field_irregular))
pts_real = np.where(field_irregular > -1e29)
pts_nan = np.where(field_irregular < -1e29)
field_irregular_nan = []
lons_irregular_nan = []
lats_irregular_nan = []
if np.size(pts_nan) > 0:
field_irregular_nan = deepcopy(field_irregular)
lons_irregular_nan = deepcopy(lons_irregular)
lats_irregular_nan = deepcopy(lats_irregular)
field_irregular_nan[:] = 0
field_irregular_nan[pts_nan] = 1
field_irregular_real = deepcopy(field_irregular[pts_real])
lons_irregular_real = deepcopy(lons_irregular[pts_real])
lats_irregular_real = deepcopy(lats_irregular[pts_real])
if not irregular:
# Flip latitudes and field if latitudes are in descending order
if np.ndim(y) == 1:
if y[0] > y[-1]:
y = y[::-1]
field = np.flipud(field)
# Plotting a sub-area of the grid produces stray contour labels
# in polar plots. Subsample the latitudes to remove this problem
if plotvars.proj == 'npstere' and np.ndim(y) == 1:
if not blockfill_ugrid and not blockfill_2d:
if irregular:
pts = np.where(lats_irregular > plotvars.boundinglat - 5)
pts = np.array(pts).flatten()
lons_irregular_real = lons_irregular_real[pts]
lats_irregular_real = lats_irregular_real[pts]
field_irregular_real = field_irregular_real[pts]
else:
myypos = find_pos_in_array(vals=y, val=plotvars.boundinglat)
if myypos != -1:
y = y[myypos:]
field = field[myypos:, :]
if plotvars.proj == 'spstere' and np.ndim(y) == 1:
if not blockfill_ugrid and not blockfill_2d:
if irregular:
pts = np.where(lats_irregular_real < plotvars.boundinglat + 5)
lons_irregular_real = lons_irregular_real[pts]
lats_irregular_real = lats_irregular_real[pts]
field_irregular_real = field_irregular_real[pts]
else:
myypos = find_pos_in_array(vals=y, val=plotvars.boundinglat, above=True)
if myypos != -1:
y = y[0:myypos + 1]
field = field[0:myypos + 1, :]
# Set the longitudes and latitudes
lons, lats = x, y
# Set the plot limits
if lonrange > 350:
gset(
xmin=plotvars.lonmin,
xmax=plotvars.lonmax,
ymin=plotvars.latmin,
ymax=plotvars.latmax,
user_gset=0)
else:
if user_mapset == 1:
gset(xmin=plotvars.lonmin,
xmax=plotvars.lonmax,
ymin=plotvars.latmin,
ymax=plotvars.latmax,
user_gset=0)
else:
gset(xmin=np.nanmin(lons),
xmax=np.nanmax(lons),
ymin=np.nanmin(lats),
ymax=np.nanmax(lats),
user_gset=0)
# Filled contours
if fill:
if verbose:
print('con - adding filled contours')
# Get colour scale for use in contouring
# If colour bar extensions are enabled then the colour map goes
# from 1 to ncols-2. The colours for the colour bar extensions
# are then changed on the colorbar and plot after the plot is made
colmap = cscale_get_map()
cmap = matplotlib.colors.ListedColormap(colmap)
if (plotvars.levels_extend ==
'min' or plotvars.levels_extend == 'both'):
cmap.set_under(plotvars.cs[0])
if (plotvars.levels_extend ==
'max' or plotvars.levels_extend == 'both'):
cmap.set_over(plotvars.cs[-1])
# For fast map contours add transform_first=True to contourf command
# and make lons and lats 2D
if transform_first is None and np.ndim(lons) == 1 and np.ndim(lats) == 1:
if np.size(lons) >= 400:
transform_first = True
# Fast map contours are also needed when clevs is a integer
if type(clevs) == int and plotvars.plot_type == 1 and plotvars.proj == 'cyl':
transform_first=True
if transform_first:
if np.ndim(lons) == 1 and np.ndim(lats) == 1:
lons, lats = np.meshgrid(lons, lats)
# Filled colour contours
if not irregular or orca is True:
plotvars.image = mymap.contourf(lons, lats, field * fmult, clevs,
extend=plotvars.levels_extend,
cmap=cmap, norm=plotvars.norm,
alpha=alpha, transform=ccrs.PlateCarree(),
zorder=zorder, transform_first=transform_first)
else:
if np.size(field_irregular_real) > 0:
plotvars.image = mymap.tricontourf(lons_irregular_real, lats_irregular_real, field_irregular_real * fmult,
clevs, extend=plotvars.levels_extend,
cmap=cmap, norm=plotvars.norm,
alpha=alpha, transform=ccrs.PlateCarree(),
zorder=zorder)
# Block fill
if blockfill and not blockfill_ugrid:
if verbose:
print('con - adding blockfill')
two_d = False
if np.ndim(x) == 2 and np.ndim(y) == 2:
two_d = True
if isinstance(f, cf.Field):
if f.ref('grid_mapping_name:transverse_mercator', default=False):
# Special case for transverse mercator
bfill(f=f, clevs=clevs, lonlat=False, alpha=alpha, fast=blockfill_fast, zorder=zorder)
#elif orca:
elif two_d:
#bfill(f=f, clevs=clevs, lonlat=False, alpha=alpha, fast=blockfill_fast,zorder=zorder)
#bfill(x=x, y=y, f=field * fmult, clevs=clevs, lonlat=False, alpha=alpha,\
# fast=blockfill_fast, zorder=zorder, orca=True)
bfill(x=x, y=y, f=field * fmult, clevs=clevs, lonlat=False, alpha=alpha,\
fast=blockfill_fast, zorder=zorder)
else:
if f.coord('X').has_bounds() and f.coord('Y').has_bounds():
xpts = np.squeeze(f.coord('X').bounds.array[:, 0])
ypts = np.squeeze(f.coord('Y').bounds.array[:, 0])
# Add last longitude point
xpts = np.append(xpts, f.coord('X').bounds.array[-1, 1])
# Add last latitude point
ypts = np.append(ypts, f.coord('Y').bounds.array[-1, 1])
bfill(f=field_orig * fmult, x=xpts, y=ypts, clevs=clevs,