-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathacs_plotting_maps.py
4239 lines (3623 loc) · 156 KB
/
acs_plotting_maps.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
"""Standardising Australia Hazard Maps
This module plots maps to consistently present climate hazards for Australia.
It is code is designed to work with hh5 analysis3-24.04 venv"""
import datetime
import os
# import packages used in this workflow
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import image, cm
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from glob import glob
# import colormap packages
import cmaps
from matplotlib.colors import ListedColormap, BoundaryNorm, LinearSegmentedColormap
from shapely.geometry import box
import shapely
# set tolerance for simplifying geometries
# this tolerance is approx 110 m resolution for lat lon data
tolerance=0.001
import matplotlib as mpl
mpl.rcParams['hatch.linewidth'] = 0.3
plt.rcParams['savefig.facecolor']='white'
# define some standard imput for the maps
projection = ccrs.LambertConformal(
central_latitude=-24.75,
central_longitude=134.0,
cutoff=30,
standard_parallels=(-10, -40),
)
from pathlib import Path
# Set title sizes according to ACS comms guidance (in pt)
fontsize_title = 14
fontsize_subtitle = 12
fontsize_cbar = 10
fontsize_footnote = 8
padding = 5
# # Suggested colormaps and scales
# Using suggested colormaps and scales will improve the consistency across teams
# producing similar variables. This will support comparison across different plots.
# - see many colormaps here:
# https://www.ncl.ucar.edu/Document/Graphics/color_table_gallery.shtml
# This suggested colormaps are matched with possible variables to plot.
# This includes color maps for the total amount and anomalies
#ipcc colormaps from github.com/IPCC-WG1/colormaps/
cmap_mustard = LinearSegmentedColormap.from_list(
"mustard",
["#5a4c10", "#977f1b","#d3b125", "#e2c85e", "#eddd9b",]
)
cmap_mustard.set_bad(color="lightgrey")
cmap_BuGnPi = LinearSegmentedColormap.from_list('BuGnPi', np.vstack((cm.BrBG(np.linspace(0.55, 1, 30)[::-1]),[0.9,0.9,0.9,0.9], cm.pink_r(np.linspace(0.25, 1, 30)) )))
cmap_BuGnPi_r = LinearSegmentedColormap.from_list('BuGnPi_r',cmap_BuGnPi(np.linspace(0,1,256))[::-1])
cmap_dir = f"{Path(__file__).parent}/continuous_colormaps_rgb_0-1"
cmap_dict = {
"sst": cmaps.cmocean_tempo,
"sst_anom": cmaps.cmocean_balance_r,
"mhw_days": cm.YlOrRd,
"mhw_intensity": cm.hot_r,
"hot_r": cm.hot_r,
"surface_pH": cm.YlGnBu,
"surface_aragonite_sat": cmaps.cmocean_delta,
"tas": cm.Spectral_r,
"tas_anom": cm.RdBu_r,
"tas_anom_1": cm.seismic,
"tas_deciles_bwr": cm.bwr,
"EHF_days": cm.YlOrRd,
"EHF_days_1": cm.YlOrBr,
"EHF_duration": cm.hot_r,
"AFDRS_category": ListedColormap(["white", "green", "orange", "red", "darkred"]),
"ffdi_category": ListedColormap(
["green", "blue", "yellow", "orange", "red", "darkred"]
),
"fire_climate": ListedColormap(
[ "#84a19b", "#e0d7c6", "#486136", "#737932", "#a18a6e", ]
),
"fire_climate_alternative": ListedColormap(
[ "#355834", "#F1F5F2", "#14281D", "#6E633D", "#C2A878", ]
),
'tasmax_bom': ListedColormap(
[ '#E3F4FB','#C8DEE8','#91C4EA','#56B6DC','#00A2AC','#30996C',
'#7FC69A','#B9DA88','#DCE799', '#FCE850','#EACD44','#FED98E',
'#F89E64','#E67754','#D24241', '#AD283B','#832D57','#A2667A','#AB9487']
), #not colourblind safe
'tasmax' : ListedColormap(
['#014636','#016c59','#02818a','#3690c0','#67a9cf','#a6bddb',
'#d0d1e6','#ece2f0','#fff7fb','#ffffcc','#ffeda0','#fed976',
'#feb24c','#fd8d3c','#fc4e2a','#e31a1c','#bd0026','#800026',
'#510019','#2E000E']
),
"pr": cm.YlGnBu,
"pr_1": cmaps.cmocean_deep,
"pr_days": cm.Blues,
"pr_GMT_drywet": cmaps.GMT_drywet,
"pr_anom": cm.BrBG,
"pr_anom_1": cmaps.cmocean_curl,
"pr_anom_12lev": cmaps.precip_diff_12lev,
"pr_chance_extremes": cmaps.cmocean_solar_r,
"tc_days": cm.RdPu,
"tc_intensity": cm.PuRd,
"tc_days_anom": cm.PuOr,
"tc_intensity_anom": cm.PiYG_r,
"xts_freq": cmaps.cmocean_dense,
"xts_intensity": cmaps.cmocean_matter,
"xts_freq_anom": cmaps.cmocean_balance_r,
"xts_intensity_anom": cmaps.cmocean_curl_r,
"drought_severity": cm.RdYlGn, # not colorblind friendly
"drought_severity_r": cm.RdYlGn_r, # not colorblind friendly
"drought_duration": cmaps.hotres,
"drought_duration_r": cmaps.hotres_r,
"aridity": cmap_mustard,
# "aridity_anom": cmaps.NEO_div_vegetation_a, # not colorblind friendly
# "aridity_anom_r": cmaps.NEO_div_vegetation_a_r, # not colorblind friendly
"aridity_anom": cmap_BuGnPi_r,
"aridity_anom_r": cmap_BuGnPi,
"BrBu": LinearSegmentedColormap.from_list("BrBu", ["#3c320a", "#d3b125", "lightgrey", "royalblue", "navy"]),
"BuBr": LinearSegmentedColormap.from_list("BuBr", ["navy", "royalblue", "lightgrey", "#d3b125", "#3c320a"]),
"anom_BlueYellowRed": cmaps.BlueYellowRed,
"anom_BlueYellowRed_r": cmaps.BlueYellowRed_r,
"anom": cmaps.BlueWhiteOrangeRed,
"anom_r": cmaps.BlueWhiteOrangeRed_r,
"anom_b2r": cmaps.cmp_b2r,
"anom_b2r_r": cmaps.cmp_b2r_r,
"anom_coolwarm": cmaps.MPL_coolwarm,
"anom_coolwarm_r": cmaps.MPL_coolwarm_r,
"anom_deciles": cm.bwr,
"anom_deciles_r": cm.bwr_r,
"anom_veg_1": cmaps.NEO_div_vegetation_a, # not colorblind friendly
"anom_veg_1_r": cmaps.NEO_div_vegetation_a_r, # not colorblind friendly
"BuGnRd": cmaps.temp1,
"rh_19lev": cmaps.rh_19lev,
"sunshine_9lev": cmaps.sunshine_9lev,
"sunshine_diff_12lev": cmaps.sunshine_diff_12lev,
"inferno": cm.inferno,
"Oranges": cm.Oranges,
"Oranges_r": cm.Oranges_r,
"OrRd": cm.OrRd,
"Greens": cm.Greens,
"topo": cmaps.OceanLakeLandSnow,
"gmt_relief": cmaps.GMT_relief,
"ipcc_chem_div": LinearSegmentedColormap.from_list('colormap', np.loadtxt(f"{cmap_dir}/chem_div.txt")),
"ipcc_chem_seq": LinearSegmentedColormap.from_list('colormap', np.loadtxt(f"{cmap_dir}/chem_seq.txt")),
"ipcc_cryo_div": LinearSegmentedColormap.from_list('colormap', np.loadtxt(f"{cmap_dir}/cryo_div.txt")),
"ipcc_cryo_seq": LinearSegmentedColormap.from_list('colormap', np.loadtxt(f"{cmap_dir}/cryo_seq.txt")),
"ipcc_misc_div": LinearSegmentedColormap.from_list('colormap', np.loadtxt(f"{cmap_dir}/misc_div.txt")),
"ipcc_misc_seq_1": LinearSegmentedColormap.from_list('colormap', np.loadtxt(f"{cmap_dir}/misc_seq_1.txt")),
"ipcc_misc_seq_2": LinearSegmentedColormap.from_list('colormap', np.loadtxt(f"{cmap_dir}/misc_seq_2.txt")),
"ipcc_misc_seq_3": LinearSegmentedColormap.from_list('colormap', np.loadtxt(f"{cmap_dir}/misc_seq_3.txt")),
"ipcc_prec_div": LinearSegmentedColormap.from_list('colormap', np.loadtxt(f"{cmap_dir}/prec_div.txt")),
"ipcc_prec_seq": LinearSegmentedColormap.from_list('colormap', np.loadtxt(f"{cmap_dir}/prec_seq.txt")),
"ipcc_slev_div": LinearSegmentedColormap.from_list('colormap', np.loadtxt(f"{cmap_dir}/slev_div.txt")),
"ipcc_slev_seq": LinearSegmentedColormap.from_list('colormap', np.loadtxt(f"{cmap_dir}/slev_seq.txt")),
"ipcc_temp_div": LinearSegmentedColormap.from_list('colormap', np.loadtxt(f"{cmap_dir}/temp_div.txt")),
"ipcc_temp_seq": LinearSegmentedColormap.from_list('colormap', np.loadtxt(f"{cmap_dir}/temp_seq.txt")),
"ipcc_wind_div": LinearSegmentedColormap.from_list('colormap', np.loadtxt(f"{cmap_dir}/wind_div.txt")),
"ipcc_wind_seq": LinearSegmentedColormap.from_list('colormap', np.loadtxt(f"{cmap_dir}/wind_seq.txt")),
"acs_geophysical_biochemical_div1": cm.PRGn,
"acs_geophysical_biochemical_div2": cm.PiYG,
"acs_geophysical_biochemical_seq1": cm.BuGn,
"acs_geophysical_biochemical_seq2": cm.RdPu,
"acs_geophysical_biochemical_seq3": cm.Greys,
"acs_precipitation_div1": cm.RdBu,
"acs_precipitation_div2": cm.BrBG,
"acs_precipitation_seq1": cm.Blues,
"acs_precipitation_seq2": cm.GnBu,
"acs_precipitation_seq3": cm.PuBuGn,
"acs_temperature_div1": cm.RdYlBu_r,
"acs_temperature_div2": cm.PuOr_r,
"acs_temperature_seq1": cm.YlOrRd,
"acs_temperature_seq2": cm.Reds,
"acs_temperature_seq3": cm.YlGnBu,
}
# Here are some suggestions for the ticks/ scale for some variables.
# Some scales are taken from climate maps on bom.gov.au/climate
tick_dict = {
"pr_annual": [0, 50, 100, 200, 300, 400, 600, 1000, 1500, 2000, 3000, 6000],
"pr_6mon": [0, 50, 100, 200, 300, 400, 600, 900, 1200, 1800, 2400, 6000],
"pr_3mon": [0, 10, 25, 50, 100, 200, 300, 400, 600, 800, 1200, 2500],
"pr_mon": [0, 1, 5, 10, 25, 50, 100, 200, 300, 400, 600, 1200],
"pr_hour": [0, 1, 2, 5, 10, 15, 20, 30, 50, 75, 100, 200,],
"pr_days": [0, 2, 3, 5, 10, 20, 30, 40, 50, 75, 100, 125, 150, 175],
"pr_anom_mon": [ -1000, -400, -200, -100, -50, -25, -10, 0, 10, 25, 50, 100, 200, 400, 1000,],
"pr_anom_3mon": [ -2000, -600, -400, -200, -100, -50, -25, 0, 25, 50, 100, 200, 400, 600, 2000,],
"pr_anom_6mon": [ -3000, -1200, -800, -400, -200, -100, -50, 0, 50, 100, 200, 400, 800, 1200, 3000,],
"pr_anom_ann": [ -4000, -1800, -1200, -800, -400, -200, -100, 0, 100, 200, 400, 800, 1200, 1800, 4000,],
"pr_diff_mon": [ -1000, -400, -200, -100, -50, -25, -10, 10, 25, 50, 100, 200, 400, 1000,],
"pr_diff_ann": [ -3000, -1800, -1200, -800, -400, -200, -100, 100, 200, 400, 800, 1200, 1800, 3000,],
"frost_days": [0, 10, 20, 30, 40, 50, 75, 100, 150, 300],
"frost_days_mon": [0, 2, 5, 10, 15, 20, 25, 31],
"tas": np.arange(-9, 52, 3),
"tas_anom_day": np.arange(-14, 14.1, 2),
"tas_anom_mon": np.arange(-7, 7.1, 1),
"tas_anom_ann": np.arange(-3.5, 3.6, 0.5),
"apparent_tas": np.arange(-6, 42, 3),
"percent": np.arange(0, 101, 10),
"xts_freq": [0.00, 0.005, 0.01, 0.02, 0.03, 0.05, 0.07, 0.10, 0.12, 0.15],
"fire_climate_ticks": [ 100, 101, 102, 103, 104, ],
"fire_climate_labels": [
"Tropical Savanna",
"Arid grass \nand woodland",
"Wet Forest",
"Dry Forest",
"Grassland",
],
"aridity_index_ticks": [0.0, 0.05, 0.2, 0.5, 0.65],
"aridity_index_labels": ["Hyper-arid", "Arid", "Semi-arid", "Dry sub-humid"],
}
# # Load the State and Region shape files
class RegionShapefiles:
"""Load and return a shapefile based on its name."""
def __init__(self, path, shapefiles):
"""Create an instance of the RegionShapefiles class.
Parameters
----------
path : str
The path to the shapefiles directory.
shapefiles : list
A list of shapefile names to load.
"""
self.path = path
self.shapefiles = shapefiles
self._regions_dict = {}
def __call__(self, name):
"""Retrieve the shapefile for the given region name.
Parameters
----------
name : str
The name of the region to retrieve.
Returns
-------
GeoDataFrame or GeoSeries
The shapefile data for the specified region.
"""
if name not in self._regions_dict:
if name in self.shapefiles:
regions = gpd.read_file(glob(f"{self.path}/{name}/*.shp")[0]).to_crs(crs = "GDA2020")
regions[["geometry"]] =shapely.simplify(regions[["geometry"]], tolerance)
self._regions_dict[name] = regions
elif name == "not_australia":
# Define a white mask for the area outside of Australian land
# We will use this to hide data outside the Australian land borders.
# note that this is not a data mask,
# the data under the masked area is still loaded and computed, but not visualised
base_name = name[4:] # Remove 'not_' prefix
base_region = self(base_name).copy().to_crs(crs = "GDA2020")
base_region[["geometry"]] =shapely.simplify(base_region[["geometry"]], tolerance)
# This mask is a rectangular box around the maximum land extent of Australia
# with a buffer of 20 degrees on every side,
# with the Australian land area cut out, only the ocean is hidden.
not_region = gpd.GeoSeries(
data=[
box(*box(*base_region.total_bounds).buffer(20).bounds
).difference(base_region["geometry"].values[0])
],
crs=ccrs.PlateCarree(),
)
self._regions_dict[name] = not_region
else:
raise ValueError(f"Shapefile for region '{name}' not found.")
return self._regions_dict[name]
def __getitem__(self, name):
if name not in self._regions_dict:
self(name)
return self._regions_dict[name]
def __setitem__(self, name, value):
self._regions_dict[name] = value
def keys(self):
return self._regions_dict.keys()
def __len__(self):
return len(self._regions_dict)
def __repr__(self):
return repr(self._regions_dict)
def update(self, *args, **kwargs):
self._regions_dict.update(*args, **kwargs)
# Define the path and shapefiles
# These will be used for state boundaries, LGAs, NRM, etc
PATH = "/g/data/ia39/aus-ref-clim-data-nci/shapefiles/data"
shapefile_list = ["aus_local_gov",
"aus_states_territories",
"australia",
"broadacre_regions",
"NCRA_Marine_region",
"ncra_regions",
"NCRA_regions_coastal_waters_GDA94",
"nrm_regions",
"plantations"]
# Create an instance of the RegionShapefiles class
regions_dict = RegionShapefiles(PATH, shapefile_list)
# define a white mask for the area outside of Australian land
# We will use this to hide data outside of the Australian land borders.
# note that this is not a data mask,
# the data under the masked area is still loaded and computed, but not visualised
australia = regions_dict["australia"]
# This mask is a rectangular box around the maximum land extent of Australia
# with a buffer of 10 degrees on every side,
# with the Australian land area cut out so only the ocean is hidden.
not_australia = regions_dict["not_australia"]
def crop_cmap_center(cmap, ticks, mid, extend=None):
"""
This function is used to align the centre of a colormap (cmap)
to a specified midpoint (mid). Allows the cmap to be normalised on
the specified ticks with cbar_extend arrows taken into account.
Intended for divergent colormaps that show anomalies that are mostly
all positive (or negative), for example, temperature anomalies in
future climate projections.
The shorter side of the colormap is cropped and not stretched - unlike
matplotlib.colors.TwoSlopeNorm.
Parameters
-------
cmap: matplotlib colormap
original colormap to crop
ticks: list or np.array
list or array of numerical ticks
mid: float
Set the midpoint value of the colormap.
For example, 0.0
extend: {'neither', 'both', 'min', 'max'}
Make pointed end(s) for out-of-range values (unless 'neither').
These are set for a given colormap using the colormap set_under
and set_over methods.
Returns
-------
matplotlib.colors.LinearSegmentedColormap
"""
ticks=np.array(ticks)
# number of color segments:
below = (ticks<mid).sum()
above = (ticks>mid).sum()
if extend =="both" or extend == "max":
above=above+1
if extend =="both" or extend =="min":
below=below+1
#total segments
N = below+above
#porportion below mid point
prop_below_mid = below/(max(below,above)+1)
# propotion above mid point
prop_above_mid = above/(max(below,above)+1)
bounds = np.linspace(0.5*(1-prop_below_mid),
0.5*(1+prop_above_mid),
N+1)
between_bounds = [(n0+n1)/2 for n0, n1 in zip(bounds[:-1], bounds[1:])]
new_cmap_list = cmap(between_bounds)
new_cmap = LinearSegmentedColormap.from_list("new_cmap",new_cmap_list, N)
return new_cmap
# Define subfunctions for different parts of the plotting
# so that they can be reused for single panel and multi panel plots
def plot_data(regions=None,
data=None,
station_df = None,
markersize=None,
stippling=None,
shading=None,
xlim=(114, 162),
ylim=(-43.5, -7.5),
cmap=cm.Greens,
cbar_extend="both",
ticks=None,
tick_interval=1,
tick_labels=None,
contourf=False,
contour=False,
ax=None,
subtitle = "",
subtitle_xy = None,
facecolor="none",
edgecolor="k",
area_linewidth=0.3,
coastlines=False,
mask_not_australia = True,
mask_australia=False,
agcd_mask=False,
select_area = None,
vcentre=None,
):
"""This function takes one axis and plots the hazard data to one map of Australia.
This function takes gridded "data" and/or a "station_df" dataframe and "regions" shapefiles
to visualise hazard data from a 2D Xarray data array and plots the data on a map
of Australia with the regions outlines.
Parameters
----------
regions: geopandas.GeoDataFrame
region geometries for regions/states/catchments etc
data: xr.DataArray
a 2D xarray DataArray which has already computed the
average, sum, anomaly, metric or index you wish to visualise.
This function is resolution agnostic.
station_df: pd.DataFrame, optional
a pandas.DataFrame with columns ["lon", "lat", variable].
If station_df is given, then variable values are represented as dots on
the mapaccordingg to the lat lon coordinates and coloured according to
cmap colors and ticks.
markersize: int, optional
default None. If None then the markersize will adjust to the size of the
figure and the number of markers in the plot such that when there are
many markers and the figure is small, the markersize is smaller.
stippling: xr.DataArray, optional
a True/False mask to define regions of stippling hatching.
Intended to show information such as "model agreement for direction of change".
shading: xr.DataArray, optional
A list of True/False masks to define regions of dotted stippling
for each subplot. ["....."]
Intended to obscure data outside the area of interest.
xlim: tuple, optional
longitude min and max of the plot area.
default is cropped to Australian continent xlim=(114, 162)
ylim: tuple, optional
latitude min and max of the plot area.
default is cropped to Australian continent ylim=(-43.5, -7.5),
cmap:
defines the colormap used for the data.
See cmap_dict for suggested colormaps.
Default cmap set to cm.Greens.
Please choose appropriate colormaps for your data.
cbar_extend: one of {'neither', 'both', 'min', 'max'}.
eg "both" changes the ends of the colorbar to arrows to indicate that
values are possible outside the scale shown.
If contour or contourf is True, then cbar_extend will be overridden to "none".
ticks: list or arraylike
Define the ticks on the colorbar. Define any number of intervals.
This will make the color for each interval one discrete color,
instead of a smooth color gradient.
If None, linear ticks will be auto-generated to fit the provided data.
tick_interval: int
Default 1
For showing every second tick label, set tick_interval=2
tick_labels: list
Labels for categorical data.
If tick_labels is used, then pcolormesh is used to plot data
and does not allow contour or contourf to be used.
Tick labels will correspond to the ticks.
contourf: bool
if True then the gridded data is visualised as smoothed filled contours.
Default is False.
Use with caution when plotting data with negative and positive values;
Check output for NaNs and misaligned values.
High resolution data is slow to compute.
contour: bool
if True then the gridded data is visualised as smoothed unfilled grey contours.
Default is False.
High resolution data is slow to compute.
Using both contourf and contour results in smooth filled contours
with grey outlines between the color levels.
ax: matplotlib.axes.Axes
Axes object of existing figure to put the plotting.
subtitle: str
default ""
Intended to label global warming levels for subplots eg. "GWL 1.2"
subtitle_xy: tuple, optional
(x, y) location of subtitle relative to transAxes.
defines the top left location for the subtitle.
facecolor: color
color of land when plotting the regions without climate data and select_area is None.
facecolor recommendations include "white", "lightgrey", "none".
default is "none"
edgecolor: color
defines the color of the state/region borders.
edgecolor recommendations include "black" and "white".
mask_not_australia: boolean
decides whether or not the area outside of Australian land is hidden
under white shape.
Default is True.
mask_australia: boolean
decides whether or not Australian land is hidden under white shape.
Eg, use when plotting ocean only.
Default is False.
agcd_mask: boolean
If True, applies a ~5km mask for data-sparse inland areas of Australia.
Default is False.
area_linewidth: float, optional
the width of state/region borders only. All other linewidths are hardcoded.
coastlines: boolean
If True, add cartopy coastlines for all coasts (not just Australia).
Default is False.
select_area: list
If None, then don't add region borders geometries.
vcentre: float, eg 0
default is None.
Align centre of colormap to this value.
Intended for using a divergent colormap with uneven number of ticks
around the centre, eg for future temperature anomalies with a larger
positive range compared to the negative range.
Returns
-------
ax, norm, cont, middle_ticks
ax, the matplotlib.axes.Axes with the plot
norm, the normalisation for the colormap and plotted contours according to ticks
cont, the QuadContourSet or QuadMesh of the plotted gridded data
middle_ticks, the location to label categorical tick labels
"""
if vcentre is not None:
cmap = crop_cmap_center(cmap, ticks, vcentre, extend=cbar_extend)
ax.set_extent([xlim[0], xlim[1], ylim[0], ylim[1]])
middle_ticks=[]
# set norm from ticks and cbar_extend
if ticks is None:
norm = None
else:
# if ticks are labelled or if there is one more tick than tick labels,
# do the usual normalisation
if tick_labels is None or (len(tick_labels) == len(ticks) - 1):
norm = BoundaryNorm(ticks, cmap.N+1, extend = cbar_extend)
if tick_labels is not None:
middle_ticks = [
(ticks[i + 1] + ticks[i]) / 2 for i in range(len(ticks) - 1)
]
else:
middle_ticks = []
else:
middle_ticks = [
(ticks[i + 1] + ticks[i]) / 2 for i in range(len(ticks) - 1)
]
outside_bound_first = [ticks[0] - (ticks[1] - ticks[0]) / 2]
outside_bound_last = [ticks[-1] + (ticks[-1] - ticks[-2]) / 2]
bounds = outside_bound_first + middle_ticks + outside_bound_last
norm = BoundaryNorm(bounds, cmap.N, extend = "neither")
# define cont
if data is None:
cont=None
else:
data = data.squeeze()
if agcd_mask:
# mask data where observations are sparse
directory = "/g/data/ia39/aus-ref-clim-data-nci/shapefiles/masks/AGCD-05i"
agcd = xr.open_dataset(f"{directory}/mask-fraction_agcd_v1-0-2_precip_weight_1960_2022.nc").fraction
data = data.where(agcd>=0.8)
# plot the hazard data
if contourf and tick_labels is None:
if data.max()>=0 and data.min()<=0:
print("Using contourf to plot data. Use with caution and check output for data crossing zero")
cont = ax.contourf(
data.lon,
data.lat,
data,
cmap=cmap,
norm=norm,
levels=ticks,
extend=cbar_extend,
zorder=2,
transform=ccrs.PlateCarree(),
)
else:
cont = ax.pcolormesh(
data.lon,
data.lat,
data,
cmap=cmap,
norm=norm,
zorder=2,
transform=ccrs.PlateCarree(),
)
if contour and tick_labels is None:
ax.contour(
data.lon,
data.lat,
data,
colors="grey",
norm=norm,
levels=ticks,
extend=cbar_extend,
linewidths=0.2,
zorder=3,
transform=ccrs.PlateCarree(),
)
# for station data
if station_df is not None:
# assuming columns are named "lon", "lat", variable,
gdf = gpd.GeoDataFrame(
station_df, geometry=gpd.points_from_xy(station_df.lon, station_df.lat), crs=ccrs.PlateCarree()
)
var = gdf.columns[[2]][0]
# norm = BoundaryNorm(ticks, cmap.N, extend=cbar_extend)
cont = ax.scatter(x=station_df.lon,
y=station_df.lat,
s=markersize,
c=station_df[var],
# edgecolors="k",
alpha = 0.8,
zorder=7,
transform=ccrs.PlateCarree(),
cmap= cmap,
norm = norm)
if stippling is not None:
ax.contourf(stippling.lon,
stippling.lat,
stippling,
alpha=0,
hatches = ["","xxxxxx"],
zorder=4,
transform=ccrs.PlateCarree(),
)
if shading is not None:
ax.contourf(shading.lon,
shading.lat,
shading,
alpha=0,
hatches = ["","....."],
zorder=4,
transform=ccrs.PlateCarree(),
)
# cover area outside australia land area eg mask ocean
if mask_not_australia:
# inside the shape, fill white
ax.add_geometries(
not_australia,
crs=ccrs.PlateCarree(),
facecolor="white",
linewidth=0,
zorder=5,
)
# cover australia land area eg for ocean data
if mask_australia:
# inside the shape, fill white
ax.add_geometries(
australia["geometry"],
crs=ccrs.PlateCarree(),
facecolor="lightgrey",
linewidth=0.3,
edgecolor="k",
zorder=4,
)
if select_area is None:
# add region borders unless you have selected area
ax.add_geometries(
regions["geometry"],
crs=ccrs.PlateCarree(),
facecolor=facecolor,
edgecolor=edgecolor,
linewidth=area_linewidth,
zorder=6,
)
if coastlines:
try:
ax.add_feature(cfeature.BORDERS, zorder=5, linewidth=area_linewidth*0.8,)
except:
print("could not download borders")
ax.coastlines(resolution = "10m", zorder=5, linewidth=area_linewidth*0.8,)
# subtitle
if subtitle_xy is None:
subtitle_xy = (0.02, 0.02)
ax.text(
x=subtitle_xy[0],
y=subtitle_xy[1],
s=subtitle,
fontsize=fontsize_subtitle,
horizontalalignment="left",
verticalalignment="bottom",
transform=ax.transAxes,
zorder=10,
wrap=True,
)
return ax, norm, cont, middle_ticks
def plot_cbar(cont=None,
norm=None,
ax=None,
cbar_extend=None,
cbar_label=None,
ticks=None,
tick_interval=1,
tick_labels=None,
middle_ticks=[],
cax_bounds =None,
contour=False,
location=None,
rotation=None,
):
"""This function defines and plots the colorbar.
It takes crucial information from the plot_data function.
Parameters
----------
cont:
output from matplotlib plotting
norm:
normalisation
ax: matplotlib.axes.Axes
Axes to put the colorbar
cbar_extend: one of {'neither', 'both', 'min', 'max'}.
eg "both" changes the ends of the colorbar to arrows to indicate that
values are possible outside the scale shown.
cbar_label: str
Title for colorbar. Include name of metric and [units]
ticks: list or array
numerical location of ticks
tick_interval: int
Default 1
For showing every second tick label, set tick_interval=2
tick_labels: list
If categorical data, these labels go inbetween the numerical bounds set by ticks
middle_ticks: list
If categorical data, this specifies the location of the tick labels.
Typically in the middle of the bounds set by ticks
cax_bounds: [left, bottom, width, height]
Colorbar axes bounds relative to ax
contour: bool
draw lines on colorbar if True
Default is False
location: {"top", "bottom", "left", "right"}
location of the colorbar. Defaults to right.
rotation: [-360,360]
rotation of tick labels in degrees. Set to 0 for horizontal.
Returns
-------
matplotlib.colorbar
"""
if cax_bounds is not None:
cax = ax.inset_axes(cax_bounds)
else:
cax=None
cbar = None
if norm is None:
return cbar
if tick_labels is None:
cbar = plt.colorbar(
cont,
ax=ax,
extend=cbar_extend,
cax=cax,
ticks=ticks,
norm=norm,
location=location,
fraction=0.046,
pad=0.04
)
# only label ticks at specified tick intervals
[l.set_visible(False) for (i,l) in enumerate(cbar.ax.xaxis.get_ticklabels()) if i % tick_interval != 0]
[l.set_visible(False) for (i,l) in enumerate(cbar.ax.yaxis.get_ticklabels()) if i % tick_interval != 0]
else:
# for categorical data
cbar = plt.colorbar(
cont,
ax=ax,
extend='neither',
cax=cax,
ticks=ticks,
norm=norm,
drawedges=True,
location=location,
fraction=0.046,
pad=0.04
)
if location=="bottom":
if len(ticks) == len(tick_labels):
cbar.ax.set_xticks(ticks, tick_labels, wrap=True, verticalalignment="top")
elif len(middle_ticks) == len(tick_labels):
cbar.ax.set_xticks(middle_ticks, tick_labels, wrap=True, verticalalignment="top")
else:
if len(ticks) == len(tick_labels):
cbar.ax.set_yticks(ticks, tick_labels, wrap=True)
elif len(middle_ticks) == len(tick_labels):
cbar.ax.set_yticks(middle_ticks, tick_labels, wrap=True)
cbar.ax.tick_params(labelsize=8)
if contour and tick_labels is None:
cbar.add_lines(cont)
# Label colorbar
if cbar is not None:
cbar.ax.set_title(cbar_label,
zorder=10,
loc="center",
fontsize=fontsize_cbar,
fontstretch="normal",
verticalalignment="bottom",
wrap=True)
cbar.ax.tick_params(rotation=rotation)
return cbar
def plot_select_area(select_area=None,
ax=None,
xlim=None,
ylim=None,
regions=None,
land_shadow=False,
area_linewidth=0.3,
):
"""This function takes a list of named areas to plot and adjusts
the limits of the axis.
Parameters
----------
select_area: list
list of selected areas to plot. Must be name of area in regions.
ax: matplotlib.axes.Axes
axis
xlim:
longitude limits to use if select_area is None
ylim:
latitude limits to use if select_area is None
regions: geopandas.GeoDataFrame
region border data, must contain a column name with "NAME" in it
to select those areas.
land_shadow: bool
whether or not to shade in the non-selected areas. Can help
visualise land-sea borders.
area_linewidth: float
default 0.3
linewidth of area edges. Larger values have thicker borders.
Returns
-------
matplotlib.axes.Axes
"""
if select_area is None:
ax.set_extent([xlim[0], xlim[1], ylim[0], ylim[1]], crs=ccrs.PlateCarree())
pass
else:
assert isinstance(select_area, list), "select_area must be a list"
# select state
name_column = [name for name in regions.columns if "NAME" in name.upper()][0]
area = regions.loc[regions[name_column].isin(select_area)]
area= area.to_crs(crs = "GDA2020")
area[["geometry"]] =shapely.simplify(area[["geometry"]], tolerance)
map_total_bounds = area.total_bounds
minx, miny, maxx, maxy = map_total_bounds
mid_x = (minx + maxx) / 2
mid_y = (miny + maxy) / 2
max_range = np.max([(maxy - miny), (maxx - minx)])
buffer = 0.1 * max_range
not_area = gpd.GeoSeries(
data=[
box(*box(*map_total_bounds).buffer(10 * buffer).bounds).difference(
area.dissolve()["geometry"].values[0]
)
],
crs=ccrs.PlateCarree(),
)
# mask outside selected area
if land_shadow:
# show land as light grey
ax.add_geometries(not_area,
crs=ccrs.PlateCarree(),
facecolor="lightgrey",
linewidth=area_linewidth,
edgecolor="k",
zorder=4)
else:
# mask white
ax.add_geometries(not_area,
crs=ccrs.PlateCarree(),
facecolor="white",
linewidth=area_linewidth,
edgecolor="k",
zorder=4)
ax.set_extent([mid_x - 0.6 * max_range,
mid_x + 0.8 * max_range,
mid_y - 0.7 * max_range,
mid_y + 0.7 * max_range],
crs=ccrs.PlateCarree())
return ax
def plot_titles(title="title",
date_range = "DD Mon YYYY to DD Mon YYYY",
baseline = None,
dataset_name= None,
issued_date=None,
watermark= None,
watermark_color="r",
ax=None,
text_xy = None,
title_ha = "left",
show_copyright = True,):
"""
Set the plot title and axis labels
Parameters
----------
title: str
Main text. Size 14, bold. Location set by text_xy["title"].
Horizontal alignment set by title_ha.
date_range: str
Text under title. Size 12. Horizontal alignment set by title_ha.
Intended for data date range. Also can be used for any subtitle.
Default "DD Mon YYYY to DD Mon YYYY" to indicate prefered date format.