Source Code

Backend Class

   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
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
#cython: linetrace=True
#distutils: define_macros=CYTHON_TRACE_NOGIL=1
#cython: language_level=3
"""
AMSIMP Backend Class. For information about this class is described below.

Copyright (C) 2020 AMSIMP

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program 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 General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <https://www.gnu.org/licenses/>.
"""

# -----------------------------------------------------------------------------------------#

# Importing Dependencies
import os, sys
from datetime import datetime
import socket
from amsimp.download cimport Download
from amsimp.download import Download
import numpy as np
from astropy import constants as constant
from astropy import units
from astropy.units.quantity import Quantity
cimport numpy as np
from cpython cimport bool
import matplotlib.pyplot as plt
from matplotlib import ticker
import cartopy.crs as ccrs
from cartopy.util import add_cyclic_point
import iris
from iris.coords import DimCoord
from iris.cube import Cube
from iris.coord_systems import GeogCS

# -----------------------------------------------------------------------------------------#

cdef class Backend(Download):
    """
    AMSIMP Backend Class - This is the base class for AMSIMP, as such, all
    other classes within AMSIMP inherit this class, either directly or
    indirectly.

    For methods to be included in this class they must meet one of the following
    criteria:
    (1) they are considered essential components of the software.
    (2) the output of these methods are generated solely from the methods
    classified as (1).
    (3) these methods import data from the Global Data Assimilation System.
    This data is retrieved from the AMSIMP data repository, which is updated 
    every hour.
    (4) they don't classify nicely into any other class.
    (5) they offer a visualisation of any of the methods found in this class.
    
    Below is a list of the methods included within this class, with a short
    description of their intended purpose and a bracketed number signifying
    which of the above criteria they meet. Please see the relevant class methods
    for more information. Please note that the unit information with AMSIMP
    is provided by the unit module within astropy.

    latitude_lines ~ generates a NumPy array of latitude lines (1). The unit
    of measurement is degrees.
    longitude_lines ~ generates a NumPy array of longitude lines (1). The unit
    of measurement is degrees.
    pressure_surfaces ~ generates a NumPy array of constant pressure
    surfaces (1). The unit of measurement is hectopascals. This is the
    isobaric coordinate system.

    coriolis_parameter ~ generates a NumPy arrray of the Coriolis parameter
    at various latitudes of the Earth's surface (2). The unit of measurement
    is radians per second.

    geopotential_height ~ outputs a NumPy array of geopotential height (3). The
    unit of measurement is metres.
    relative_humidity ~ outputs a NumPy array of relative humidity (3). The unit
    of measurement is percentage.
    temperature ~ outputs a NumPy array of temperature (3). The unit of
    measurement is Kelvin.
    remove_all_files ~ this method deletes any file created by the software,
    whilst downloading data from the AMSIMP Initial Atmospheric Conditions
    Data Repository (3). It does not, however, delete any forecast generated by,
    the, simulate, method. Please note that you may specifcy whether you would like
    to save such files throught the utilisation of the parameter, remove_files,
    when the class is initialised.
    
    pressure_thickness ~ outputs a NumPy array of atmospheric pressure
    thickness (4). The unit of measurement is metres.
    troposphere_boundaryline ~ generates a NumPy array of the
    troposphere - stratosphere boundary line (4). The unit of measurement is
    metres.
    
    longitude_contourf ~ generates a contour plot for a desired atmospheric
    process, with the axes being latitude, and longitude (5).
    altitude_contourf ~ generates a contour plot for a desired atmospheric
    process, with the axes being latitude, and altitude (5).
    thickness_contourf ~ generates a contour plot for pressure thickness,
    with the axes being latitude and longitude. This plot is then layed
    on top of a EckertIII global projection (5).
    """

    # Define units of measurement for AMSIMP.
    units = units

    def __cinit__(
            self,
            int delta_latitude=5,
            int delta_longitude=5,
            forecast_length=72, 
            delta_t=2, 
            bool ai=True, 
            data_size=150, 
            epochs=200, 
            input_date=None, 
            bool input_data=False,
            psurfaces=None,
            lat=None,
            lon=None,
            height=None, 
            temp=None, 
            rh=None, 
            u=None, 
            v=None,
            dict constants={
                "sidereal_day": (23 + (56 / 60)) * 3600,
                "angular_rotation_rate": ((2 * np.pi) / ((23 + (56 / 60)) * 3600)),
                "planet_radius": constant.R_earth.value,
                "planet_mass": constant.M_earth.value,
                "specific_heat_capacity_psurface": 718,
                "gravitational_acceleration": 9.80665,
                "planet": "Earth"
            }
        ):
        """
        The parameters, delta_latitude and delta_longitude, defines the
        horizontal resolution between grid points within the software. The
        parameter values must be between 1 and 10 degrees if you are utilising
        the default initial conditions included within the software. Defaults
        to a value of 5 degrees.

        The parameter, input_date, allows the end-user to specify the initial
        date at which the forecast is generated from. Please note that
        currently you cannot specify dates before Janurary 2020. Please also
        note that the latest data available to the software is typically
        from three days ago. This is due to the rate at which the NOAA
        updates the data on the Global Data Assimilation System website. This
        feature is intended for situations where the default initial conditions
        included within the software are being used.

        The parameter, input_data, is a boolean value, and when set to True,
        it will allow the end-user to provide their own initial conditions
        to the software. Please note that this feature has not been fully
        tested as of this moment, and such is experimental for the time
        being.

        The parameter, ai, is a boolean value, and when set to True, a 
        recurrent neural network will be trained on the amount of days 
        specified by the parameter, data_size. A prediction by the
        recurrent neural network will be generated, and will be 
        combined with a traditional physical model of the atmosphere.
        It can also be used independently. 
        """
        # Make the aforementioned variables available else where in the class.
        self.delta_latitude = delta_latitude
        self.delta_longitude = delta_longitude
        self.input_date = input_date
        self.input_data = input_data

        # The date at which the initial conditions was gathered (i.e. how
        # recent the data is).
        if self.input_date == None:
            data_date_file = self.download(
                "https://github.com/amsimp/initial-conditions/raw/master/date.npy", 
                bar=False,
            )
            data_date = np.load(data_date_file)
            os.remove(data_date_file)
            date = datetime(
                int(data_date[2]),
                int(data_date[1]),
                int(data_date[0]),
                int(data_date[3]),
            )
            self.date = date
        else:        
            # Ensure input_date is of the type "datetime.datetime".
            if type(self.input_date) != datetime:
                raise Exception(
                    "input_date must of the type 'datetime.datetime'. The value of input_date was: {}".format(
                        self.input_date
                    )
                )
            date = self.input_date
            self.date = date

        # Ensure input_data is a boolean value.
        if not isinstance(self.input_data, bool):
            raise Exception(
                "input_data must be a boolean value. The value of input_data was: {}".format(
                    self.input_data
                )
            )

        # Function to ensure that the input data is 3 dimensional.
        def dimension(input_variable):
            if np.ndim(input_variable) != 3:
                raise Exception(
                    "All input data variables (height, rh, temp, u, v) must be a 3 dimensional."
                )

        # Function to check if a function is a NumPy array, or a list.
        def type_check(input_variable):
            if type(input_variable) == Quantity:
                pass
            elif type(input_variable) == np.ndarray:
                pass
            elif type(input_variable) == list:
                pass
            else:
                raise Exception(
                    "All input data variables must be either NumPy arrays, or lists."
                )
        
        # Function to check if grid input variables are 1 dimensional.
        def dimension_grid(input_variable):
            if np.ndim(input_variable) != 1:
                raise Exception(
                    "All grid input variables (psurfaces, lat, lon) must be a 1 dimensional."
                )

        # Function to convert input lists into NumPy arrays.
        def list_to_numpy(input_variable):
            if type(input_variable) == list:
                input_variable = np.asarray(input_variable)
            
            return input_variable

        if self.input_data:
            # Ensure input data variables are either NumPy arrays, or lists.
            type_check(psurfaces)
            type_check(lat)
            type_check(lon)
            type_check(height)
            type_check(rh)
            type_check(temp)
            type_check(u)
            type_check(v)

            # Check if grid input variables are 1 dimensional.
            dimension_grid(psurfaces)
            dimension_grid(lat)
            dimension_grid(lon)

            # Check if input data is 3 dimensional.
            dimension(height)
            dimension(rh)
            dimension(temp)
            dimension(u)
            dimension(v)

            # Convert input data lists to NumPy arrays.
            list_to_numpy(psurfaces)
            list_to_numpy(lat)
            list_to_numpy(lon)
            list_to_numpy(height)
            list_to_numpy(rh)
            list_to_numpy(temp)
            list_to_numpy(u)
            list_to_numpy(v)

            # Add units to input variables.
            # psurfaces variable.
            if type(psurfaces) != Quantity:
                psurfaces = psurfaces * units.hPa
            # lat variable.
            if type(lat) != Quantity:
                lat = lat * units.deg
            # lon variable.
            if type(height) != Quantity:
                lon = lon * units.deg
            # height variable.
            if type(height) != Quantity:
                height = height * units.m
            # rh variable.
            if type(rh) != Quantity:
                rh = rh * units.percent
            # temp variable.
            if type(temp) != Quantity:
               temp = temp * units.K
            # u variable.
            if type(u) != Quantity:
               u = u * (units.m / units.s)
            # v variable.
            if type(v) != Quantity:
               v = v * (units.m / units.s)

        # Ensure self.delta_latitude is between 1 and 10 degrees.
        if not input_data:
            if self.delta_latitude > 10 or self.delta_latitude < 1:
                raise Exception(
                    "delta_latitude must be a positive integer between 1 and 10. The value of delta_latitude was: {}".format(
                        self.delta_latitude
                    )
                )

            # Ensure self.delta_longitude is between 1 and 10 degrees.
            if self.delta_longitude > 10 or self.delta_longitude < 1:
                raise Exception(
                    "delta_longitude must be a positive integer between 1 and 10. The value of delta_longitude was: {}".format(
                        self.delta_longitude
                    )
                )

            folder = "https://github.com/amsimp/initial-conditions/raw/master/initial_conditions/"

            # Define date.
            year = self.date.year
            month = self.date.month
            day = self.date.day
            hour = self.date.hour

            # Adds zero before single digit numbers.
            if day < 10:
                day = "0" + str(day)

            if month < 10:
                month =  "0" + str(month)

            if hour < 10:
                hour = "0" + str(hour)

            # Converts integers to strings.
            day = str(day)
            month = str(month)
            year = str(year)
            hour = str(hour)

            folder += (
                year + "/" + month + "/" + day + "/" + hour + "/"
            )
            # The url to the NumPy pressure surfaces file stored on the AMSIMP
            # Initial Conditions Data repository.
            url = folder + "initial_conditions.nc"

            # Download the NumPy file and store the NumPy array into a variable.
            try:
                cube = iris.load("initial_conditions.nc")
            except OSError:  
                fname = self.download(url)
                cube = iris.load(fname)
            
            height = cube.extract("geopotential_height")[0]
            temp = cube.extract("air_temperature")[0]
            rh = cube.extract("relative_humidity")[0]
            u = cube.extract("x_wind")[0]
            v = cube.extract("y_wind")[0]

        # Constants dictionary.
        # Sidereal day.
        if type(constants["sidereal_day"]) != Quantity:
            constants["sidereal_day"] = constants["sidereal_day"] * units.s
        # Angular rotation rate.
        if type(constants["angular_rotation_rate"]) != Quantity:
            constants["angular_rotation_rate"] = constants["angular_rotation_rate"] * (
                units.rad / units.s
            )
        # Planet radius.
        if type(constants["planet_radius"]) != Quantity:
            constants["planet_radius"] = constants["planet_radius"] * units.m
        # Planet mass.
        if type(constants["planet_mass"]) != Quantity:
            constants["planet_mass"] = constants["planet_mass"] * units.kg
        # Specific heat capacity on a constant pressure surface.
        if type(constants["specific_heat_capacity_psurface"]) != Quantity:
            constants["specific_heat_capacity_psurface"] = constants[
                "specific_heat_capacity_psurface"
            ] * (units.J / (units.kg * units.K))
        # Gravitational acceleration.
        if type(constants["gravitational_acceleration"]) != Quantity: 
            constants["gravitational_acceleration"] = constants["gravitational_acceleration"] * (
                units.m / (units.s ** 2)
            )
        # Planet
        planet = constants["planet"].lower()
        planet = planet.replace(" ", "")
        planet = planet.capitalize()
        # Gas constant
        R = 287.05799596 * (units.J / (units.kg * units.K))
        constants["gas_constant"] = R
        # Universal gravitational constant.
        G = constant.G
        G = G.value * G.unit
        constants["universal_gravitational_constant"] = G

        # Make the input data variables available else where in the class.
        self.psurfaces = psurfaces
        self.lat = lat
        self.lon = lon
        self.input_height = height
        self.input_rh = rh
        self.input_temp = temp
        self.input_u = u
        self.input_v = v
        self.constants = constants

        # Predefined Constants.
        # Angular rotation rate of the planet.
        self.sidereal_day = constants["sidereal_day"]
        self.Upomega = constants["angular_rotation_rate"]
        # Mean radius of the planet.
        self.a = self.constants["planet_radius"]
        # Mass of the planet.
        self.M = self.constants["planet_mass"]
        # The specific heat capacity on a constant pressure surface.
        self.c_p = self.constants["specific_heat_capacity_psurface"]
        # Gravitational acceleration.
        self.g = self.constants["gravitational_acceleration"]
        # Planet
        self.planet = self.constants["planet"]
        # Gas Constant.
        self.R = constants["gas_constant"]
        # Universal gravitational constant.
        self.G = constants["universal_gravitational_constant"]

        # Function to check for an internet connection.
        def is_connected():
            try:
                host = socket.gethostbyname("www.github.com")
                s = socket.create_connection((host, 80), 2)
                s.close()
                return True
            except OSError:
                pass
            return False

        # Check for an internet connection.
        if not input_data:
            if not is_connected():
                raise Exception(
                    "You must connect to the internet in order to utilise AMSIMP."
                    + " Apologies for any inconvenience caused."
                )

        # RNN.
        self.epochs = epochs
        self.data_size = data_size

        # Ensure epochs is an integer value.
        if not isinstance(self.epochs, int):
            raise Exception(
                "epochs must be a integer value. The value of epochs was: {}".format(
                    self.ai
                )
            )

        # Ensure epochs is a natural number.
        if not self.epochs > 0:
            raise Exception(
                "epochs must be a integer value. The value of epochs was: {}".format(
                    self.ai
                )
            )

        # Ensure data_size is an integer value.
        if not isinstance(self.data_size, int):
            raise Exception(
                "data_size must be a integer value. The value of data_size was: {}".format(
                    self.ai
                )
            )

        # Ensure data_size is a natural number and is greater than 14.
        if not self.data_size > 14:
            raise Exception(
                "data_size must be a integer value. The value of data_size was: {}".format(
                    self.ai
                )
            )

# -----------------------------------------------------------------------------------------#

    cpdef np.ndarray latitude_lines(self):
        """
        Generates a NumPy array of latitude lines.
        """
        cdef float i
        if not self.input_data:
            latitude_lines = [
                i
                for i in np.arange(
                    -89, 90, self.delta_latitude
                )
            ]

            # Convert list to NumPy array and add the unit of measurement.
            latitude_lines = np.asarray(latitude_lines) * units.deg
        else:
            latitude_lines = self.lat

        return latitude_lines
    
    cpdef np.ndarray longitude_lines(self):
        """
        Generates a NumPy array of longitude lines.
        """
        cdef float i
        if not self.input_data:
            longitude_lines = [
                i
                for i in np.arange(
                    0, 360, self.delta_longitude
                )
            ]

            # Convert list to NumPy array and add the unit of measurement.
            longitude_lines = np.asarray(longitude_lines) * units.deg
        else:
            longitude_lines = self.lon

        return longitude_lines

    cpdef np.ndarray pressure_surfaces(self, dim_3d=False):
        """
        Generates a NumPy array of the constant pressure surfaces. This
        is the isobaric coordinate system.
        """
        if not self.input_data:
            pressure_surfaces = self.input_rh.coords('pressure')[0].points
            pressure_surfaces = pressure_surfaces[::-1] * self.units.Pa
            pressure_surfaces = pressure_surfaces.to(self.units.hPa)
            pressure_surfaces = pressure_surfaces.value
        else:
            pressure_surfaces = self.psurfaces.value

        # Convert Pressure Array into 3D Array.
        if dim_3d:
            list_pressure = []
            for p in pressure_surfaces:
                p = np.zeros((
                    len(self.latitude_lines()), len(self.longitude_lines())
                )) + p
                
                p = p.tolist()
                list_pressure.append(p)
            pressure_surfaces = np.asarray(list_pressure)

        pressure_surfaces *= self.units.hPa
        return pressure_surfaces

# -----------------------------------------------------------------------------------------#

    cpdef np.ndarray gradient_longitude(self, parameter=None):
        """
        Explain here.
        """
        # Determine gradient with respect to longitude.
        cdef np.ndarray longitude = np.radians(self.longitude_lines().value)

        # Define variable types.
        cdef np.ndarray parameter_central, parameter_boundaries
        cdef np.ndarray parameter_lgmt, parameter_rgmt
        cdef np.ndarray lon_lgmt, lon_rgmt, lon_boundaries

        # Compute gradient.
        # Central points.
        parameter_central = np.gradient(
            parameter, longitude, axis=2
        )

        # Boundaries.
        # Define atmospheric parameter boundaries.
        parameter_lgmt = parameter[:, :, -3:].value
        parameter_rgmt = parameter[:, :, :3].value
        # Define longitude boundaries.
        lon_lgmt = longitude[-3:]
        lon_rgmt = longitude[:3]
        # Concatenate atmospheric parameter boundaries.
        parameter_boundaries = np.concatenate(
            (parameter_lgmt, parameter_rgmt), axis=2
        ) * parameter.unit
        # Concatenate longitude boundaries.
        lon_boundaries = np.concatenate((lon_lgmt, lon_rgmt))
        # Compute gradient at boundaries.
        parameter_boundaries = np.gradient(
            parameter_boundaries, lon_boundaries, axis=2
        )
        # Redefine parameter as computed gradient.
        parameter = parameter_central
        parameter[:, :, -1] = parameter_boundaries[:, :, 2]
        parameter[:, :, 0] = parameter_boundaries[:, :, 3]

        # Make 1d latitude array into a 3d latitude array.
        cdef np.ndarray latitude = np.radians(self.latitude_lines())
        latitude = self.make_3dimensional_array(
            parameter=latitude, axis=1
        )

        # Handle longitudinal variation with respect to latitude.
        gradient_longitude = (
            (1 / (self.a * np.cos(latitude.value))) * parameter
        )

        return gradient_longitude

    cpdef np.ndarray gradient_latitude(self, parameter=None):
        """
        Explain here.
        """
        # Define variable.
        cdef np.ndarray latitude = np.radians(self.latitude_lines().value)

        # Determine gradient with respect to latitude.
        gradient_latitude = np.gradient(
            parameter.value, self.a.value * latitude, axis=1
        ) * (parameter.unit / units.m)

        return gradient_latitude
    
    cpdef np.ndarray gradient_pressure(self, parameter=None):
        """
        Explain here.
        """
        # Define variable.
        cdef np.ndarray pressure = self.pressure_surfaces()

         # Determine gradient with respect to pressure.
        gradient_pressure = np.gradient(
            parameter.value, pressure.value, axis=0
        ) * (parameter.unit / units.hPa)

        return gradient_pressure

    cpdef np.ndarray make_3dimensional_array(self, parameter=None, axis=1):
        """
        Explain here.
        """
        if axis == 0:
            a = self.latitude_lines().value
            b = self.longitude_lines().value
        elif axis == 1:
            a = self.pressure_surfaces().value
            b = self.longitude_lines().value
        elif axis == 2:
            a = self.pressure_surfaces().value
            b = self.latitude_lines().value

        list_parameter = []
        for param in parameter:
            unit = param.unit
            param = param.value + np.zeros((len(a), len(b)))
            param = param.tolist()
            list_parameter.append(param)
        parameter = np.asarray(list_parameter) * unit

        if axis == 1:
            parameter = np.transpose(parameter, (1, 0, 2))
        elif axis == 2:
            parameter = np.transpose(parameter, (1, 2, 0))
        
        return parameter

    cpdef dict retrieve_constants(self):
        """
        Explain here.
        """
        return self.constants

# -----------------------------------------------------------------------------------------

    cpdef np.ndarray coriolis_parameter(self):
        """
        Generates a NumPy arrray of the Coriolis parameter at various latitudes
        of the Earth's surface.
        
        The Coriolis parameter is defined as two times the angular rotation of
        the Earth by the sin of the latitude you are interested in.

        Equation:
            f_0 = 2 \* Upomega * sin(\phi)
        """
        coriolis_parameter = (
            2 * self.Upomega * np.sin(np.radians(self.latitude_lines()))
        )

        return coriolis_parameter

# -----------------------------------------------------------------------------------------#

    cpdef np.ndarray geopotential_height(self):
        """
        This method imports geopotential height data from a NetCDF file, which is
        located in the AMSIMP Initial Conditions Data Repository on GitHub.
        The data stored within this repo is updated every six hours by amsimp-bot.
        Following which, it outputs a NumPy array in the shape of
        (len(pressure_surfaces), len(latitude_lines), len(longitude_lines)).

        Geopotential Height is the height above sea level of a pressure level.
        For example, if a station reports that the 500 hPa height at its
        location is 5600 m, it means that the level of the atmosphere over
        that station at which the atmospheric pressure is 500 hPa is 5600
        meters above sea level.

        I strongly recommend storing the output into a variable in order to
        prevent the need to repeatly self.download the file. For more information,
        visit https://github.com/amsimp/initial-conditions.
        """
        if not self.input_data:
            # Input data.
            height = self.input_height

            pressure = self.pressure_surfaces().to(self.units.Pa)
            # Grid.
            grid_points = [
                ('pressure',  pressure.value),
                ('latitude',  self.latitude_lines().value),
                ('longitude', self.longitude_lines().value),                
            ]

            # Interpolation
            height = height.interpolate(
                grid_points, iris.analysis.Linear()
            )

            # Get data.
            geopotential_height = height.data
            geopotential_height = np.asarray(geopotential_height.tolist())

            geopotential_height *= self.units.m
        else:
            geopotential_height = self.input_height
        
        return geopotential_height

    cpdef np.ndarray relative_humidity(self):
        """
        This method imports relative humidity data from a NetCDF file, which is
        located in the AMSIMP Initial Conditions Data Repository on GitHub.
        The data stored within this repo is updated every six hours by amsimp-bot.
        Following which, it outputs a NumPy array in the shape of
        (len(pressure_surfaces), len(latitude_lines), len(longitude_lines)).
        
        Relative Humidity is the amount of water vapour present in air expressed
        as a percentage of the amount needed for saturation at the same temperature.

        I strongly recommend storing the output into a variable in order to
        prevent the need to repeatly self.download the file. For more information,
        visit https://github.com/amsimp/initial-conditions.
        """
        if not self.input_data:
            # Input data.
            rh = self.input_rh

            pressure = self.pressure_surfaces().to(self.units.Pa)
            # Grid.
            grid_points = [
                ('pressure',  pressure.value),
                ('latitude',  self.latitude_lines().value),
                ('longitude', self.longitude_lines().value),                
            ]

            # Interpolation
            rh = rh.interpolate(
                grid_points, iris.analysis.Linear()
            )

            # Get data.
            relative_humidity = rh.data
            relative_humidity = np.asarray(relative_humidity.tolist())

            relative_humidity *= self.units.percent
        else:
            relative_humidity = self.input_rh
        
        return relative_humidity

    cpdef np.ndarray temperature(self):
        """
        This method imports temperature data from a NetCDF file, which is
        located in the AMSIMP Initial Conditions Data Repository on GitHub.
        The data stored within this repo is updated every six hours by amsimp-bot.
        Following which, it outputs a NumPy array in the shape of
        (len(pressure_surfaces), len(latitude_lines), len(longitude_lines)).

        Temperature is defined as the mean kinetic energy density of molecular
        motion.

        I strongly recommend storing the output into a variable in order to
        prevent the need to repeatly self.download the file. For more information,
        visit https://github.com/amsimp/initial-conditions.
        """
        if not self.input_data:
            # Input data.
            temp = self.input_temp

            pressure = self.pressure_surfaces().to(self.units.Pa)
            # Grid.
            grid_points = [
                ('pressure',  pressure.value),
                ('latitude',  self.latitude_lines().value),
                ('longitude', self.longitude_lines().value),                
            ]

            # Interpolation
            temp = temp.interpolate(
                grid_points, iris.analysis.Linear()
            )

            # Get data.
            temperature = temp.data
            temperature = np.asarray(temperature.tolist())

            temperature *= self.units.K
        else:
            temperature = self.input_temp
        
        return temperature

    cpdef exit(self):
        """
        This method deletes any file created by the software,
        whilst downloading data from the AMSIMP Initial Atmospheric
        Conditions Data Repository. It does not, however, delete any
        forecast generated and saved by the simualte method.

        This command will also cause the Python interpreter to close
        and exit.
        """
        # Initial atmospheric conditions file.
        try:
            os.remove("initial_conditions.nc")
        except OSError:
            pass
        
        # Close Python.
        sys.exit()

# -----------------------------------------------------------------------------------------#

    cpdef np.ndarray pressure_thickness(self, p1=1000, p2=500):
        """
        Generates a NumPy array of pressure thickness between two defined
        constant pressure surfaces.

        Pressure thickness is defined as the distance between two
        pressure surfaces. Pressure thickness is determined through the
        hypsometric equation.

        Equation:
            h = z_2 - z_1
        """
        # Ensure p1 is greater than p2.
        if p1 < p2:
            raise Exception("Please note that p1 must be greater than p2.")

        # Find the nearest constant pressure surface to p1 and p2 in pressure.
        cdef int indx_p1 = (np.abs(self.pressure_surfaces().value - p1)).argmin()
        cdef int indx_p2 = (np.abs(self.pressure_surfaces().value - p2)).argmin()

        pressure_thickness = (
            self.geopotential_height()[indx_p2] - self.geopotential_height()[indx_p1]
        )

        return pressure_thickness

    cpdef np.ndarray troposphere_boundaryline(self):
        """
        Generates a NumPy array of the troposphere - stratosphere
        boundary line in the shape (len(longitude_lines), len(latitude_lines).
        This is calculated by looking at the vertical temperature profile in
        the method, amsimp.Backend.temperature().

        Explain here.
        """
        cdef np.ndarray temperature = self.temperature().value
        grad_temperature = np.gradient(temperature)
        cdef np.ndarray gradient_temperature = grad_temperature[0]
        gradient_temperature = np.transpose(gradient_temperature, (2, 1, 0))

        cdef np.ndarray psurface = self.pressure_surfaces().value

        cdef list trop_strat_line = []
        cdef list lat_trop_strat_line
        cdef np.ndarray temp, t
        cdef int n
        cdef float t_n, p
        for temp in gradient_temperature:
            lat_trop_strat_line = []
            for t in temp:
                n = 0
                while n < len(t):
                    t_n = t[n]
                    if t[n] >= 0:
                        p = psurface[n]
                        if p < 400:
                            lat_trop_strat_line.append(p)
                            n = len(t)
                    n += 1
            trop_strat_line.append(lat_trop_strat_line)

        troposphere_boundaryline = np.asarray(trop_strat_line)
        troposphere_boundaryline = np.transpose(troposphere_boundaryline, (1, 0)) 
        troposphere_boundaryline *= units.hPa
        return troposphere_boundaryline

# -----------------------------------------------------------------------------------------#

    def longitude_contourf(self, which="air_temperature", psurface=1000):
        """
        Plots a desired atmospheric process on a contour plot, with the axes
        being latitude and longitude. This plot is then layed on top of a
        EckertIII global projection. 
        
        For the raw data, please see the other methods found in this class.

        If you would like a particular atmospheric process to
        be added to this method, either create an issue on
        GitHub or send an email to support@amsimp.com.
        """
        if self.planet == "Earth":        
            # Ensure, which, is a string.
            if not isinstance(which, str):
                raise Exception(
                    "which must be a string of the name of the atmospheric parameter of interest."
                )

            # Index of the nearest pressure surface in amsimp.Backend.pressure_surfaces()
            indx_psurface = (np.abs(self.pressure_surfaces().value - psurface)).argmin()
            
            # Defines the axes, and the data.
            latitude, longitude = np.meshgrid(
                self.latitude_lines().value,
                self.longitude_lines().value
            )

            if which == "temperature" or which == "air_temperature":
                data = self.temperature()[indx_psurface, :, :]
                data_type = "Air Temperature"
                unit = " (K)"
            elif which == "geopotential_height" or which == "height":
                data = self.geopotential_height()[indx_psurface, :, :]
                data_type = "Geopotential Height"
                unit = " (m)"
            elif which == "density" or which == "atmospheric_density":
                data = self.density()[indx_psurface, :, :]
                data_type = "Atmospheric Density"
                unit = " ($\\frac{kg}{m^3}$)"
            elif which == "humidity" or which == "relative_humidity":
                data = self.relative_humidity()[indx_psurface, :, :]
                data_type = "Relative Humidity"
                unit = " (%)"
            elif which == "virtual_temperature":
                data = self.virtual_temperature()[indx_psurface, :, :]
                data_type = "Virtual Temperature"
                unit = " (K)"
            elif which == "vapor_pressure":
                data = self.vapor_pressure()[indx_psurface, :, :]
                data_type = "Vapor Pressure"
                unit = " (hPa)"
            elif which == "potential_temperature":
                data = self.potential_temperature()[indx_psurface, :, :]
                data_type = "Potential Temperature"
                unit = " (K)"
            elif which == "precipitable_water" or which == "precipitable_water_vapor":
                data = self.precipitable_water()
                data_type = "Precipitable Water Vapor"
                unit = " (mm)"
            elif which == "zonal_wind":
                data = self.wind()[0][indx_psurface, :, :]
                data_type = "Zonal Wind"
                unit = " ($\\frac{m}{s}$)"
            elif which == "meridional_wind":
                data = self.wind()[1][indx_psurface, :, :]
                data_type = "Meridional Wind"
                unit = " ($\\frac{m}{s}$)"
            elif which == "wind":
                data = self.wind()
                data = np.sqrt(data[0]**2 + data[1]**2)[indx_psurface, :, :]
                data_type = "Wind"
                unit = " ($\\frac{m}{s}$)"
            elif which == "mixing_ratio":
                data = self.mixing_ratio()[indx_psurface, :, :]
                data_type = "Mixing Ratio"
                unit = " (Dimensionless)"
            else:
                raise Exception(
                    "Invalid keyword. which must be a string of an atmospheric parameter included with AMSIMP."
                )

            if psurface < self.pressure_surfaces()[-1].value or psurface > self.pressure_surfaces()[0].value:
                raise Exception(
                    "psurface must be a real number within the isobaric boundaries. The value of psurface was: {}".format(
                        psurface
                    )
                )

            # EckertIII projection details.
            ax = plt.axes(projection=ccrs.EckertIII())
            ax.set_global()
            ax.coastlines()
            ax.gridlines()

            # Contourf plotting.
            minimum = data.min()
            maximum = data.max()
            levels = np.linspace(minimum, maximum, 21)
            data, lon = add_cyclic_point(data, coord=self.longitude_lines().value)
            data, lat = add_cyclic_point(
                np.transpose(data), coord=self.latitude_lines().value
            )
            data = np.transpose(data)
            contour = plt.contourf(
                lon,
                lat,
                data,
                transform=ccrs.PlateCarree(),
                levels=levels,
            )

            # Add SALT.
            plt.xlabel("Latitude ($\phi$)")
            plt.ylabel("Longitude ($\lambda$)")
            if not self.input_data:
                title = (
                    data_type + " ("
                    + str(self.date.year) + '-' + str(self.date.month) + '-'
                    + str(self.date.day) + " " + str(self.date.hour)
                    + ":00 h)"
                )
            else:
                title = data_type
        
            if which != "surface_temperature" and which != "precipitable_water" and which != "precipitable_water_vapor":
                title = (
                    title + " (" + "Pressure Surface = " + str(self.pressure_surfaces()[indx_psurface]) +")"
                )

            plt.title(title)

            # Colorbar creation.
            colorbar = plt.colorbar()
            tick_locator = ticker.MaxNLocator(nbins=15)
            colorbar.locator = tick_locator
            colorbar.update_ticks()
            colorbar.set_label(
                data_type + unit
            )

            plt.show()
            plt.close()
        else:
            raise NotImplementedError(
                "Visualisations for planetary bodies other than Earth is not currently implemented."
            )
    
    def psurface_contourf(self, which="air_temperature", central_long=352.3079):
        """
        Plots a desired atmospheric process on a contour plot,
        with the axes being latitude, and pressure surfaces. 
        
        For the raw data, please see the other methods found in
        this class.
        
        If you would like a particular atmospheric process to
        be added to this method, either create an issue on
        GitHub or send an email to support@amsimp.com.
        """
        # Ensure, which, is a string.
        if not isinstance(which, str):
            raise Exception(
                "which must be a string of the name of the atmospheric parameter of interest."
            )
        
        # Ensure central_long is between 0 and 359.
        if central_long < 0 or central_long > 359:
            raise Exception(
                "central_long must be a real number between 0 and 359. The value of central_long was: {}".format(
                    central_long
                )
            )
        
        # Index of the nearest central_long in amsimp.Backend.longtitude_lines()
        indx_long = (np.abs(self.longitude_lines().value - central_long)).argmin()

        # Defines the axes, and the data.
        latitude, pressure_surfaces = np.meshgrid(
            self.latitude_lines().value, self.pressure_surfaces()
        )
        
        if which == "temperature" or which == "air_temperature":
            data = self.temperature()[:, :, indx_long]
            data_type = "Air Temperature"
            unit = " (K)"
        elif which == "density" or which == "atmospheric_density":
            data = self.density()[:, :, indx_long]
            data_type = "Atmospheric Density"
            unit = " ($\\frac{kg}{m^3}$)"
        elif which == "humidity" or which == "relative_humidity":
            data = self.relative_humidity()[:, :, indx_long]
            data_type = "Relative Humidity"
            unit = " (%)"
        elif which == "virtual_temperature":
            data = self.virtual_temperature()[:, :, indx_long]
            data_type = "Virtual Temperature"
            unit = " (K)"
        elif which == "vapor_pressure":
            data = self.vapor_pressure()[:, :, indx_long]
            data_type = "Vapor Pressure"
            unit = " (hPa)"
        elif which == "potential_temperature":
            data = self.potential_temperature()[:, :, indx_long]
            data_type = "Potential Temperature"
            unit = " (K)"
        elif which == "zonal_wind":
            data = self.wind()[0][:, :, indx_long]
            data_type = "Zonal Wind"
            unit = " ($\\frac{m}{s}$)"
        elif which == "meridional_wind":
            data = self.wind()[1][:, :, indx_long]
            data_type = "Meridional Wind"
            unit = " ($\\frac{m}{s}$)"
        elif which == "mixing_ratio":
            data = self.mixing_ratio()[:, :, indx_long]
            data_type = "Mixing Ratio"
            unit = " (Dimensionless)"
        else:
            raise Exception(
                "Invalid keyword. which must be a string of an atmospheric parameter included with AMSIMP."
            )

        # Contourf plotting.
        minimum = data.min()
        maximum = data.max()
        levels = np.linspace(minimum, maximum, 21)
        plt.contourf(
            latitude,
            pressure_surfaces,
            data,
            levels=levels,
        )

        # Add SALT.
        plt.xlabel("Latitude ($\phi$)")
        plt.ylabel("Pressure (hPa)")
        plt.yscale('log')
        plt.gca().invert_yaxis()
        if not self.input_data:
            plt.title(
                data_type + " ("
                + str(self.date.year) + '-' + str(self.date.month) + '-'
                + str(self.date.day) + " " + str(self.date.hour)
                + ":00 h, Longitude = "
                + str(np.round(self.longitude_lines()[indx_long], 2)) + ")"
            )
        else:
            plt.title(
                data_type + " ("
                + "Longitude = "
                + str(np.round(self.longitude_lines()[indx_long], 2)) + ")"
            )

        # Colorbar creation.
        colorbar = plt.colorbar()
        tick_locator = ticker.MaxNLocator(nbins=15)
        colorbar.locator = tick_locator
        colorbar.update_ticks()
        colorbar.set_label(
            data_type + unit
        )

        # Average boundary line between the troposphere and the stratosphere.
        troposphere_boundaryline = self.troposphere_boundaryline()
        avg_tropstratline = np.mean(troposphere_boundaryline) + np.zeros(
            len(troposphere_boundaryline)
        )

        # Plot average boundary line on the contour plot.
        plt.plot(
            latitude[1],
            avg_tropstratline,
            color="black",
            linestyle="dashed",
            label="Troposphere - Stratosphere Boundary Line",
        )
        plt.legend(loc=0)

        plt.show()
        plt.close()
    
    def thickness_contourf(self, p1=1000, p2=500):
        """
        Plots pressure thickness on a contour plot, with the axes being
        latitude and longitude. This plot is then layed on top of a EckertIII
        global projection. 
        
        For the raw data, please use the amsimp.Backend.pressure_thickness()
        method.
        """
        if self.planet == "Earth":
            # Defines the axes, and the data.
            latitude, longitude = self.latitude_lines().value, self.longitude_lines().value

            pressure_thickness = self.pressure_thickness(p1=p1, p2=p2)

            # EckertIII projection details.
            ax = plt.axes(projection=ccrs.EckertIII())
            ax.set_global()
            ax.coastlines()
            ax.gridlines()

            # Contourf plotting.
            minimum = pressure_thickness.min()
            maximum = pressure_thickness.max()
            levels = np.linspace(minimum, maximum, 21)
            pressure_thickness, longitude = add_cyclic_point(
                pressure_thickness, coord=longitude
            )
            pressure_thickness, latitude = add_cyclic_point(
                np.transpose(pressure_thickness), coord=latitude
            )
            pressure_thickness = np.transpose(pressure_thickness)
            contour = plt.contourf(
                longitude,
                latitude,
                pressure_thickness,
                transform=ccrs.PlateCarree(),
                levels=levels,
            )

            # Index of the rain / snow line
            indx_snowline = (np.abs(levels.value - 5400)).argmin()
            contour.collections[indx_snowline].set_color('black')
            contour.collections[indx_snowline].set_linewidth(1) 

            # Add SALT.
            plt.xlabel("Latitude ($\phi$)")
            plt.ylabel("Longitude ($\lambda$)")
            if not self.input_data:
                plt.title("Pressure Thickness ("
                    + str(self.date.year) + '-' + str(self.date.month) + '-'
                    + str(self.date.day) + " " + str(self.date.hour) + ":00 h" + ")"
                )
            else:
                plt.title("Pressure Thickness")

            # Colorbar creation.
            colorbar = plt.colorbar()
            tick_locator = ticker.MaxNLocator(nbins=15)
            colorbar.locator = tick_locator
            colorbar.update_ticks()
            colorbar.set_label(
                "Pressure Thickness (" + str(p1) + " hPa - " + str(p2) + " hPa) (m)"
            )

            # Footnote
            if p1 == 1000 and p2 == 500:
                plt.figtext(
                    0.99,
                    0.01,
                    "Rain / Snow Line is marked by the black line (5,400 m).",
                    horizontalalignment="right",
                )

            plt.show()
            plt.close()
        else:
            raise NotImplementedError(
                "Visualisations for planetary bodies other than Earth is not currently implemented."
            )

Wind Class

  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
#cython: linetrace=True
#distutils: define_macros=CYTHON_TRACE_NOGIL=1
#cython: language_level=3
"""
AMSIMP Wind Class. For information about this class is described below.

Copyright (C) 2020 AMSIMP

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program 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 General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <https://www.gnu.org/licenses/>.
"""

# -----------------------------------------------------------------------------------------#

# Importing Dependencies
import os
import iris
import cartopy.crs as ccrs
from cartopy.util import add_cyclic_point
import matplotlib.pyplot as plt
from matplotlib import ticker
import numpy as np
from amsimp.moist cimport Moist
from amsimp.moist import Moist
cimport numpy as np

# -----------------------------------------------------------------------------------------#

cdef class Wind(Moist):
    """
    AMSIMP Wind Class - This class is concerned with calculating numerical
    values for wind in the troposphere and the stratosphere.

    Below is a list of the methods included within this class, with a short
    description of their intended purpose. Please see the relevant class methods
    for more information.

    wind ~ generates a NumPy array of the components of wind. The unit of
    measurement is metres per second.
    static_stability ~ generates a NumPy array of the the static
    stability within a vertical profile.
    verical_motion ~ generates a NumPy array of the vertical motion
    with the atmosphere, by integrating the isobaric version of the
    mass continunity equation.
    wind_contourf ~ generates a vertical profile of geostrophic wind,
    and outputs this as a contour plot.

    globe ~ generates a geostrophic wind contour plot, adds wind vectors to
    that said plot, and overlays both on a Nearside Projection of the Earth.
    """

    cpdef np.ndarray static_stability(self):
        """
        Generates a NumPy array of the the static stability within
        a vertical profile. 
        
        Static stability is defined as the stability of the
        atmosphere in hydrostatic equilibrium with respect to vertical
        displacements.

        Equation:
            sigma = - frac{R T}{p theta} frac{partial theta}{partial p}
        """
        cdef np.ndarray theta = self.potential_temperature(moist=False)
        cdef np.ndarray temperature = self.temperature()
        cdef np.ndarray pressure = self.pressure_surfaces(dim_3d=True)

        static_stability = - self.R * (
                temperature / (pressure * theta)
            ) * self.gradient_p(
                parameter=theta
        )
        return static_stability

    cpdef tuple wind(self):
        """
        This method imports wind data from a GRIB file, which is
        located in the AMSIMP Initial Conditions Data Repository on GitHub.
        The data stored within this repo is updated every six hours by amsimp-bot.
        Following which, it outputs a NumPy array in the shape of
        (len(pressure_surfaces), len(latitude_lines), len(longitude_lines)).
        
        Wind is the flow of gases on a large scale. Wind consists of
        the bulk movement of air.
        
        I strongly recommend storing the output into a variable in order to
        prevent the need to repeatly download the file. For more information,
        visit https://github.com/amsimp/initial-conditions.

        For geostrophic wind, please see the method, wind.
        """
        if not self.input_data:
            # Input data.
            # Zonal.
            u = self.input_u
            # Meridional.
            v = self.input_v

            pressure = self.pressure_surfaces().to(self.units.Pa)
            # Grid.
            grid_points = [
                ('pressure',  pressure.value),
                ('latitude',  self.latitude_lines().value),
                ('longitude', self.longitude_lines().value),                
            ]

            # Interpolation
            # Zonal.
            u = u.interpolate(
                grid_points, iris.analysis.Linear()
            )
            # Meridional.
            v = v.interpolate(
                grid_points, iris.analysis.Linear()
            )

            # Get data.
            # Zonal.
            u = u.data
            u = np.asarray(u.tolist())
            # Meridional.
            v = v.data
            v = np.asarray(v.tolist())

            u *= self.units.m / self.units.s
            v *= self.units.m / self.units.s
        else:
            # Zonal.
            u = self.input_u
            # Meridional.
            v = self.input_v

        return u, v

    cpdef np.ndarray vertical_motion(self):
        """
        Generates a NumPy array of the vertical motion within
        the atmosphere, by integrating the isobaric version of
        the mass continunity equation. 
        
        Typical large-scale vertical motions in the atmosphere
        are of the order of 0.01 – 0.1 m s−1.

        Equation:
            frac{\partial u}{\partial x} + frac{\partial v}{\partial y} = - frac{\partial \omega}{\partial p}
        """
        # Define variables.
        cdef np.ndarray latitude = np.radians(self.latitude_lines()).value
        cdef np.ndarray longitude = np.radians(self.longitude_lines()).value
        cdef np.ndarray pressure = self.pressure_surfaces(dim_3d=True)
        cdef np.ndarray u, v
        u, v = self.wind()
        
        # Determine the LHS pf the equation by calculating the derivatives.
        # Change in meridional wind with respect to latitude.
        v_dy = self.gradient_latitude(parameter=v)

        # Change in zonal wind with respect to longitude.
        u_dx = self.gradient_longitude(parameter=u)

        LHS = u_dx + v_dy
        LHS *= -1

        # Integrate the continunity equation.
        cdef list vertical_motion_list = []
        cdef int n = 0
        cdef int len_pressure = len(pressure)
        cdef omega, omega_unit
        while n < len_pressure:
            p1 = n
            p2 = n+2
            y = LHS[p1:p2, :, :]
            p = pressure[p1:p2, :, :]
            omega = np.trapz(
                y=y, x=p, axis=0
            )
            omega_unit = omega.unit
            vertical_motion_list.append(omega.value)

            n += 1

        # Convert list to NumPy array.
        vertical_motion = np.asarray(vertical_motion_list)

        # Add units.
        vertical_motion *= omega_unit

        # Ensure the shape is correct.
        if np.shape(vertical_motion) != np.shape(u):
            raise Exception("Unable to determine vertical motion"
                + " at this time. Please contact the developer for futher"
                + " assistance.")

        return vertical_motion

# -----------------------------------------------------------------------------------------#

    def globe(self, central_lat=53.1424, central_long=352.3079, psurface=1000, which_wind=0):
        """
        This particular method adds wind vectors to a contour plot. It also
        overlays both of these elements onto a Nearside Perspective projection
        of the Earth.

        By default, the perspective view is looking directly down at the city
        of Dublin in the country of Ireland. If which_wind is set to 0, it will
        plot the non-geostophic wind data, otherwise, if it is set to 1, it will
        plot the geostrophic wind data.

        Note:
        The NumPy method, seterr, is used to suppress a weird RunTime warning
        error.
        """
        if self.planet == "Earth":
            # Ensure central_lat is between -90 and 90.
            if central_lat < -90 or central_lat > 90:
                raise Exception(
                    "central_lat must be a real number between -90 and 90. The value of central_lat was: {}".format(
                        central_lat
                    )
                )

            # Ensure central_long is between 0 and 359.
            if central_long < 0 or central_long > 359:
                raise Exception(
                    "central_long must be a real number between 0 and 359. The value of central_long was: {}".format(
                        central_long
                    )
                )
            
            # Ensure psurface is between 1000 and 1 hPa above sea level.
            if psurface < 1 or psurface > 1000:
                raise Exception(
                    "psurface must be a real number between 1 and 1,000. The value of psurface was: {}".format(
                        psurface
                    )
                )
            
            # Index of the nearest alt in amsimp.Backend.altitude_level()
            indx_psurface = (np.abs(self.pressure_surfaces().value - psurface)).argmin()

            # Ignore NumPy errors.
            np.seterr(all="ignore")

            # Define the axes, and the data.
            latitude = self.latitude_lines().value
            longitude = self.longitude_lines().value

            wind = self.wind()
            u = wind[0][indx_psurface, :, :]
            v = wind[1][indx_psurface, :, :]
            title = "Wind"

            u_norm = u / np.sqrt(u ** 2 + v ** 2)
            v_norm = v / np.sqrt(u ** 2 + v ** 2)

            wind = np.sqrt(u ** 2 + v ** 2)

            ax = plt.axes(
                projection=ccrs.NearsidePerspective(
                    central_longitude=central_long, central_latitude=central_lat
                )
            )

            # Add latitudinal and longitudinal grid lines, as well as, 
            # coastlines to the globe.
            ax.set_global()
            ax.coastlines()
            ax.gridlines()

            # Contour plotting.
            minimum = wind.min()
            maximum = wind.max()
            levels = np.linspace(minimum, maximum, 21)
            wind, lon = add_cyclic_point(wind, coord=longitude)
            wind, lat = add_cyclic_point(np.transpose(wind), coord=latitude)
            wind = np.transpose(wind)
            contourf = plt.contourf(
                lon,
                lat,
                wind,
                transform=ccrs.PlateCarree(),
                levels=levels,
            )
            
            # Reduce density of wind vectors.
            skip_lat = latitude.shape[0] / 23
            skip_lat = int(np.round(skip_lat))
            skip_lon = longitude.shape[0] / 23
            skip_lon = int(np.round(skip_lat))

            # Add geostrophic wind vectors.
            plt.quiver(
                longitude[::skip_lon],
                latitude[::skip_lat],
                u_norm.value[::skip_lat, ::skip_lon],
                v_norm.value[::skip_lat, ::skip_lon],
                transform=ccrs.PlateCarree(),   
            )

            # Add SALT.
            plt.xlabel("Latitude ($\phi$)")
            plt.ylabel("Longitude ($\lambda$)")
            if not self.input_data:
                plt.title(title + " ("
                    + str(self.date.year) + '-' + str(self.date.month) + '-'
                    + str(self.date.day) + " " + str(self.date.hour)
                    + ":00 h, Pressure Surface = "
                    + str(self.pressure_surfaces()[indx_psurface]) +")"
                )
            else:
                plt.title(title + " ("
                    + "Pressure Surface = "
                    + str(self.pressure_surfaces()[indx_psurface]) +")"
                )

            # Add colorbar.
            colorbar = plt.colorbar(contourf)
            tick_locator = ticker.MaxNLocator(nbins=15)
            colorbar.locator = tick_locator
            colorbar.update_ticks()
            colorbar.set_label("Velocity ($\\frac{m}{s}$)")

            # Footnote
            if which_wind == 1:
                plt.figtext(
                    0.99,
                    0.01,
                    "Note: Geostrophic balance does not hold near the equator.",
                    horizontalalignment="right",
                )

            plt.show()
            plt.close()
        else:
            raise NotImplementedError(
                "Visualisations for planetary bodies other than Earth is not currently implemented."
            )

Moist Class

  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
#cython: linetrace=True
#distutils: define_macros=CYTHON_TRACE_NOGIL=1
#cython: language_level=3
"""
AMSIMP Moist Thermodynamics Class. For information about this class is described below.

Copyright (C) 2020 AMSIMP

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program 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 General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <https://www.gnu.org/licenses/>.
"""

# -----------------------------------------------------------------------------------------#

# Importing Dependencies
import numpy as np
from scipy.integrate import quad
import cartopy.crs as ccrs
from cartopy.util import add_cyclic_point
import matplotlib.pyplot as plt
from matplotlib import ticker
from amsimp.backend cimport Backend
from amsimp.backend import Backend
cimport numpy as np
from cpython cimport bool

# -----------------------------------------------------------------------------------------#


cdef class Moist(Backend):
    """
    AMSIMP Moist Thermodynamics Class - This class is concerned with incorpating
    moisture into the atmospheric model. This is done through the utilisation of
    Mosit Thermodynamics. Atmospheric thermodynamics is the study of heat-to-work
    transformations (and their reverse) that take place in the earth's atmosphere
    and manifest as weather or climate. Atmospheric thermodynamics use the laws
    of classical thermodynamics, to describe and explain such phenomena as
    the properties of moist air

    vapor_pressure ~ generates a NumPy array of saturated vapor pressure.
    virtual_temperature ~ generates a NumPy array of the virtual temperature.
    The unit of measurement is Kelvin.

    density ~ generates a NumPy array of the atmospheric density. The unit of
    measurement is kilograms per cubic metric.
    exner_function ~ generates a NumPy array of the Exner function (4). This
    method has no unit of measurement, i.e. it is dimensionless.
    mixing_ratio ~ generates a NumPy array of the mixing ratio (water vapor).
    This method has no unit of measurement, i.e. it is dimensionless.
    potential_temperature ~ outputs a NumPy array of potential temperature.
    The unit of measurement is Kelvin.
    
    precipitable_water ~ generates a NumPy array of precipitable water
    vapor. The unit of measurement is millimetres.
    """

    cpdef np.ndarray vapor_pressure(self):
        """
        Generates a NumPy array of vapor pressure. 
        
        Vapor pressure, in meteorology, is the partial pressure of water vapor.
        The partial pressure of water vapor is the pressure that water vapor
        contributes to the total atmospheric pressure.

        Equation (Saturated Vapor Pressure):
            e_s = 6.112 \* \exp((17.67 * T) / (T + 243.15))

        Equation (Vapor Pressure):
            e = frac{r_h \* e_s}{100}
        """
        # Convert temperature in Kelvin to degrees centigrade.
        cdef np.ndarray temperature = self.temperature().value - 273.15

        # Saturated water vapor pressure
        sat_vapor_pressure = 0.61121 * np.exp(
            (
                18.678 - (temperature / 234.5)
            ) * (temperature / (257.14 + temperature)
            )
        )

        # Add units of measurement.
        sat_vapor_pressure *= self.units.kPa
        sat_vapor_pressure = sat_vapor_pressure.to(self.units.hPa)

        # Vapor pressure, accounting for relative humidity.
        vapor_pressure = self.relative_humidity().value * sat_vapor_pressure
        vapor_pressure /= 100

        return vapor_pressure

    cpdef np.ndarray virtual_temperature(self):
        """
        Generates a NumPy array of the virtual temperature. 
        
        The virtual temperature is the temperature at which dry air
        would have the same density as the moist air, at a given
        pressure. In other words, two air samples with the same
        virtual temperature have the same density, regardless
        of their actual temperature or relative humidity.
        
        Equation:
        T_v = frac{T}{1 - frac{0.378 e}{p}}
        """
        virtual_temperature = self.temperature() / (
            1 - (
            self.vapor_pressure() / self.pressure_surfaces(dim_3d=True)
            ) * (1 - 0.622)
        )

        return virtual_temperature

# -----------------------------------------------------------------------------------------#

    cpdef np.ndarray density(self):
        """
        Generates a NumPy array of atmospheric pressure utilising the Ideal
        Gas Law.

        The atmospheric density is the mass of the atmosphere per unit
        volume. The ideal gas equation is the equation of state for the
        atmosphere, and is defined as an equation relating temperature,
        pressure, and specific volume of a system in theromodynamic
        equilibrium.

        Equation:
            rho = frac{p}{R * T_v}
        """
        cdef np.ndarray pressure_surfaces = self.pressure_surfaces(dim_3d=True)

        density = pressure_surfaces / (self.virtual_temperature() * self.R)

        # Switch to the appropriate SI units.
        density = density.si
        return density

    cpdef np.ndarray exner_function(self):
        """
        Generates a NumPy array of the exner function. 
        
        The Exner function can be viewed as non-dimensionalized pressure.

        Equation:
            \Pi = frac{T}{theta}
        """
        exner_function = self.virtual_temperature() / self.potential_temperature()

        return exner_function

    cpdef np.ndarray mixing_ratio(self):
        """"
        Generates a NumPy array of the mixing ratio.

        The mixing ratio is the ratio of the mass of a variable atmospheric
        constituent to the mass of dry air. In this particular case, it refers
        to water vapor.

        Equation:
            q = frac{0.622 * e}{p - e}
        """
        mixing_ratio = (
            0.622 * self.vapor_pressure()
        ) / (
                self.pressure_surfaces(dim_3d=True) - self.vapor_pressure()
        )

        return mixing_ratio

    cpdef np.ndarray potential_temperature(self, moist=False):
        """
        Generates a NumPy array of potential temperature. 
        
        The potential temperature of a parcel of fluid at pressure P is the
        temperature that the parcel would attain if adiabatically brought
        to a standard reference pressure.

        Equation:
            theta = T \* (frac{P}{P_0}) ** (-R / c_p)
        """
        # Ensure moist is a boolean value.
        if not isinstance(moist, bool):
            raise Exception(
                "moist must be a boolean value. The value of moist was: {}".format(
                    moist
                )
            )

        # Determine whether to calculate wet-bulb potential temperature, or
        # potential temperature.
        cdef np.ndarray temperature
        if moist:
            temperature = self.virtual_temperature().value
        else:
            temperature = self.temperature().value

        cdef np.ndarray pressure_surfaces = self.pressure_surfaces(dim_3d=True).value
        cdef float R = self.R.value
        cdef float c_p = self.c_p.value

        cdef list list_potentialtemperature = []
        cdef int n = 0
        cdef int len_psurfaces = len(pressure_surfaces)
        while n < len_psurfaces:
            theta = temperature[n] * (
                (pressure_surfaces[n] / pressure_surfaces[0]) ** (-R / c_p)
            )
            
            list_potentialtemperature.append(theta)
            
            n += 1

        potential_temperature = np.asarray(list_potentialtemperature)
        potential_temperature *= self.units.K
        return potential_temperature

# -----------------------------------------------------------------------------------------#

    def __mixing_ratio(self, pressure, vapor_pressure):
        """
        This method is solely utilised for integration in the
        amsimp.Water.precipitable_water() method.
        """
        y = (0.622 * vapor_pressure) / (pressure - vapor_pressure)

        return y

    cpdef np.ndarray precipitable_water(self, sum_pwv=True):
        """
        Generates a NumPy array of saturated precipitable water vapor.
        
        Precipitable water is the total atmospheric water vapor contained in a
        vertical column of unit cross-sectional area extending between any two
        specified levels. For a contour plot of this data, please use the
        amsimp.Water.contourf() method.

        Equation:
            PW = frac{1}{rho \* g} \* \int r dp
        """
        # Defining some variables.
        cdef np.ndarray pressure = self.pressure_surfaces(dim_3d=True).to(self.units.Pa).value
        pressure = np.transpose(pressure, (2, 1, 0))
        cdef np.ndarray vapor_pressure = self.vapor_pressure().to(self.units.Pa)
        vapor_pressure = np.transpose(vapor_pressure.value, (2, 1, 0))
        cdef g = self.g
        cdef rho_w = 997 * (self.units.kg / self.units.m ** 3)

        # Integrate the mixing ratio with respect to pressure between the
        # pressure boundaries of p1, and p2.
        cdef list list_precipitablewater = []
        cdef list pwv_lat, pwv_alt
        cdef tuple integration_term
        cdef float p1, p2, e, Pwv_intergrationterm
        cdef P_wv
        cdef int len_pressurelong = len(pressure)
        cdef int len_pressurelat = len(pressure[0])
        cdef int len_pressurealt = len(pressure[0][0]) - 1
        cdef int i, n, k
        i = 0
        while i < len_pressurelong:
            pwv_lat = []

            n = 0
            while n < len_pressurelat:
                pwv_alt = []

                k = 0
                while k < len_pressurealt:
                    p1 = pressure[i][n][k]
                    p2 = pressure[i][n][k + 1]
                    e = (vapor_pressure[i][n][k] + vapor_pressure[i][n][k + 1]) / 2

                    integration_term = quad(self.__mixing_ratio, p1, p2, args=(e,))
                    Pwv_intergrationterm = integration_term[0]
                    pwv_alt.append(Pwv_intergrationterm)

                    k += 1

                if sum_pwv: 
                    P_wv = np.sum(pwv_alt)
                else:
                    P_wv = pwv_alt
                
                pwv_lat.append(P_wv)                   

                n += 1
            
            list_precipitablewater.append(pwv_lat)
            i += 1

        precipitable_water = np.asarray(list_precipitablewater) * self.units.Pa

        if np.shape(precipitable_water) == (len_pressurelong, len_pressurelat):
            precipitable_water = np.transpose(precipitable_water)
        else:
            precipitable_water = np.transpose(precipitable_water, (2, 1, 0))

        # Multiplication term by integration.
        precipitable_water *= -1 / (rho_w * g)

        precipitable_water = precipitable_water.to(self.units.mm)
        return precipitable_water

Dynamics Class

   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
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
#cython: linetrace=True
#distutils: define_macros=CYTHON_TRACE_NOGIL=1
#cython: language_level=3
"""
AMSIMP Recurrent Neural Network and Dynamics Class. For information about
this class is described below.

Copyright (C) 2020 AMSIMP

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program 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 General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <https://www.gnu.org/licenses/>.
"""

# -----------------------------------------------------------------------------------------#

# Importing Dependencies
from datetime import timedelta, datetime
import os, sys
from sklearn.preprocessing import MinMaxScaler
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense, RepeatVector, TimeDistributed
from tensorflow.keras.layers import Bidirectional
from tensorflow.keras.optimizers import Adam
import tensorflow as tf
import matplotlib.pyplot as plt
from matplotlib import style, ticker, gridspec
import numpy as np
from astropy import constants as constant
import cartopy.crs as ccrs
from cartopy.util import add_cyclic_point
from cpython cimport bool
from amsimp.wind cimport Wind
from amsimp.wind import Wind
from astropy.units.quantity import Quantity
import iris
from iris.coords import DimCoord, AuxCoord
from iris.cube import Cube, CubeList
from iris import save
from progress.bar import IncrementalBar

# -----------------------------------------------------------------------------------------#

cdef class RNN(Wind):
    """
    AMSIMP Recurrent Neural Network Class - This class is concerned with
    determining the evolution of the state of the atmosphere through
    the utilisation of a recurrent neural network. Weather forecasting
    has traditionally been done by physical models of the atmosphere,
    which are unstable to perturbations, and thus are inaccurate for
    large periods of time. Since machine learning techniques are more
    robust to perturbations, it would be logical to combine a neural network
    with a physical model. Weather forecasting is a sequential data problem,
    therefore, a recurrent neural network is the most suitable option for
    this task. For this particular version of the software, the
    recurrent neural network will be trained solely on temperature and
    relative humidity data. Geopotential height will be added in a future release.
    A seed is set in order to ensure reproducibility.

    Below is a list of the methods included within this class, with a short
    description of their intended purpose. Please see the relevant class methods
    for more information:

    download_historical_data ~ generates a NumPy array of the previous state
    of the atmosphere. The number of days of data which is generated can
    be specified through the utilisation of the parameter, data_size, when
    the class is initialised.
    standardise_data ~ generates a NumPy array which standardises the data
    generated by the download_historical_data method. 
    preprocess_data ~ generates the training dataset for the recurrent
    neural network.
    model_prediction ~ generates a NumPy array of evolution of the
    future state of the atmosphere. This is the prediction generated
    by the recurrent neural network by training on the generated and
    preprocessed training data.
    """

    # Feature Scaling 
    temp_sc = MinMaxScaler(feature_range=(0,1))
    rh_sc = MinMaxScaler(feature_range=(0,1))
    height_sc = MinMaxScaler(feature_range=(0,1))

    # Ensure reproducibility.
    tf.random.set_seed(13)

    def download_historical_data(self):
        """
        Generates a NumPy array of the previous state of the atmosphere.
        The number of days of data which is generated can be specified
        through the utilisation of the parameter, data_size, when the
        class is initialised.

        This data is retrieved from the historical dataset contained
        on the AMSIMP Initial Atmospheric Conditions Repository.
        """
        # Folder containing historical data on GitHub.
        folder = "https://github.com/amsimp/initial-conditions/raw/master/initial_conditions/"

        # Data lists.
        # Temperature.
        cdef list T_list = []

        # Geopotential Height.
        cdef list geo_list = []

        # Relative Humidity.
        cdef list rh_list = []

        # Define variable types.
        cdef np.ndarray T, height, rh

        # Beginning date of the dataset.
        cdef date = self.date
        date = date + timedelta(days=-self.data_size)

        # Download progresss.
        max_bar = self.data_size * 4
        bar = IncrementalBar('Downloading Historical Data', max=max_bar)

        # Check if the parameter, input_data, is set to True.
        if self.input_data:
            raise NotImplementedError(
                "Currently, it is not possible to train the neural network on user defined datasets."
            ) 

        # Remove initial conditions file.
        try:
            os.remove("initial_conditions.nc")
        except OSError:
            pass

        for i in range(max_bar):
            # Configure the Wind class, so, that it aligns with the
            # paramaters defined by the user.
            config = Wind(
                delta_latitude=self.delta_latitude,
                delta_longitude=self.delta_longitude,
                input_date=date,
            )

            # Redefine NumPy arrays.
            # Geopotential Height.
            height = config.geopotential_height().value
            height = height.flatten()
            # Temperature.
            T = config.temperature().value
            T = T.flatten()
            # Relative Humidity.
            rh = config.relative_humidity().value
            rh = rh.flatten()

            # Append to list.
            # Geopotential Height.
            geo_list.append(height)
            # Temperature.
            T_list.append(T)
            # Relative Humidity.
            rh_list.append(rh)

            # Add six hours to date to get the next dataset.
            date = date + timedelta(hours=+6)

            # Remove download file.
            os.remove("initial_conditions.nc")

            # Increment progress bar.
            bar.next()

        # Convert lists to NumPy arrays.
        # Geopotential Height.
        geopotential_height = np.asarray(geo_list)
        # Temperature.
        temperature = np.asarray(T_list)
        # Relative Humidity.
        relative_humidity = np.asarray(rh_list)

        # Finish bar.
        bar.finish()

        # Output.s
        output = (temperature, relative_humidity, geopotential_height)
        return output        

    def standardise_data(self):
        """
        Generates a NumPy array which standardises the data
        generated by the download_historical_data method. It is important
        to scale features before training a neural network.
        
        Normalisation is a common way of doing this scaling by
        subtracting the mean and dividing by the standard deviation
        of each feature. In order for the most optimal performance,
        the method ``MinMaxScaler" from the library, scikit-learn,
        is utilised within the software.
        """
        # Define atmospheric parameters.
        historical_data = self.download_historical_data()
        temperature = historical_data[0]
        relative_humidity = historical_data[1]
        geopotential_height = historical_data[2]

        # Standardise the data.
        # Temperature.
        temperature = self.temp_sc.fit_transform(temperature)
        # Relative Humidity.
        relative_humidity = self.rh_sc.fit_transform(relative_humidity)
        # Geopotential Height.
        geopotential_height = self.height_sc.fit_transform(geopotential_height)

        # Output.
        output = (temperature, relative_humidity, geopotential_height)
        return output

    def preprocess_data(
        self, dataset, past_history, future_target
    ):
        """
        Generates the training dataset for the recurrent neural network.

        The data set in question is updated every six hours by the National Oceanic
        and Atmospheric Administration. This means for a single day,
        there will be four observations. The goal for the software will be to,
        first predict the relevant atmospheric parameter in seven days time
        given the last thirty days of data and combine this RNN prediction with the
        physical model prediction in an attempt to make a more accurate prediction
        overall. In order to make such predictions, it is necessary to create a
        window of the last 120 (30 x 4) observations to train the model.
        """
        X, y = list(), list()
        for i in range(len(dataset)):
            # Find the end.
            end_ix = i + past_history
            out_end_ix = end_ix + future_target
            
            # Determine if we are beyond the dataset.
            if out_end_ix > len(dataset):
                break
            
            # Gather the input and output components.
            seq_x, seq_y = dataset[i:end_ix, :], dataset[end_ix:out_end_ix, :]

            # Append to list.
            X.append(seq_x)
            y.append(seq_y)

        return np.asarray(X), np.asarray(y)

    def model_prediction(self):
        """
        Generates a NumPy array of evolution of the future state of the atmosphere.

        This is the prediction generated by the recurrent neural network by training
        on the generated and preprocessed training data. Please see the other
        methods within this class for preprocessing information.
        """
        # Dataset.
        dataset = self.standardise_data()

        # Temperature.
        temperature = dataset[0]
        # Relative Humidity.
        relative_humidity = dataset[1]
        # Geopotential Height.
        geopotential_height = dataset[2]

        # The network is shown data from the last 15 days.
        past_history = 20 * 4

        # The network predicts the next 7 days worth of steps.
        future_target = 7 * 4

        # The dataset is preprocessed.
        # Temperature.
        x_temp, y_temp = self.preprocess_data(
            temperature, past_history, future_target
        )
        # Relative Humidity.
        x_rh, y_rh = self.preprocess_data(
            relative_humidity, past_history, future_target
        )
        # Geopotential Height.
        x_height, y_height = self.preprocess_data(
            geopotential_height, past_history, future_target
        )

        # Number of features.
        features = x_temp.shape[2]

        # Create, and train models.
        # Temperature model.
        # Optimiser.
        opt_temp = Adam(lr=1e-6, decay=1e-10, clipvalue=0.6)
        # Create.
        temp_model = Sequential()
        temp_model.add(
            LSTM(
                400, activation='tanh', input_shape=(past_history, features)
            )
        )
        temp_model.add(RepeatVector(future_target))
        temp_model.add(LSTM(400, activation='tanh', return_sequences=True))
        temp_model.add(LSTM(400, activation='relu', return_sequences=True))
        temp_model.add(LSTM(400, activation='relu', return_sequences=True))
        temp_model.add(TimeDistributed(Dense(features)))
        temp_model.compile(optimizer=opt_temp, loss='mse', metrics=['mean_absolute_error'])
        # Train.
        temp_model.fit(
            x_temp, y_temp, epochs=self.epochs, batch_size=20
        )

        # Relative Humidity model.
        # Optimiser.
        opt_rh = Adam(lr=1e-6, decay=1e-10, clipvalue=0.6)
        # Create.
        rh_model = Sequential()
        rh_model.add(
            LSTM(
                400, activation='tanh', input_shape=(past_history, features)
            )
        )
        rh_model.add(RepeatVector(future_target))
        rh_model.add(LSTM(400, activation='tanh', return_sequences=True))
        rh_model.add(LSTM(400, activation='relu', return_sequences=True))
        rh_model.add(LSTM(400, activation='relu', return_sequences=True))
        rh_model.add(TimeDistributed(Dense(features)))
        rh_model.compile(optimizer=opt_rh, loss='mse', metrics=['mean_absolute_error'])
        # Train.
        rh_model.fit(
            x_rh, y_rh, epochs=self.epochs, batch_size=20
        )

        # Geopotential height model.
        # Optimiser.
        opt_height = Adam(lr=1e-6, decay=1e-10, clipvalue=0.6)
        # Create.
        height_model = Sequential()
        height_model.add(
            LSTM(
                400, activation='tanh', input_shape=(past_history, features)
            )
        )
        height_model.add(RepeatVector(future_target))
        height_model.add(LSTM(400, activation='tanh', return_sequences=True))
        height_model.add(LSTM(400, activation='relu', return_sequences=True))
        height_model.add(LSTM(400, activation='relu', return_sequences=True))
        height_model.add(TimeDistributed(Dense(features)))
        height_model.compile(optimizer=opt_height, loss='mse', metrics=['mean_absolute_error'])
        # Train.
        height_model.fit(
            x_height, y_height, epochs=self.epochs, batch_size=20
        )

        # Prediction.
        # Set up inputs.
        predict_temp_input = temperature[-past_history:, :]
        predict_temp_input = predict_temp_input.reshape(
            1, predict_temp_input.shape[0], predict_temp_input.shape[1]
        )
        predict_rh_input = relative_humidity[-past_history:, :]
        predict_rh_input = predict_rh_input.reshape(
            1, predict_rh_input.shape[0], predict_rh_input.shape[1]
        )
        predict_height_input = geopotential_height[-past_history:, :]
        predict_height_input = predict_height_input.reshape(
            1, predict_height_input.shape[0], predict_height_input.shape[1]
        )
        # Make predictions.
        predict_temp = temp_model.predict(predict_temp_input)
        predict_temp = predict_temp[0]
        predict_rh = rh_model.predict(predict_rh_input)
        predict_rh = predict_rh[0]
        predict_height = rh_model.predict(predict_height_input)
        predict_height = predict_height[0]

        # Invert normalisation.
        predict_temp = self.temp_sc.inverse_transform(predict_temp)
        predict_rh = self.rh_sc.inverse_transform(predict_rh)
        predict_height = self.rh_sc.inverse_transform(predict_height)

        # Reshape into 3d arrays.
        # Dimensions.
        len_time = len(predict_temp)
        len_pressure = len(self.pressure_surfaces())
        len_lat = len(self.latitude_lines())
        len_lon = len(self.longitude_lines())
        # Reshape.
        predict_temp = predict_temp.reshape(len_time, len_pressure, len_lat, len_lon)
        predict_rh = predict_rh.reshape(len_time, len_pressure, len_lat, len_lon)
        predict_height = predict_height.reshape(len_time, len_pressure, len_lat, len_lon)

        return predict_temp, predict_rh, predict_height

cdef class Dynamics(RNN):
    """
    AMSIMP Dynamics Class - Also, known as Motus Aeris @ AMSIMP. This class
    generates a simulation of tropospheric and stratsopheric dynamics. 
    Predictions are made by numerically solving the isobaric version of the
    Primitive Equations (they are coupled set of nonlinear PDEs). The initial conditions
    are defined in the class methods of Water, Moist, and Backend. For more information
    on the initial conditions, please see those classes.
    
    Below is a list of the methods included within this class, with a short
    description of their intended purpose. Please see the relevant class methods
    for more information.

    forecast_period ~ generates the period of time over which the forecast will
    be generated for. Please also note that this method also outputs the change in time
    used. This is used in the prognostic equations. 
    simulate ~ generates a simulation of tropospheric and stratsopheric dynamics.
    visualise ~ please explain here.
    """

    def __cinit__(
            self,
            int delta_latitude=5,
            int delta_longitude=5,
            forecast_length=72, 
            delta_t=2, 
            bool ai=True, 
            data_size=150, 
            epochs=200, 
            input_date=None, 
            bool input_data=False,
            psurfaces=None,
            lat=None,
            lon=None,
            height=None, 
            temp=None, 
            rh=None, 
            u=None, 
            v=None,
            dict constants={
                "sidereal_day": (23 + (56 / 60)) * 3600,
                "angular_rotation_rate": ((2 * np.pi) / ((23 + (56 / 60)) * 3600)),
                "planet_radius": constant.R_earth.value,
                "planet_mass": constant.M_earth.value,
                "specific_heat_capacity_psurface": 718,
                "gravitational_acceleration": 9.80665,
                "planet": "Earth"
            }
        ):
        """
        The parameter, forecast_length, defines the length of the 
        simulation (defined in hours). Defaults to a value of 72.

        The parameter, delta_t, defines the change with respect to time
        (in minutes) defined in the simulation. Defaults to a value of 2.
        
        For more information, please refer to amsimp.Backend.__cinit__()
        method.
        """
        # Add units to forecast length variable.
        if type(forecast_length) != Quantity:
            forecast_length = forecast_length * self.units.h
        # Add units to delta_t variable.
        if type(delta_t) != Quantity:
            delta_t = delta_t * self.units.min

        # Declare class variables.
        super().__init__(delta_latitude)
        super().__init__(delta_longitude)
        self.forecast_length = forecast_length
        self.delta_t = delta_t
        self.ai = ai
        super().__init__(input_date)
        super().__init__(input_data)
        super().__init__(height)
        super().__init__(temp)
        super().__init__(rh)
        super().__init__(u)
        super().__init__(v)
        super().__init__(constants)
        
        # Ensure self.forecast_length is greater than, or equal to 1.
        if self.forecast_length.value <= 0:
            raise Exception(
                "forecast_length must be a positive number greater than, or equal to 1. The value of forecast_length was: {}".format(
                    self.forecast_length
                )
            )
        
        # Ensure ai is a boolean value.
        if not isinstance(self.ai, bool):
            raise Exception(
                "ai must be a boolean value. The value of ai was: {}".format(
                    self.ai
                )
            )

    cpdef tuple forecast_period(self):
        """
        Generates the period of time over which the forecast will be generated
        for. 
        
        Explain more.
        """
        segments = int((self.forecast_length / self.delta_t).si.value) + 1

        # Define the forecast period.
        forecast_period = np.linspace(
            0, self.forecast_length.value, segments
        ) * self.forecast_length.unit

        # Convert to seconds.
        delta_t = self.delta_t.to(self.units.s)

        return forecast_period, delta_t

    def __perturbations_errorcheck(self, input_perturbation):
        """
        Explain here.
        """
        if input_perturbation != None:                
            # Check if first index of tuple is a function.
            if not callable(input_perturbation):
                raise Exception(
                    "perturbations must be callable functions."
                )

    def __interpolation_cube(self, input_cube, grid_points):
        """
        Explain here.
        """
        output = input_cube.interpolate(grid_points, iris.analysis.Linear())
        
        return output

    cpdef simulate(
            self, 
            bool save_file=False,
            perturbations_temperature=None,
            perturbations_zonalwind=None,
            perturbations_meridionalwind=None,
            perturbations_mixingratio=None
        ):
        """
        Generates a simulation of tropospheric and stratsopheric dynamics.
        Predictions are made by numerically solving the Primitive
        Equations. Depending on the parameter specifed in the initialisation
        of the class, a recurrent neural network may be incorpated in
        the output.
        
        The Primitive Equations are a set of nonlinear partial differential
        equations that are used to approximate global atmospheric flow and
        are used in most atmospheric models. They consist of three main sets
        of balance equations: a continuity equation, conservation of
        momentum, and a thermal energy equation.

        The Lax–Friedrichs method is used to numerically solve the Primitive
        Equations within the software. It is a numerical method for the
        solution of hyperbolic partial differential equations based on
        finite differences. The method can be described as the 
        FTCS (forward in time, centered in space) scheme with a numerical
        dissipation term of 1/2. One can view the Lax–Friedrichs method as an
        alternative to Godunov's scheme, where one avoids solving a Riemann
        problem at each cell interface

        The parameter, save_file, may be utilised to save the output of
        this class. The output will be saved as a NetCDF file. These
        files can be opened by using the Iris library, which can
        be downloaded via Anaconda Cloud.
        """
        # Error checking.
        np.seterr(all='raise')

        # Ensure save_file is a boolean value.
        if not isinstance(save_file, bool):
            raise Exception(
                "save_file must be a boolean value. The value of save_file was: {}".format(
                    save_file
                )
            )
        
        # Perturbations error checking.
        # Air temperature.
        self.__perturbations_errorcheck(perturbations_temperature)
        # Zonal wind.
        self.__perturbations_errorcheck(perturbations_zonalwind)
        # Meridional wind.
        self.__perturbations_errorcheck(perturbations_meridionalwind)
        # Mixing Ratio.
        self.__perturbations_errorcheck(perturbations_mixingratio)

        # Define variables that do not vary with respect to time.
        cdef np.ndarray latitude = self.latitude_lines().value
        cdef np.ndarray longitude = self.longitude_lines().value
        cdef np.ndarray pressure = self.pressure_surfaces()
        cdef np.ndarray pressure_3d = self.pressure_surfaces(dim_3d=True)
        cdef np.ndarray interpolate_pressure = np.insert(
            pressure, 0, pressure[0] + np.abs(pressure[0] - pressure[1])
        )
        cdef time = self.forecast_period()[0]
        
        # The Coriolis parameter at various latitudes of the Earth's surface,
        # under various approximations.
        # No approximation.
        cdef np.ndarray f = self.coriolis_parameter().value / self.units.s
        f = self.make_3dimensional_array(parameter=f, axis=1)

        # Define the change with respect to time.
        cdef delta_t = self.forecast_period()[1]
        # Forecast length.
        cdef forecast_length = self.forecast_length.to(self.units.s)
        cdef int t = 0

        # Define initial conditions.
        # Gravitational Acceleration.
        cdef g = self.g
        # Geopotential Height.
        cdef np.ndarray height = self.geopotential_height()
        # Wind.
        cdef tuple wind = self.wind()
        # Zonal Wind.
        cdef np.ndarray u = wind[0]
        # Meridional Wind.
        cdef np.ndarray v = wind[1]
        # Vertical Motion.
        cdef np.ndarray omega = self.vertical_motion()
        # Temperature.
        # Air Temperature.
        cdef np.ndarray T = self.temperature()
        # Virtual Temperature.
        cdef np.ndarray T_v = self.virtual_temperature()
        # Relative Humidity.
        cdef np.ndarray rh = self.relative_humidity()
        # Mixing Ratio.
        cdef np.ndarray q = self.mixing_ratio()
        # Precipitable Water.
        cdef np.ndarray pwv = self.precipitable_water()

        # Prediction from Recurrent Neural Network.
        cdef np.ndarray prediction_ai_temp, prediction_ai_height, prediction_ai_rh
        cdef int iterator_ai
        if self.ai:
            prediction_ai = self.model_prediction()
            prediction_ai_height = prediction_ai[2] * self.units.m
            prediction_ai_temp = prediction_ai[0] * self.units.K
            prediction_ai_rh = prediction_ai[1] * self.units.percent
            iterator_ai = 0

        # Define extra variable types.
        # Numy Arrays.
        cdef np.ndarray T_n, q_n, u_n, v_n, height_n, A, B, C, D, E
        cdef np.ndarray RHS, e, T_c, sat_vapor_pressure, geopotential_height
        cdef np.ndarray mean_u, mean_v, mean_T, mean_q, temperature
        cdef np.ndarray virtual_temperature, zonal_wind, meridional_wind
        cdef np.ndarray static_stability, relative_humidity, mixing_ratio
        cdef np.ndarray precipitable_water, Tv_layers, z_0, z_1, z_2, h
        cdef np.ndarray Tv_layer, log_pressure_layer
        # Booleans.
        cdef bool break_loop
        # Ints / Floats.
        cdef int len_p = len(pressure)
        cdef int i
        # Tuples.
        cdef tuple shape, shape_2d

        # Create a bar to determine progress.
        max_bar = len(time)
        bar = IncrementalBar('Progress', max=max_bar)
        # Start progress bar.
        bar.next()

        # Define NumPy array for ouputs.
        shape = (len(time), len(pressure), len(latitude), len(longitude))
        shape_2d = (len(time), len(latitude), len(longitude))
        # Geopotential Height.
        geopotential_height = np.zeros(shape) * height.unit
        geopotential_height[0, :, :, :] = height
        # Temperature.
        # Air Temperature.
        temperature = np.zeros(shape) * T.unit
        temperature[0, :, :, :] = T
        T_n = T.copy()
        # Virtual Temperature.
        virtual_temperature = np.zeros(shape) * T_v.unit
        virtual_temperature[0, :, :, :] = T_v 
        # Geostrophic Wind.
        # Zonal Wind.
        zonal_wind = np.zeros(shape) * u.unit
        zonal_wind[0, :, :, :] = u
        u_n = u.copy()
        # Meridional Wind.
        meridional_wind = np.zeros(shape) * v.unit
        meridional_wind[0, :, :, :] = v
        v_n = v.copy()
        # Vertical Motion.
        vertical_motion = np.zeros(shape) * omega.unit
        vertical_motion[0, :, :, :] = omega
        # Relative Humidity.
        relative_humidity = np.zeros(shape) * rh.unit
        relative_humidity[0, :, :, :] = rh
        # Mixing Ratio.
        mixing_ratio = np.zeros(shape) * q.unit
        mixing_ratio[0, :, :, :] = q
        q_n = q.copy()
        # Precipitable Water.
        precipitable_water = np.zeros(shape_2d) * pwv.unit
        precipitable_water[0, :, :] = pwv

        # Define the initial state for the perturbation functions.
        config = Wind(
            delta_latitude=self.delta_latitude,
            delta_longitude=self.delta_longitude,
            input_data=True,
            psurfaces=self.pressure_surfaces(),
            lat=self.latitude_lines(),
            lon=self.longitude_lines(),
            height=height, 
            temp=T, 
            rh=rh,
            u=u,
            v=v,
            constants=self.constants 
        )

        # Define the coordinates for the cubes. 
        # Latitude.
        lat = DimCoord(
            latitude,
            standard_name='latitude',
            units='degrees'
        )
        # Longitude
        lon = DimCoord(
            longitude,
            standard_name='longitude', 
            units='degrees'
        )
        # Pressure Surfaces.
        p = DimCoord(
            pressure,
            long_name='pressure', 
            units='hPa'
        )
        # Time.
        forecast_period = DimCoord(
            time,
            standard_name='forecast_period', 
            units='hours'
        )
        # Forecast reference time.
        ref_time = AuxCoord(
            self.date.strftime("%Y-%m-%d %H:%M:%S"),
            standard_name='forecast_reference_time'
        )

        # Determine geopotential height of lowest constant pressure surface for
        # the integration of the Hydrostatic Equation. 
        # Define grid to interpolate onto.
        grid_points = [
            ('pressure',  interpolate_pressure.value),
            ('latitude',  latitude),
            ('longitude', longitude),
        ]

        # Define cube.
        z = Cube(height.value,
            standard_name='geopotential_height',
            units='m',
            dim_coords_and_dims=[
                (p, 0), (lat, 1), (lon, 2)
            ],
        )

        # Interpolation of geopotential height based on new grid.
        z = z.interpolate(grid_points, iris.analysis.Linear())
        z_0 = np.asarray(z.data.tolist())[0] * height.unit

        cdef int nt = 1
        try:
            while t < forecast_length.value:
                # Wind
                # Zonal Wind.
                # Determine each term in the zonal momentum equation.
                A = self.gradient_longitude(parameter=-u*u)
                B = self.gradient_latitude(parameter=-v*u)
                C = self.gradient_pressure(parameter=-omega*u)
                D = self.gradient_longitude(parameter=-g*height)
                E = f * v

                # Add any perturbations of zonal wind defined by the user.
                if perturbations_zonalwind != None:
                    perturbations = perturbations_zonalwind(config)
                else:
                    perturbations = np.zeros(u.value.shape) * E.unit

                # Sum the RHS terms and multiple by the time step.
                RHS = (A + B + C + D + E + perturbations) * delta_t
                mean_u = (
                    (
                        u[2:, 1:-1, 1:-1] + u[:-2, 1:-1, 1:-1]
                    ) + (
                        u[1:-1, 2:, 1:-1] + u[1:-1, :-2, 1:-1]
                    ) + (
                        u[1:-1, 1:-1, 2:] + u[1:-1, 1:-1, :-2]
                    )
                ) / 6
                u_n[1:-1, 1:-1, 1:-1] = mean_u + RHS[1:-1, 1:-1, 1:-1]

                # Meridional Wind.
                # Determine each term in the meridional momentum equation.
                A = self.gradient_longitude(parameter=-u*v)
                B = self.gradient_latitude(parameter=-v*v)
                C = self.gradient_pressure(parameter=-omega*v)
                D = self.gradient_latitude(parameter=-g*height)
                E = -f * u

                # Add any perturbations of meridional wind defined by the user.
                if perturbations_meridionalwind != None:
                    perturbations = perturbations_meridionalwind(config)
                else:
                    perturbations = np.zeros(v.value.shape) * E.unit

                # Sum the RHS terms and multiple by the time step.
                RHS = (A + B + C + D + E + perturbations) * delta_t
                mean_v = (
                    (
                        v[2:, 1:-1, 1:-1] + v[:-2, 1:-1, 1:-1]
                    ) + (
                        v[1:-1, 2:, 1:-1] + v[1:-1, :-2, 1:-1]
                    ) + (
                        v[1:-1, 1:-1, 2:] + v[1:-1, 1:-1, :-2]
                    )
                ) / 6
                v_n[1:-1, 1:-1, 1:-1] = mean_v + RHS[1:-1, 1:-1, 1:-1]

                # Temperature.
                # Air Temperature.
                # Thermodynamic Equation.
                A = self.gradient_longitude(parameter=-u*T)
                B = self.gradient_latitude(parameter=-v*T)
                C = self.gradient_pressure(parameter=-omega*T)
                D = omega * ((self.R * T) / (pressure_3d * self.c_p))

                # Add any perturbations of air temperature defined by the user.
                if perturbations_temperature != None:
                    perturbations = perturbations_temperature(config)
                else:
                    perturbations = np.zeros(T.value.shape) * D.unit
                
                # Sum the RHS terms and multiple by the time step.
                RHS = (A + B + C + D + perturbations) * delta_t
                mean_T = (
                    (
                        T[2:, 1:-1, 1:-1] + T[:-2, 1:-1, 1:-1]
                    ) + (
                        T[1:-1, 2:, 1:-1] + T[1:-1, :-2, 1:-1]
                    ) + (
                        T[1:-1, 1:-1, 2:] + T[1:-1, 1:-1, :-2]
                    )
                ) / 6
                T_n[1:-1, 1:-1, 1:-1] = mean_T + RHS[1:-1, 1:-1, 1:-1]

                # Mixing ratio
                # The scalar gradients of the mixing ratio.
                dq_dx = self.gradient_longitude(parameter=q)
                dq_dy = self.gradient_latitude(parameter=q)
                dq_dp = self.gradient_pressure(parameter=q)

                # Advect mixing ratio via wind.
                A = self.gradient_longitude(parameter=-u*q)
                B = self.gradient_latitude(parameter=-v*q)
                C = self.gradient_pressure(parameter=-omega*q)

                # Add any perturbations of mixing ratio defined by the user.
                if perturbations_mixingratio != None:
                    perturbations = perturbations_mixingratio(config)
                else:
                    perturbations = np.zeros(q.value.shape) * C.unit
                
                # Sum the RHS terms and multiple by the time step.
                RHS = (A + B + C + perturbations) * delta_t
                mean_q = (
                    (
                        q[2:, 1:-1, 1:-1] + q[:-2, 1:-1, 1:-1]
                    ) + (
                        q[1:-1, 2:, 1:-1] + q[1:-1, :-2, 1:-1]
                    ) + (
                        q[1:-1, 1:-1, 2:] + q[1:-1, 1:-1, :-2]
                    )
                ) / 6
                q_n[1:-1, 1:-1, 1:-1] = mean_q + RHS[1:-1, 1:-1, 1:-1]

                # Vapor pressure.
                e = pressure_3d * q_n / (0.622 + q_n)

                # Convert temperature in Kelvin to degrees centigrade.
                T_c = T_n.value - 273.15
                # Saturated vapor pressure.
                sat_vapor_pressure = 0.61121 * np.exp(
                    (
                        18.678 - (T_c / 234.5)
                    ) * (T_c / (257.14 + T_c)
                    )
                ) 
                # Add units of measurement.
                sat_vapor_pressure *= self.units.kPa
                sat_vapor_pressure = sat_vapor_pressure.to(self.units.hPa)

                # Relative Humidity.
                rh = (e.value / sat_vapor_pressure.value) * 100
                rh[rh > 100] = 100
                rh[rh < 0] = 0
                rh *= self.units.percent

                # Virtual Temperature.
                T_v = T_n / (
                    1 - (
                    e / pressure_3d
                    ) * (1 - 0.622)
                )

                # Geopotential Height (Hydrostatic Balance).
                # Define variables.
                z_1 = z_0
                height_n = np.zeros(height.value.shape) * height.unit

                # Interpolation of virtual temperature.
                # Define cube.
                cube = Cube(T_v.value,
                    standard_name='virtual_temperature',
                    units='K',
                    dim_coords_and_dims=[
                        (p, 0), (lat, 1), (lon, 2)
                    ],
                )
                # Interpolate.
                cube = cube.interpolate(grid_points, iris.analysis.Linear())
                # Convert to NumPy array
                Tv_layers = np.asarray(cube.data.tolist()) * T_v.unit

                for i in range(len_p):
                    # Virtual temperature layer.
                    Tv_layer = Tv_layers[i:i+2, :, :]

                    # Log pressure of layer.
                    log_pressure_layer = np.log(interpolate_pressure[i:i+2].value)

                    # Determine thickness of layer.
                    h = (
                        -self.R / self.g * np.trapz(
                            Tv_layer, x=log_pressure_layer, axis=0
                        )
                    )

                    # Determine geopotential height of constant pressure
                    # surface.
                    z_2 = z_1 + h

                    # Redefine bottom pressure surface.
                    z_1 = z_2

                    # Add to NumPy array.
                    height_n[i, :, :] = z_2

                # Recurrent Neural Network.
                if self.ai:
                    # Apply predictions every six hours.
                    if t != 0 and t % 21600 == 0:
                        # Geopotential Height.
                        height = (height + prediction_ai_height[iterator_ai]) / 2
                        # Temperature. 
                        T = (T + prediction_ai_temp[iterator_ai]) / 2
                        # Relative Humidity.
                        rh = (rh + prediction_ai_rh[iterator_ai]) / 2
                        # Iterator.
                        iterator_ai += 1
                
                # Configure the Wind class, so, that it aligns with the
                # paramaters defined by the user.
                break_loop = False
                while not break_loop:
                    try:
                        config = Wind(
                            delta_latitude=self.delta_latitude,
                            delta_longitude=self.delta_longitude,
                            input_data=True,
                            psurfaces=self.pressure_surfaces(),
                            lat=self.latitude_lines(),
                            lon=self.longitude_lines(),
                            height=height_n, 
                            temp=T_n, 
                            rh=rh,
                            u=u_n,
                            v=v_n,
                            constants=self.constants 
                        )
                        break_loop = True
                    except:
                        pass

                # Define other prognostic variables.
                # Precipitable Water.
                pwv = config.precipitable_water()

                # Add predictions to NumPy arrays.
                # Geopotential Height.
                geopotential_height[nt, :, :, :] = height_n
                height = height_n
                # Geostrophic Wind.
                # Zonal Wind.
                zonal_wind[nt, :, :, :] = u_n
                u = u_n
                # Meridional Wind.
                meridional_wind[nt, :, :, :] = v_n
                v = v_n
                # Temperature.
                # Air Temperature.
                temperature[nt, :, :, :] = T_n
                T = T_n
                # Virtual Temperature.
                virtual_temperature[nt, :, :, :] = T_v
                # Relative Humidity.
                relative_humidity[nt, :, :, :] = rh
                # Mixing Ratio.
                mixing_ratio[nt, :, :, :] = q_n
                q = q_n
                # Precipitable Water.
                precipitable_water[nt, :, :] = pwv

                # Add time step.
                t += delta_t.value
                nt += 1
                bar.next()
        except KeyboardInterrupt:
            pass

        # Cubes.
        grid_points = [
            ('forecast_period', time.value), 
            ('pressure',  pressure.value),
            ('latitude',  latitude),
            ('longitude', longitude),
        ]
        grid_points_pwv = [
            ('forecast_period', time.value),
            ('latitude',  latitude),
            ('longitude', longitude),
        ]

        # Geopotential Height Cube.
        height_cube = Cube(geopotential_height[:, 1:-1, 1:-1, 1:-1],
            standard_name='geopotential_height',
            units='m',
            dim_coords_and_dims=[
                (forecast_period, 0), (p[1:-1], 1), (lat[1:-1], 2), (lon[1:-1], 3)
            ],
            attributes={
                'source': 'Motus Aeris @ AMSIMP',
            }
        )
        height_cube = self.__interpolation_cube(
            input_cube=height_cube, grid_points=grid_points
        )
        height_cube.add_aux_coord(ref_time)
        # Wind Cubes.
        # Zonal Wind Cube.
        u_cube = Cube(zonal_wind[:, 1:-1, 1:-1, 1:-1],
            standard_name='x_wind',
            units='m s-1',
            dim_coords_and_dims=[
                (forecast_period, 0), (p[1:-1], 1), (lat[1:-1], 2), (lon[1:-1], 3)
            ],
            attributes={
                'source': 'Motus Aeris @ AMSIMP',
            }
        )
        u_cube = self.__interpolation_cube(
            input_cube=u_cube, grid_points=grid_points
        )
        u_cube.add_aux_coord(ref_time)
        # Meridional Wind Cube.
        v_cube = Cube(meridional_wind[:, 1:-1, 1:-1, 1:-1],
            standard_name='y_wind',
            units='m s-1',
            dim_coords_and_dims=[
                (forecast_period, 0), (p[1:-1], 1), (lat[1:-1], 2), (lon[1:-1], 3)
            ],
            attributes={
                'source': 'Motus Aeris @ AMSIMP',
            }
        )
        v_cube = self.__interpolation_cube(
            input_cube=v_cube, grid_points=grid_points
        )
        v_cube.add_aux_coord(ref_time)
        # Temperature.
        # Air Temperature.
        T_cube = Cube(temperature[:, 1:-1, 1:-1, 1:-1],
            standard_name='air_temperature',
            units='K',
            dim_coords_and_dims=[
                (forecast_period, 0), (p[1:-1], 1), (lat[1:-1], 2), (lon[1:-1], 3)
            ],
            attributes={
                'source': 'Motus Aeris @ AMSIMP',
            }
        )
        T_cube = self.__interpolation_cube(
            input_cube=T_cube, grid_points=grid_points
        )
        T_cube.add_aux_coord(ref_time)
        # Virtual Temperature.
        Tv_cube = Cube(virtual_temperature[:, 1:-1, 1:-1, 1:-1],
            standard_name='virtual_temperature',
            units='K',
            dim_coords_and_dims=[
                (forecast_period, 0), (p[1:-1], 1), (lat[1:-1], 2), (lon[1:-1], 3)
            ],
            attributes={
                'source': 'Motus Aeris @ AMSIMP',
            }
        )
        Tv_cube = self.__interpolation_cube(
            input_cube=Tv_cube, grid_points=grid_points
        )
        Tv_cube.add_aux_coord(ref_time)
        # Relative Humidity.
        rh_cube = Cube(relative_humidity[:, 1:-1, 1:-1, 1:-1],
            standard_name='relative_humidity',
            units='%',
            dim_coords_and_dims=[
                (forecast_period, 0), (p[1:-1], 1), (lat[1:-1], 2), (lon[1:-1], 3)
            ],
            attributes={
                'source': 'Motus Aeris @ AMSIMP',
            }
        )
        rh_cube = self.__interpolation_cube(
            input_cube=rh_cube, grid_points=grid_points
        )
        rh_cube.add_aux_coord(ref_time)
        # Mixing Ratio.
        q_cube = Cube(mixing_ratio[:, 1:-1, 1:-1, 1:-1],
            long_name='mixing_ratio',
            dim_coords_and_dims=[
                (forecast_period, 0), (p[1:-1], 1), (lat[1:-1], 2), (lon[1:-1], 3)
            ],
            attributes={
                'source': 'Motus Aeris @ AMSIMP',
            }
        )
        q_cube = self.__interpolation_cube(
            input_cube=q_cube, grid_points=grid_points
        )
        q_cube.add_aux_coord(ref_time)
        # Precipitable Water.
        pwv_cube = Cube(precipitable_water[:, 1:-1, 1:-1],
            long_name='precipitable_water',
            units='mm',
            dim_coords_and_dims=[
                (forecast_period, 0), (lat[1:-1], 1), (lon[1:-1], 2)
            ],
            attributes={
                'source': 'Motus Aeris @ AMSIMP',
            }
        )
        pwv_cube = self.__interpolation_cube(
            input_cube=pwv_cube, grid_points=grid_points_pwv
        )
        pwv_cube.add_aux_coord(ref_time)

        # Create Cube list of output parameters.
        output = CubeList([
            height_cube,
            u_cube,
            v_cube,
            T_cube,
            Tv_cube,
            rh_cube,
            q_cube,
            pwv_cube
        ])

        # If specified, save the forecast in the file format, .nc.
        if save_file:
            # Establish the file name.
            filename = 'motusaeris_amsimp_' + str(self.date.year)
            filename += str(self.date.month) + str(self.date.day)
            filename += str(self.date.hour) + '.nc'

            # Save.
            save(output, filename)

        # Finish progress bar.
        bar.finish()

        return output

    def visualise(
            self,
            data=None,
            plot=[
                "air_temperature",
                "geopotential_height",
                "precipitable_water",
                "relative_humidity"
            ],
            psurface=1000
        ):
        """
        Explain here.

        Need to fix.
        """
        # Declare variable types.
        # NumPy arrays
        cdef np.ndarray time, lat, lon, longitude, data1, data2, data3, data4
        cdef np.ndarray level1, level2, level3, level4
        # Floats.
        cdef float min1, min2, min3, min4, max1, max2, max3, max4

        if self.planet == "Earth":
            # Error checking.
            # Ensure a dataset is provided.
            if data == None:
                raise Exception("The dataset for simulation must be defined.")

            # Pressure Surface.
            try:
                pressure = data[0].coords("pressure")[0].points
            except:
                pressure = data[1].coords("pressure")[0].points
            if psurface < pressure[-1] or psurface > pressure[0]:
                raise Exception(
                    "psurface must be a real number within the isobaric boundaries. The value of psurface was: {}".format(
                        psurface
                    )
                )

            # Index of the nearest pressure surface in amsimp.Backend.pressure_surfaces()
            indx_psurface = (np.abs(pressure - psurface)).argmin()

            # Define the forecast period.
            time = data[0].coords("forecast_period")[0].points
            time_unit = str(data[0].coords("forecast_period")[0].units)

            # Grid.
            lat = data[0].coords("latitude")[0].points
            lon = data[0].coords("longitude")[0].points
            longitude = lon

            # Style of graph.
            style.use("fivethirtyeight")

            # Define layout.
            gs = gridspec.GridSpec(2, 2)
            fig = plt.figure(figsize=(18.5, 7.5))
            fig.subplots_adjust(hspace=0.340, bottom=0.105, top=0.905)
            plt.ion()

            # Graph 1.
            ax1 = plt.subplot(gs[0, 0], projection=ccrs.EckertIII())
            label1 = plot[0]
            data1 = data.extract(label1)[0].data
            if data1.ndim == 3: 
                data1 = data1;
            else: 
                data1 = data1[:, indx_psurface, :, :];
            data1 = np.asarray(data1.tolist())
            unit1 = " (" + str(data.extract(label1)[0].metadata.units) + ")"
            label1 = label1.replace("_", " ").title()
            # Graph 2.
            ax2 = plt.subplot(gs[1, 0], projection=ccrs.EckertIII())
            label2 = plot[1]
            data2 = data.extract(label2)[0].data
            if data2.ndim == 3:
                data2 = data2;
            else: 
                data2 = data2[:, indx_psurface, :, :];
            data2 = np.asarray(data2.tolist())
            unit2 = " (" + str(data.extract(label2)[0].metadata.units) + ")"
            label2 = label2.replace("_", " ").title()
            # Graph 3.
            ax3 = plt.subplot(gs[0, 1], projection=ccrs.EckertIII())
            label3 = plot[2]
            data3 = data.extract(label3)[0].data
            if data3.ndim == 3:
                data3 = data3;
            else:
                data3 = data3[:, indx_psurface, :, :];
            data3 = np.asarray(data3.tolist())
            unit3 = " (" + str(data.extract(label3)[0].metadata.units) + ")"
            label3 = label3.replace("_", " ").title()
            # Graph 4.
            ax4 = plt.subplot(gs[1, 1], projection=ccrs.EckertIII())
            label4 = plot[3]
            data4 = data.extract(label4)[0].data[:, indx_psurface, :, :]
            if data4.ndim == 3: 
                data4 = data4; 
            else: 
                data4 = data4[:, indx_psurface, :, :];
            data4 = np.asarray(data4.tolist())
            unit4 = " (" + str(data.extract(label4)[0].metadata.units) + ")"
            label4 = label4.replace("_", " ").title()

            # Get date.
            date = data[0].coords("forecast_reference_time")[0].points[0]

            for i in range(len(time)):
                # Plot 1.
                # EckertIII projection details.
                ax1.set_global()
                ax1.coastlines()
                ax1.gridlines()

                # Add SALT to the graph.
                ax1.set_xlabel("Longitude ($\lambda$)")
                ax1.set_ylabel("Latitude ($\phi$)")
                ax1.set_title(label1)

                # Contour plot.
                cmap1 = plt.get_cmap("hot")
                min1 = np.min(data1)
                max1 = np.max(data1)
                level1 = np.linspace(min1, max1, 21)
                data1_plot, lon = add_cyclic_point(data1, coord=longitude)
                contour1 = ax1.contourf(
                    lon,
                    lat,
                    data1_plot[i, :, :],
                    cmap=cmap1,
                    levels=level1,
                    transform=ccrs.PlateCarree()
                )

                # Checks for a colorbar.
                if i == 0:
                    cb1 = fig.colorbar(contour1, ax=ax1)
                    tick_locator = ticker.MaxNLocator(nbins=10)
                    cb1.locator = tick_locator
                    cb1.update_ticks()
                    cb1.set_label(unit1)

                # Plot 2.
                cmap2 = plt.get_cmap("jet")
                min2 = np.min(data2)
                max2 = np.max(data2)
                level2 = np.linspace(min2, max2, 21)
                data2_plot, lon = add_cyclic_point(data2, coord=longitude)
                contour2 = ax2.contourf(
                    lon,
                    lat, 
                    data2_plot[i, :, :], 
                    cmap=cmap2, 
                    levels=level2,
                    transform=ccrs.PlateCarree()
                )

                # EckertIII projection details.
                ax2.set_global()
                ax2.coastlines()
                ax2.gridlines()

                # Checks for a colorbar.
                if i == 0:
                    cb2 = fig.colorbar(contour2, ax=ax2)
                    tick_locator = ticker.MaxNLocator(nbins=10)
                    cb2.locator = tick_locator
                    cb2.update_ticks()
                    cb2.set_label(unit2)

                # Add SALT to the graph.
                ax2.set_xlabel("Longitude ($\lambda$)")
                ax2.set_ylabel("Latitude ($\phi$)")
                ax2.set_title(label2)

                # Plot 3.
                cmap3 = plt.get_cmap("seismic")
                min3 = np.min(data3)
                max3 = np.max(data3)
                level3 = np.linspace(min3, max3, 21)
                data3_plot, lon = add_cyclic_point(data3, coord=longitude)
                contour3 = ax3.contourf(
                    lon, 
                    lat, 
                    data3_plot[i, :, :], 
                    cmap=cmap3, 
                    levels=level3,
                    transform=ccrs.PlateCarree()
                )

                # EckertIII projection details.
                ax3.set_global()
                ax3.coastlines()
                ax3.gridlines()

                # Checks for a colorbar.
                if i == 0:
                    cb3 = fig.colorbar(contour3, ax=ax3)
                    tick_locator = ticker.MaxNLocator(nbins=10)
                    cb3.locator = tick_locator
                    cb3.update_ticks()
                    cb3.set_label(unit3)

                # Add SALT to the graph.
                ax3.set_xlabel("Longitude ($\lambda$)")
                ax3.set_ylabel("Latitude ($\phi$)")
                ax3.set_title(label3)

                # Plot 4.
                min4 = np.min(data4)
                max4 = np.max(data4)
                level4 = np.linspace(min4, max4, 21)
                data4_plot, lon = add_cyclic_point(data4, coord=longitude)
                contour4 = ax4.contourf(
                    lon, 
                    lat, 
                    data4_plot[i, :, :], 
                    levels=level4,
                    transform=ccrs.PlateCarree()
                )

                # EckertIII projection details.
                ax4.set_global()
                ax4.coastlines()
                ax4.gridlines()

                # Checks for a colorbar.
                if i == 0:
                    cb4 = fig.colorbar(contour4, ax=ax4)
                    tick_locator = ticker.MaxNLocator(nbins=10)
                    cb4.locator = tick_locator
                    cb4.update_ticks()
                    cb4.set_label(unit4)

                # Add SALT to the graph.
                ax4.set_xlabel("Longitude ($\lambda$)")
                ax4.set_ylabel("Latitude ($\phi$)")
                ax4.set_title(label4)

                # Title
                title = (
                    "Motus Aeris @ AMSIMP (" + date + ", +" + str(np.round(time[i], 2)) + time_unit + ", " + str(pressure[indx_psurface]) + " hPa)"
                )
                fig.suptitle(title)

                # Check if folder needs to created, and if so create one.
                try:
                    os.mkdir('visualisations')
                except OSError:
                    pass

                # Save image to folder.
                plt.savefig('visualisations/timestep_' + str(i), dpi=300)

                # Displaying simualtion.
                plt.show()
                plt.pause(0.01)
                if i < (len(time) - 1):
                    ax1.clear()
                    ax2.clear()
                    ax3.clear()
                    ax4.clear()
                else:
                    # Compile photos into video using ffmeg.
                    os.system("ffmpeg -framerate 30 -i visualisations/timestep_%d.png -vf format=yuv420p visualise.mp4")
                    
                    # Remove images.
                    for i in range(len(time)):
                        os.remove("visualisations/timestep_" + str(i) + ".png")

                    # Remove folder.
                    try:
                        os.rmdir('visualisations')
                    except OSError:
                        pass
        else:
            raise NotImplementedError(
                "Visualisations for planetary bodies other than Earth is not currently implemented."
            )