Data and code with comprehensive comments for articles on The Seasonal Variation of the DC Global Electric Circuit
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 

8.7 KiB

Data preprocessing for further calculations

Import libraries

In [2]:
import datetime as dt

import numpy as np
import pandas as pd
In [ ]:
# also available at https://eee.ipfran.ru/files/seasonal-variation-2024/
# attention: the files are very large (~ 350 GB totally)
src_path = "../shared_files/eee_public_files/seasonal-variation-2024/"

Preprocessing WRF T2m data

In [42]:
# available numbers of simulated days for analysis
wrf_N_days = 4992
inmcm_N_days = 3650
In [43]:
# dates corresponding to the indices (0 axis) of the data arrays
# note: for WRF dates correspond to real dates

wrf_dt_indicies = np.array(
    [dt.date(1980, 1, 1) + dt.timedelta(i * 3) for i in range(wrf_N_days)]
)
inmcm_dt_indicies = np.array(
    [dt.date(2022, 1, 1) + dt.timedelta(i % 365) for i in range(inmcm_N_days)]
)
In [147]:
wrf_T2_data = np.load(f"{src_path}/T2-MAP-FULL.npy")[:wrf_N_days]
wrf_T2_data.shape
Out[147]:
(4992, 180, 360)
In [149]:
# mean surface air temperature values for different latitudes and months
wrf_mon_T2 = np.zeros((180, 12))

for month_idx in range(12):
    monthly_indicies = [
        i for i, date in enumerate(wrf_dt_indicies) if date.month == month_idx + 1
    ]  # indicies of days available for `month_idx+1` month

    wrf_mon_T2[:, month_idx] = wrf_mean_T2[monthly_indicies].mean(axis=0)
In [164]:
np.save(f"./data/WRF/WRF_T2_LATxMON.npy",wrf_mon_T2)

INMCM and WRF IP: classic parametrization

In [72]:
wrf_daily_latitudal_ip = {}
inmcm_daily_latitudal_ip = {}

wrf_hourly_total_ip = {}
inmcm_hourly_total_ip = {}
In [73]:
for cape_thres in [800, 1000, 1200]: # J/kg
    print(cape_thres)

    # grid cell contributions to the IP (not normalised) with the shape
    # (number of days, number of hours, number of latitudes, number of longitudes)
    wrf_raw_ip_data = np.load(f"{src_path}/WRF-IP-MAP-{cape_thres}.npy")[:wrf_N_days]
    wrf_raw_ip_data = wrf_raw_ip_data[:, :24, :, :]
    wrf_raw_ip_data /= (1/240e3) * wrf_raw_ip_data.sum(axis=(-2,-1)).mean()

    wrf_daily_latitudal_ip[cape_thres] = wrf_raw_ip_data.mean(axis=1).sum(axis=-1)
    wrf_hourly_total_ip[cape_thres] = wrf_raw_ip_data.sum(axis=(-2, -1))

    inmcm_raw_ip_data = np.load(f"{src_path}/INMCM-IP-MAP-{cape_thres}.npy").reshape((inmcm_N_days, 24, 120, 180))
    inmcm_raw_ip_data /= (1/240e3) * inmcm_raw_ip_data.sum(axis=(-2,-1)).mean()
    
    inmcm_daily_latitudal_ip[cape_thres] = inmcm_raw_ip_data.mean(axis=1).sum(axis=-1)
    inmcm_hourly_total_ip[cape_thres] = inmcm_raw_ip_data.sum(axis=(-2, -1))

    del wrf_raw_ip_data
    del inmcm_raw_ip_data
800
1000
1200
In [89]:
for cape_thres in [800, 1000, 1200]:  # J/kg
    np.save(
        f"./data/INMCM/INMCM_HOURLY_TOTAL_IP_{cape_thres}.npy",
        inmcm_hourly_total_ip[cape_thres],
    )
    np.save(
        f"./data/WRF/WRF_HOURLY_TOTAL_IP_{cape_thres}.npy",
        wrf_hourly_total_ip[cape_thres],
    )

    wrf_data_LATxMON = np.zeros((180, 12))
    inmcm_data_LATxMON = np.zeros((120, 12))

    for month_idx in range(12):
        monthly_indicies = [
            i for i, date in enumerate(wrf_dt_indicies) if date.month == month_idx + 1
        ]  # indicies of days available for `month_idx+1` month

        wrf_data_MONxLAT[:, month_idx] = wrf_daily_latitudal_ip[cape_thres][monthly_indicies].mean(
            axis=0
        )
    np.save(
        f"./data/WRF/WRF_IP_{cape_thres}_LATxMON.npy",
        wrf_data_MONxLAT,
    )

    for month_idx in range(12):
        monthly_indicies = [
            i for i, date in enumerate(inmcm_dt_indicies) if date.month == month_idx + 1
        ]  # indicies of days available for `month_idx+1` month

        inmcm_data_LATxMON[:, month_idx] = inmcm_daily_latitudal_ip[cape_thres][
            monthly_indicies
        ].mean(axis=0)
    np.save(
        f"./data/INMCM/INMCM_IP_{cape_thres}_LATxMON.npy",
        inmcm_data_LATxMON,
    )

WRF IP: parametrization based on T2

In [98]:
wrf_raw_ip_data = np.load(f"WRF-IP-MAP-500-T2-25.npy")[:wrf_N_days]
wrf_raw_ip_data = wrf_raw_ip_data[:, :24, :, :]
wrf_raw_ip_data /= (1/240e3) * wrf_raw_ip_data.sum(axis=(-2,-1)).mean()

wrf_daily_latitudal_ip = wrf_raw_ip_data.mean(axis=1).sum(axis=-1)
wrf_hourly_total_ip = wrf_raw_ip_data.sum(axis=(-2, -1))

np.save(
    f"./data/WRF/WRF_HOURLY_TOTAL_IP_500_T2_25.npy",
    wrf_hourly_total_ip,
)

wrf_data_LATxMON = np.zeros((180, 12))

for month_idx in range(12):
    monthly_indicies = [
        i for i, date in enumerate(wrf_dt_indicies) if date.month == month_idx + 1
    ]

    wrf_data_MONxLAT[:, month_idx] = wrf_daily_latitudal_ip[monthly_indicies].mean(
        axis=0
    )
np.save(
    f"./data/WRF/WRF_IP_500_T2_25_LATxMON.npy",
    wrf_data_MONxLAT
)