AVG amplitude total.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7528690",
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "import numpy as np\n",
    "import scipy\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import signal\n",
    "import pandas as pd\n",
    "import scipy.sparse as sparse\n",
    "import os\n",
    "\n",
    "indir_name =  \"C:/CSES/file/\"\n",
    "\n",
    "outdir_name = \"C:/CSES/plot/\"\n",
    "\n",
    "ext = ('.h5')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7cafe9ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "def getData(data):\n",
    "    Bw = 51200 * 2 / 1024\n",
    "    data = pd.DataFrame(data)\n",
    "\n",
    "    matrix = []\n",
    "    for i in range(len(data)):\n",
    "        matrix.append(data.iloc[i])\n",
    "    matrix = np.array(matrix)\n",
    "    data_t = np.empty(shape=(matrix.shape[0], matrix.shape[1]))\n",
    "    for i in range(0, matrix.shape[0]):\n",
    "        meanX_b = np.mean(matrix[i])\n",
    "        data_t[i] = matrix[i] - meanX_b\n",
    "\n",
    "    M_b = data_t.shape[1]\n",
    "    hamming_b = signal.get_window(\"hamming\", M_b)\n",
    "    FFT_low = np.array([scipy.fft.fft(data_t[i] * hamming_b) for i in range(0, data_t.shape[0])])\n",
    "    out = np.abs(FFT_low.T[:1024]) ** 2\n",
    "    outX_b = 400 + 20 * np.log10(out / Bw)\n",
    "\n",
    "    return outX_b\n",
    "\n",
    "\n",
    "def powerSpectrum(pow):\n",
    "    powerX = 400 + 20 * np.log10(pow)\n",
    "\n",
    "    return powerX\n",
    "\n",
    "def tempArray(data):\n",
    "    temp_df_x = []\n",
    "    for i in range(len(df)):\n",
    "        if (df.iloc[i].WORKMODE == 2):\n",
    "            temp_df_x.append(df[data].iloc[i])\n",
    "        else:\n",
    "            temp_df_x.append(np.empty(np.array(df[data].iloc[i]).shape))\n",
    "            temp_df_x[i][:] = np.NaN\n",
    "    temp_df_x = np.array(temp_df_x)\n",
    "    return temp_df_x\n",
    "\n",
    "\n",
    "\n",
    "def Amplitude(arr):\n",
    "\n",
    "    mask = ~np.isnan(arr[row])\n",
    "    dataX = arr[row][mask]\n",
    "    for i in range(0, arr.shape[0]):\n",
    "        if i != row:\n",
    "            mask = ~np.isnan(arr[i])\n",
    "            arr[i][mask] = np.nan\n",
    "    return arr\n",
    "\n",
    "def dataframeBurstBuild(arrX,arrY,arrZ):\n",
    " \n",
    "    df_burst = pd.DataFrame(list(zip(arrX[row, :], arrY[row, :], arrZ[row, :], GEO_LAT, GEO_LON)),\n",
    "                             columns=['arrXb',\n",
    "                                      'arrYb',\n",
    "                                      'arrZb', 'GEO_LAT',\n",
    "                                      'GEO_LON'])\n",
    "     \n",
    "\n",
    "    df_burst['DATE2'] = DATE2\n",
    "     \n",
    "    arraysXb.append([df_burst['arrXb'].to_numpy(),\n",
    "                      OrbitNumber + '_' + DATE[0]])\n",
    "\n",
    "  \n",
    "\n",
    "    arraysYb.append([df_burst['arrYb'].to_numpy(),\n",
    "                      OrbitNumber + '_' + DATE[0]])\n",
    "\n",
    "    \n",
    "\n",
    "    arraysZb.append([df_burst['arrZb'].to_numpy(),\n",
    "                      OrbitNumber + '_' + DATE[0]])\n",
    "     \n",
    "    return arraysXb,arraysYb,arraysZb,df_burst\n",
    "\n",
    "def dataframeBuild(arrX,arrY,arrZ):\n",
    "\n",
    "    df['arrX'] = arrX.T[row, :].tolist()\n",
    "    df['arrY'] = arrY.T[row, :].tolist()\n",
    "    df['arrZ'] = arrZ.T[row, :].tolist()\n",
    "\n",
    "    arraysX.append([df['arrX'].to_numpy(), OrbitNumber + '_' + DATE[0]])\n",
    "    \n",
    "    arraysY.append([df['arrY'].to_numpy(), OrbitNumber + '_' + DATE[0]])\n",
    "   \n",
    "    arraysZ.append([df['arrZ'].to_numpy(), OrbitNumber + '_' + DATE[0]])\n",
    "    return arraysX,arraysY,arraysZ,df\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cced81d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "      \n",
    "\n",
    "def readFile(f):\n",
    "    UTC_TIME = f[\"UTC_TIME\"][()][:, 0]\n",
    "    GEO_LAT = f[\"GEO_LAT\"][()][:, 0]\n",
    "    GEO_LON = f[\"GEO_LON\"][()][:, 0]\n",
    "    #ALT = f[\"ALTITUDE\"][()][:, 0]\n",
    "    Workmode = f[\"WORKMODE\"][()][:, 0]\n",
    "    #MAG_LAT = f[\"MAG_LAT\"][()][:, 0]\n",
    "    #MAG_LON = f[\"MAG_LON\"][()][:, 0]\n",
    "    #VERSE_TIME = f[\"VERSE_TIME\"][()][:, 0]\n",
    "    A131_W = f[\"A131_W\"][()]\n",
    "    A132_W = f[\"A132_W\"][()]\n",
    "    A133_W = f[\"A133_W\"][()]\n",
    "    A131_P = f[\"A131_P\"][()]\n",
    "    A132_P = f[\"A132_P\"][()]\n",
    "    A133_P = f[\"A133_P\"][()]\n",
    "    columns = list(f.keys())\n",
    "    df = pd.DataFrame([])\n",
    "    for column in columns:\n",
    "        try:\n",
    "            data = np.array(f[column])\n",
    "            if data.shape[1] == 1:\n",
    "                df[column] = data.flatten()\n",
    "            elif column.endswith('_P'):  \n",
    "                mat = sparse.coo_matrix(data, shape=data.shape)\n",
    "                df[column] = mat.toarray().tolist()\n",
    "            elif column == \"A131_W\":\n",
    "                selected_data = np.array(data[0:len(Workmode), :])\n",
    "                mat = sparse.coo_matrix(selected_data, shape=selected_data.shape)\n",
    "                df[column] = mat.toarray().tolist()\n",
    "            elif column == \"A132_W\":\n",
    "                selected_data = np.array(data[0:len(Workmode), :])\n",
    "                mat = sparse.coo_matrix(selected_data, shape=selected_data.shape)\n",
    "                df[column] = mat.toarray().tolist()\n",
    "            elif column == \"A133_W\":\n",
    "                selected_data = np.array(data[0:len(Workmode), :])\n",
    "                mat = sparse.coo_matrix(selected_data, shape=selected_data.shape)\n",
    "                df[column] = mat.toarray().tolist()\n",
    "            else:\n",
    "                print(column + ' skipped')\n",
    "        except Exception as e:\n",
    "            pass\n",
    "    S_burst = df[df.WORKMODE == 2]\n",
    "    df['DATE_TIME'] = pd.to_datetime(df.UTC_TIME, format='%Y%m%d%H%M%S%f')\n",
    "    DATE = df.DATE_TIME.map(lambda x: x.strftime('%Y-%m-%d'))\n",
    "    DATE2 = df.DATE_TIME.map(lambda x: x.strftime('%Y-%m'))\n",
    "    \n",
    " \n",
    "    return GEO_LAT,GEO_LON, A131_W, A132_W, A133_W, A131_P, A132_P, A133_P, DATE, df, S_burst,DATE2\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "582edea5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def avgtot(array,ymin,ymax, latmin,latmax,xmin,xmax,title,namefile):\n",
    "    vals_mean = []\n",
    "    vals_std = []\n",
    "    x_scale = []\n",
    "    months=pd.period_range(start='2018-08-01', end='2020-08-01', freq='M').to_series().astype(str)\n",
    "             \n",
    "    cmap = plt.cm.get_cmap('rainbow', 12)\n",
    "    c = cmap(np.linspace(0, 1, 12))\n",
    "    c2 = np.concatenate((c,c))\n",
    "    color = np.concatenate((c2,c))\n",
    "    \n",
    "    plt.figure(figsize=(60, 45))\n",
    "    arraysXdf = pd.DataFrame(array)\n",
    "    plt.xlabel('GEO_LAT', fontsize=40)\n",
    "    plt.ylabel('Amplitude [dB]', fontsize=40)\n",
    "    plt.xticks(fontsize=30)\n",
    "    plt.yticks(np.arange(-1000, 1000, step=50),fontsize=30)\n",
    "   \n",
    "    plt.title(title, fontsize=60)\n",
    "   \n",
    "    for i, date in enumerate(months.unique()):\n",
    "\n",
    "        arraysXdf_sel = arraysXdf[arraysXdf[1].str.contains(date)]\n",
    "\n",
    "        arraysX_table = pd.DataFrame()\n",
    "        for j in arraysXdf_sel.index:\n",
    "            row = pd.DataFrame(arraysXdf_sel[0][j]).transpose()\n",
    "            arraysX_table = pd.concat([arraysX_table, row], axis=0)\n",
    "\n",
    "        arraysX_table.reset_index()\n",
    "\n",
    "        vals_mean = arraysX_table.mean().transpose().to_numpy() + (i + 1) * (-50)\n",
    "        \n",
    "        vals_std = arraysX_table.std().transpose().to_numpy()\n",
    "        \n",
    "        x_scale = (np.arange(len(vals_mean)) / (len(vals_mean) - 1)) * (latmax - latmin) + latmin\n",
    "        plt.plot(x_scale, vals_mean, label=date,color=color[i], linewidth=7.0)\n",
    "        plt.fill_between(x_scale, vals_mean - vals_std, vals_mean + vals_std, color=color[i], alpha=.2)\n",
    "        plt.ylim(ymin, ymax)\n",
    "        plt.xlim(xmin, xmax)\n",
    "\n",
    "        \n",
    "        plt.legend(bbox_to_anchor=(1, 0.5), loc='center left', fontsize=40)\n",
    "        plt.savefig( namefile, bbox_inches='tight')\n",
    "\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2ded8e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "arraysX = []\n",
    "\n",
    "arraysY = []\n",
    "\n",
    "arraysZ = []\n",
    "\n",
    "arraysXb = []\n",
    "\n",
    "arraysYb = []\n",
    "\n",
    "arraysZb = []\n",
    "\n",
    "max_global_meanX = 1000\n",
    "min_global_meanX= -1000\n",
    "max_global_meanY = 1000\n",
    "min_global_meanY= -1000\n",
    "max_global_meanZ = 1000\n",
    "min_global_meanZ= -1000\n",
    "max_global_meanXb = 1000\n",
    "min_global_meanXb= -1000\n",
    "max_global_meanYb = 1000\n",
    "min_global_meanYb= -1000\n",
    "max_global_meanZb = 1000\n",
    "min_global_meanZb= -1000\n",
    "\n",
    "dir_name = \"\"\n",
    "file_name = dir_name + indir_name\n",
    "ext = ('.h5')\n",
    "\n",
    "sampleFreq=51200\n",
    "nRow=1024\n",
    "maxFreq=sampleFreq/2\n",
    "freqRow=maxFreq/nRow\n",
    "freq_array=np.arange(100,2600,100)\n",
    "for i in freq_array:\n",
    "    for path, dirc, files in os.walk(file_name):\n",
    "        for name in files:\n",
    "            if name.endswith('.h5'):\n",
    "                    OrbitNumber = name.split(\"_\")[6]\n",
    "                    with h5py.File(str(file_name) + str(name), \"r\") as f:\n",
    "\n",
    "                        GEO_LAT, GEO_LON, A131_W, A132_W, A133_W, A131_P, A132_P, A133_P, DATE, df, S_burst, DATE2 = readFile(f)\n",
    "                        f.close()\n",
    "                \n",
    "                    row=int(i/freqRow)\n",
    "                    powerX = powerSpectrum(A131_P)\n",
    "                    powerY = powerSpectrum(A132_P)\n",
    "                    powerZ = powerSpectrum(A133_P)\n",
    "        \n",
    "                    temp_df_x=tempArray('A131_W')\n",
    "                    outX_b = getData(temp_df_x)\n",
    "\n",
    "                    temp_df_y=tempArray('A132_W')\n",
    "\n",
    "                    outY_b = getData(temp_df_y)\n",
    "                \n",
    "                    temp_df_z=tempArray('A133_W')\n",
    "\n",
    "                    outZ_b = getData(temp_df_z)\n",
    "                \n",
    "\n",
    "                    outX_b = Amplitude(outX_b)\n",
    "                    outY_b = Amplitude(outY_b)\n",
    "                    outZ_b = Amplitude(outZ_b)\n",
    "            \n",
    "                    arraysX,arraysY,arraysZ,df=dataframeBuild(powerX,powerY,powerZ)\n",
    "     \n",
    "                    arraysXb,arraysYb,arraysZb,df_burst=dataframeBurstBuild(outX_b,outY_b,outZ_b)\n",
    "    \n",
    " \n",
    "\n",
    "    avgtot(arraysX,min_global_meanX,max_global_meanX,GEO_LAT.min(),GEO_LAT.max(),GEO_LAT.min(),GEO_LAT.max(),f'Amplitude EFDX Avg {i}Hz_ all orbits of every month' ,outdir_name +f'Amplitude EFDX Avg {i}Hz_ all orbits of every month.jpg')\n",
    "    avgtot(arraysY,min_global_meanY,max_global_meanY,GEO_LAT.min(),GEO_LAT.max(),GEO_LAT.min(),GEO_LAT.max(),f'Amplitude EFDY Avg {i}Hz all orbits of every month' ,outdir_name +f'Amplitude EFDY Avg {i}Hz all orbits of every month.jpg')\n",
    "    avgtot(arraysZ,min_global_meanZ,max_global_meanZ,GEO_LAT.min(),GEO_LAT.max(),GEO_LAT.min(),GEO_LAT.max(),f'Amplitude EFDZ Avg {i}Hz all orbits of every month' ,outdir_name +f'Amplitude EFDZ Avg {i}Hz all orbits of every month.jpg')\n",
    "    avgtot(arraysXb,min_global_meanXb,max_global_meanXb,GEO_LAT.min(),GEO_LAT.max(), S_burst.GEO_LAT.min(), S_burst.GEO_LAT.max(),f'Amplitude burst EFDX Avg {i}Hz all orbits of every month' ,outdir_name +f'Amplitude burst EFDX Avg {i}Hz all orbits of every month.jpg')\n",
    "    avgtot(arraysYb,min_global_meanYb,max_global_meanYb,GEO_LAT.min(),GEO_LAT.max(), S_burst.GEO_LAT.min(), S_burst.GEO_LAT.max(),f'Amplitude burst EFDY Avg {i}Hz all orbits of every month' ,outdir_name +f'Amplitude burst EFDY Avg {i}Hz all orbits of every month.jpg')\n",
    "    avgtot(arraysZb,min_global_meanZb,max_global_meanZb,GEO_LAT.min(),GEO_LAT.max(), S_burst.GEO_LAT.min(), S_burst.GEO_LAT.max(),f'Amplitude burst EFDZ Avg {i}Hz all orbits of every month' ,outdir_name +f'Amplitude burst EFDZ Avg {i}Hz all orbits of every month.jpg')\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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}