Scripts for IP simulations

0_prepare_data.ipynb

@@ -0,0 +1,305 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "98e6e23d-5ca7-4706-b1d9-dd57b54888ef",
+   "metadata": {},
+   "source": [
+    "# Data preprocessing for further calculations"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5324ceb9-24e7-454b-87b9-ba9a717078ae",
+   "metadata": {},
+   "source": [
+    "### Import libraries"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "7b2a7f44-b0cb-4471-a0c6-e56da23caf86",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import datetime as dt\n",
+    "\n",
+    "import numpy as np\n",
+    "import pandas as pd"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "36b9f49e-32e6-4544-a9d3-f6a8ba49d867",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# also available at https://eee.ipfran.ru/files/seasonal-variation-2024/\n",
+    "# attention: the files are very large (~ 350 GB totally)\n",
+    "src_path = \"../shared_files/eee_public_files/seasonal-variation-2024/\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5e16ee8e-f3b0-4251-9691-19d7dfd4aff7",
+   "metadata": {},
+   "source": [
+    "## Preprocessing WRF T2m data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 42,
+   "id": "78a4350c-59fb-479a-b7cd-e2bf9b996d36",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# available numbers of simulated days for analysis\n",
+    "wrf_N_days = 4992\n",
+    "inmcm_N_days = 3650"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 43,
+   "id": "53cb9cc3-0e56-4da4-920b-2f071a0846fb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# dates corresponding to the indices (0 axis) of the data arrays\n",
+    "# note: for WRF dates correspond to real dates\n",
+    "\n",
+    "wrf_dt_indicies = np.array(\n",
+    "    [dt.date(1980, 1, 1) + dt.timedelta(i * 3) for i in range(wrf_N_days)]\n",
+    ")\n",
+    "inmcm_dt_indicies = np.array(\n",
+    "    [dt.date(2022, 1, 1) + dt.timedelta(i % 365) for i in range(inmcm_N_days)]\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 147,
+   "id": "0e4d0268-208f-45c7-bb0a-b5e8d3c7c1f7",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(4992, 180, 360)"
+      ]
+     },
+     "execution_count": 147,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "wrf_T2_data = np.load(f\"{src_path}/T2-MAP-FULL.npy\")[:wrf_N_days]\n",
+    "wrf_T2_data.shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 149,
+   "id": "ec569ffd-93c2-4490-8ba1-69af4fab8f23",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# mean surface air temperature values for different latitudes and months\n",
+    "wrf_mon_T2 = np.zeros((180, 12))\n",
+    "\n",
+    "for month_idx in range(12):\n",
+    "    monthly_indicies = [\n",
+    "        i for i, date in enumerate(wrf_dt_indicies) if date.month == month_idx + 1\n",
+    "    ]  # indicies of days available for `month_idx+1` month\n",
+    "\n",
+    "    wrf_mon_T2[:, month_idx] = wrf_mean_T2[monthly_indicies].mean(axis=0)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 164,
+   "id": "08302aa4-cb14-47f9-8216-a9db06ae53ef",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.save(f\"./data/WRF/WRF_T2_LATxMON.npy\",wrf_mon_T2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "46d4f093-a420-42c7-b885-a8409d9d8ee4",
+   "metadata": {},
+   "source": [
+    "### INMCM and WRF IP: classic parametrization"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 72,
+   "id": "94a603c3-982d-4c78-be1c-bb6c74b86b5b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "wrf_daily_latitudal_ip = {}\n",
+    "inmcm_daily_latitudal_ip = {}\n",
+    "\n",
+    "wrf_hourly_total_ip = {}\n",
+    "inmcm_hourly_total_ip = {}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 73,
+   "id": "d8e43c4f-59af-483c-8979-535c696abb4e",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "800\n",
+      "1000\n",
+      "1200\n"
+     ]
+    }
+   ],
+   "source": [
+    "for cape_thres in [800, 1000, 1200]: # J/kg\n",
+    "    print(cape_thres)\n",
+    "\n",
+    "    # grid cell contributions to the IP (not normalised) with the shape\n",
+    "    # (number of days, number of hours, number of latitudes, number of longitudes)\n",
+    "    wrf_raw_ip_data = np.load(f\"{src_path}/WRF-IP-MAP-{cape_thres}.npy\")[:wrf_N_days]\n",
+    "    wrf_raw_ip_data = wrf_raw_ip_data[:, :24, :, :]\n",
+    "    wrf_raw_ip_data /= (1/240e3) * wrf_raw_ip_data.sum(axis=(-2,-1)).mean()\n",
+    "\n",
+    "    wrf_daily_latitudal_ip[cape_thres] = wrf_raw_ip_data.mean(axis=1).sum(axis=-1)\n",
+    "    wrf_hourly_total_ip[cape_thres] = wrf_raw_ip_data.sum(axis=(-2, -1))\n",
+    "\n",
+    "    inmcm_raw_ip_data = np.load(f\"{src_path}/INMCM-IP-MAP-{cape_thres}.npy\").reshape((inmcm_N_days, 24, 120, 180))\n",
+    "    inmcm_raw_ip_data /= (1/240e3) * inmcm_raw_ip_data.sum(axis=(-2,-1)).mean()\n",
+    "    \n",
+    "    inmcm_daily_latitudal_ip[cape_thres] = inmcm_raw_ip_data.mean(axis=1).sum(axis=-1)\n",
+    "    inmcm_hourly_total_ip[cape_thres] = inmcm_raw_ip_data.sum(axis=(-2, -1))\n",
+    "\n",
+    "    del wrf_raw_ip_data\n",
+    "    del inmcm_raw_ip_data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 89,
+   "id": "eb28cbc7-eb0a-49be-8cc1-734bba1d06f5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for cape_thres in [800, 1000, 1200]:  # J/kg\n",
+    "    np.save(\n",
+    "        f\"./data/INMCM/INMCM_HOURLY_TOTAL_IP_{cape_thres}.npy\",\n",
+    "        inmcm_hourly_total_ip[cape_thres],\n",
+    "    )\n",
+    "    np.save(\n",
+    "        f\"./data/WRF/WRF_HOURLY_TOTAL_IP_{cape_thres}.npy\",\n",
+    "        wrf_hourly_total_ip[cape_thres],\n",
+    "    )\n",
+    "\n",
+    "    wrf_data_LATxMON = np.zeros((180, 12))\n",
+    "    inmcm_data_LATxMON = np.zeros((120, 12))\n",
+    "\n",
+    "    for month_idx in range(12):\n",
+    "        monthly_indicies = [\n",
+    "            i for i, date in enumerate(wrf_dt_indicies) if date.month == month_idx + 1\n",
+    "        ]  # indicies of days available for `month_idx+1` month\n",
+    "\n",
+    "        wrf_data_MONxLAT[:, month_idx] = wrf_daily_latitudal_ip[cape_thres][monthly_indicies].mean(\n",
+    "            axis=0\n",
+    "        )\n",
+    "    np.save(\n",
+    "        f\"./data/WRF/WRF_IP_{cape_thres}_LATxMON.npy\",\n",
+    "        wrf_data_MONxLAT,\n",
+    "    )\n",
+    "\n",
+    "    for month_idx in range(12):\n",
+    "        monthly_indicies = [\n",
+    "            i for i, date in enumerate(inmcm_dt_indicies) if date.month == month_idx + 1\n",
+    "        ]  # indicies of days available for `month_idx+1` month\n",
+    "\n",
+    "        inmcm_data_LATxMON[:, month_idx] = inmcm_daily_latitudal_ip[cape_thres][\n",
+    "            monthly_indicies\n",
+    "        ].mean(axis=0)\n",
+    "    np.save(\n",
+    "        f\"./data/INMCM/INMCM_IP_{cape_thres}_LATxMON.npy\",\n",
+    "        inmcm_data_LATxMON,\n",
+    "    )"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "91bc6d7a-393c-4078-9a6d-1955393d55f5",
+   "metadata": {},
+   "source": [
+    "### WRF IP: parametrization based on T2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 98,
+   "id": "2b6f987e-ba3c-4371-af7b-c9857a7d33d9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "wrf_raw_ip_data = np.load(f\"WRF-IP-MAP-500-T2-25.npy\")[:wrf_N_days]\n",
+    "wrf_raw_ip_data = wrf_raw_ip_data[:, :24, :, :]\n",
+    "wrf_raw_ip_data /= (1/240e3) * wrf_raw_ip_data.sum(axis=(-2,-1)).mean()\n",
+    "\n",
+    "wrf_daily_latitudal_ip = wrf_raw_ip_data.mean(axis=1).sum(axis=-1)\n",
+    "wrf_hourly_total_ip = wrf_raw_ip_data.sum(axis=(-2, -1))\n",
+    "\n",
+    "np.save(\n",
+    "    f\"./data/WRF/WRF_HOURLY_TOTAL_IP_500_T2_25.npy\",\n",
+    "    wrf_hourly_total_ip,\n",
+    ")\n",
+    "\n",
+    "wrf_data_LATxMON = np.zeros((180, 12))\n",
+    "\n",
+    "for month_idx in range(12):\n",
+    "    monthly_indicies = [\n",
+    "        i for i, date in enumerate(wrf_dt_indicies) if date.month == month_idx + 1\n",
+    "    ]\n",
+    "\n",
+    "    wrf_data_MONxLAT[:, month_idx] = wrf_daily_latitudal_ip[monthly_indicies].mean(\n",
+    "        axis=0\n",
+    "    )\n",
+    "np.save(\n",
+    "    f\"./data/WRF/WRF_IP_500_T2_25_LATxMON.npy\",\n",
+    "    wrf_data_MONxLAT\n",
+    ")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

readme.md

@@ -3,12 +3,14 @@
 * `1_Earlier_measurements_images.ipynb` plots seasonal variations from external sources
 * `2_Vostok_measurements_images.ipynb` plots seasonal variations and seasonal-dirunal diagram using new and early Vostok PG measurements
 * `3_WRF_T2_images.ipynb` plots seasonal variation of `T2m` temperature averaged across different latitude bands
+* `4_IP_simulations_temporal_images.ipynb` plots seasonal variation of simulated IP grouped by datasets and different year ranges
+* `5_IP_simulations_spatial_images.ipynb` plots seasonal variation of simulated IP grouped by latitude ranges
 
 # Detailed Description of the Scripts
 
 ## Script `1_Earlier_measurements_images.ipynb`
 
-This program contains digitized data from external sources, necessary for constructing Figure 1.
+This program contains digitized data from external sources, necessary for constructing Figure 1.1.
 
 At the beginning of the script, the necessary libraries are loaded and arrays with digitized data are declared; at the end, a graph is constructed. 
 
@@ -42,9 +44,9 @@ Both described functions are used to compute values necessary for plotting graph
 
 For both new and early Vostok data, we apply the `pass_fair_weather` function, resulting in two datasets that contain only the hours of fair weather days (`df` and `earlier_df`)
 
-### Figure 2
+### Figure 1.2
 
-To construct Figure 2, using the prepared data and helper functions, we calculate the mean values, the count of fair weather days and standard errors for three sets of data:
+To construct Figure 1.2, using the prepared data and helper functions, we calculate the mean values, the count of fair weather days and standard errors for three sets of data:
 
 1. The complete series of new Vostok data.
 2. The same series up to and including the year 2012.
@@ -52,15 +54,15 @@ To construct Figure 2, using the prepared data and helper functions, we calculat
 
 > **_Note:_** The data from this figure is saved in the temporary file `vostok_2006_2020_results.npz` for use in the second article. This helps avoid code duplication or merging code to build different entities in a single cumbersome file.
 
-### Figure 3
+### Figure 1.3
 
-To construct Figure 3, we transform the Vostok data series into a matrix of 12 months x 24 hours. To do this, we group the original dataframe of fair weather hours by months and hours, and then find the mean value for all data points taken at a specific hour of a specific month (saved in dataframe `sd_df`).
+To construct Figure 1.3, we transform the Vostok data series into a matrix of 12 months x 24 hours. To do this, we group the original dataframe of fair weather hours by months and hours, and then find the mean value for all data points taken at a specific hour of a specific month (saved in dataframe `sd_df`).
 
 For clarity, we also present slices of this diurnal-seasonal diagram at 3, 9, 15, and 21 hours UTC.
 
 > **_Note:_**  Renaming the axes of the multi-index resulting from grouping (`sd_df.index.set_names(['hour', 'month'], inplace=True)`) is not necessary for the code and can be commented out; however, it may be convenient for further work with the diurnal-seasonal dataframe `sd_df`.
 
-### Figure 4
+### Figure 1.4
 #### Removal of field anomalies associated with meteorological parameters
 
 First, we load the meteorological datasets (`temp_df`, `wind_df`, `pressure_df`), averaged by days (`vostok_daily_temp`, `vostok_daily_wind`, `vostok_daily_pressure_mm_hg`). For further analysis, we use the `meteo_df` dataframe, which is created by merging the dataframe with daily average potential gradient values (`daily_df`).
@@ -71,17 +73,23 @@ We then find the regression coefficients `temp_coeffs`, `wind_coeffs`, and `pres
 
 Using the found regression coefficients, we remove the linear relationship with meteorological parameter anomalies. The corrected PG is saved in `meteo_df["CorrectedField"]`.
 
-Finally, we construct Figure 4 using the prepared data in the same manner as was done for Figures 2 and 3.
+Finally, we construct Figure 4 using the prepared data in the same manner as was done for Figures 1.2 and 1.3.
 
 
 ## Script `3_WRF_T2_images.ipynb`
 
-This script calculates the seasonal variation of the 2m-level temperature (T2m) taken from climate modeling results (see article).
+This script calculates the seasonal variation of the 2m-level temperature (T2m) taken from climate modeling results (see article 1).
 
 In the script, temperature data averaged by longitude and by month are loaded (see data description below) from `WRF_T2_MONxLAT.npy`.
 
 Next, the temperature is averaged across latitude bands 20° S–20° N, 30° S–30° N, 40° S–40° N, and 50° S–50° N. The averaging takes into account the latitudinal area factor; degree cells at higher latitudes are summed with a diminishing coefficient. The results of the averaging (seasonal temperature variation in the specified latitude band) are displayed on a figure consisting of four panels.
 
+## Script `4_IP_simulations_temporal_images.ipynb`
+...
+
+## Script `5_IP_simulations_spatial_images.ipynb`
+...
+
 # Description of the data files
 
 * `WRF_T2_MONxLAT.npy` - a `numpy` array with the shape `(180, 12)`, containing montly averaged 2-meter level temperature (`T2m`) for each 1° latitude band across a full 360° longitude. `T2m` are calculated with the Weather Research and Forecasting model (WRF) version 4.3.