MODTRAN4_prep.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1dfaef56-c6d5-49b0-b921-26ff4d4c7f31",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%pip install prettytable\n",
    "%pip install metpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d3138918-3a5c-45f3-a660-732991b42712",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# importing required library\n",
    "from prettytable import PrettyTable\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import os\n",
    "import metpy\n",
    "from metpy.units import units\n",
    "import numpy as np\n",
    "from pathlib import Path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "149dfd0c-62f8-42ae-979e-a6d71678c404",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Get the file names\n",
    "basepath = Path('/home/jovyan/Landsat_SST_algorithm')\n",
    "atmpath = basepath / 'Data/ERA5_atmprofiles/Atmospheres/'\n",
    "sstpath = basepath / 'Data/ERA5_atmprofiles/SSTs/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "081e212b-81c2-4835-b17f-82bf23e1b57f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def open_ERA5(monthfile,path):\n",
    "    # Open all files from one month into dataframes and concatenate\n",
    "    ds = xr.open_dataset(path / monthfile)\n",
    "    ds = ds.drop_sel(longitude=np.delete(ds.longitude, np.arange(0, ds.longitude.size, 4)))\n",
    "    df = ds.to_dataframe()\n",
    "    df = df.reset_index()\n",
    "    \n",
    "    # print(ds.dims['longitude']*ds.dims['latitude']*ds.dims['time'])\n",
    "    \n",
    "    # Needs to be all levels from one profile, then next profile, then concatenate\n",
    "    df = df.sort_values(by=['longitude', 'latitude','time'])\n",
    "    \n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c89db371-4463-4101-857d-d38874c1de85",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create MODTRAN atm correction input files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "7cf69402-e684-4e49-92d8-d84e161e6498",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['era5_daily_202001.nc', 'era5_daily_202101.nc']\n",
      "['era5_SST_202101.nc', 'era5_SST_202001.nc']\n",
      "Max SST: 2.794342041015625, Min: -1.6900634765625, Size: 1629\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_117/2498949025.py:48: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  is_sst.loc[:,'sst'] = is_sst['sst'] - 273.15\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['era5_daily_202102.nc', 'era5_daily_202002.nc']\n",
      "['era5_SST_202002.nc', 'era5_SST_202102.nc']\n",
      "Max SST: 2.923004150390625, Min: -1.743988037109375, Size: 1627\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_117/2498949025.py:48: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  is_sst.loc[:,'sst'] = is_sst['sst'] - 273.15\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['era5_daily_202003.nc', 'era5_daily_202103.nc']\n",
      "['era5_SST_202103.nc', 'era5_SST_202003.nc']\n",
      "Max SST: 3.113189697265625, Min: -1.85552978515625, Size: 1630\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_117/2498949025.py:48: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  is_sst.loc[:,'sst'] = is_sst['sst'] - 273.15\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['era5_daily_202104.nc', 'era5_daily_202004.nc']\n",
      "['era5_SST_202004.nc', 'era5_SST_202104.nc']\n",
      "Max SST: 1.909576416015625, Min: -1.690521240234375, Size: 1632\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_117/2498949025.py:48: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  is_sst.loc[:,'sst'] = is_sst['sst'] - 273.15\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['era5_daily_202005.nc', 'era5_daily_202105.nc']\n",
      "['era5_SST_202105.nc', 'era5_SST_202005.nc']\n",
      "Max SST: 1.647613525390625, Min: -1.690277099609375, Size: 1632\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_117/2498949025.py:48: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  is_sst.loc[:,'sst'] = is_sst['sst'] - 273.15\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['era5_daily_202106.nc', 'era5_daily_202006.nc']\n",
      "['era5_SST_202006.nc', 'era5_SST_202106.nc']\n",
      "Max SST: 1.391754150390625, Min: -1.692962646484375, Size: 1632\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_117/2498949025.py:48: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  is_sst.loc[:,'sst'] = is_sst['sst'] - 273.15\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['era5_daily_202007.nc', 'era5_daily_202107.nc']\n",
      "['era5_SST_202107.nc', 'era5_SST_202007.nc']\n",
      "Max SST: 1.225738525390625, Min: -1.7010498046875, Size: 1632\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_117/2498949025.py:48: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  is_sst.loc[:,'sst'] = is_sst['sst'] - 273.15\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['era5_daily_202108.nc', 'era5_daily_202008.nc']\n",
      "['era5_SST_202008.nc', 'era5_SST_202108.nc']\n",
      "Max SST: 0.786529541015625, Min: -1.693695068359375, Size: 1632\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_117/2498949025.py:48: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  is_sst.loc[:,'sst'] = is_sst['sst'] - 273.15\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['era5_daily_202009.nc', 'era5_daily_202109.nc']\n",
      "['era5_SST_202109.nc', 'era5_SST_202009.nc']\n",
      "Max SST: 0.901275634765625, Min: -1.6927490234375, Size: 1632\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_117/2498949025.py:48: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  is_sst.loc[:,'sst'] = is_sst['sst'] - 273.15\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['era5_daily_202010.nc', 'era5_daily_202110.nc']\n",
      "['era5_SST_202110.nc', 'era5_SST_202010.nc']\n",
      "Max SST: 0.891998291015625, Min: -1.709320068359375, Size: 1632\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_117/2498949025.py:48: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  is_sst.loc[:,'sst'] = is_sst['sst'] - 273.15\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['era5_daily_202111.nc', 'era5_daily_202011.nc']\n",
      "['era5_SST_202011.nc', 'era5_SST_202111.nc']\n",
      "Max SST: 0.677398681640625, Min: -1.7630615234375, Size: 1631\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_117/2498949025.py:48: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  is_sst.loc[:,'sst'] = is_sst['sst'] - 273.15\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['era5_daily_202012.nc', 'era5_daily_202112.nc']\n",
      "['era5_SST_202112.nc', 'era5_SST_202012.nc']\n",
      "Max SST: 1.318267822265625, Min: -1.691741943359375, Size: 1629\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_117/2498949025.py:48: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  is_sst.loc[:,'sst'] = is_sst['sst'] - 273.15\n"
     ]
    }
   ],
   "source": [
    "months = ['01','02','03','04','05','06','07','08','09','10','11','12']\n",
    "\n",
    "dir_list = os.listdir(atmpath)\n",
    "sst_list = os.listdir(sstpath)\n",
    "\n",
    "for i in months:\n",
    "    \n",
    "    # Prep atmospheric profiles\n",
    "    # Find all files for a month\n",
    "    monthfiles = [file for file in dir_list if file.endswith(f'{i}.nc')]\n",
    "    print (monthfiles)\n",
    "    \n",
    "    dxs = []\n",
    "    \n",
    "    # Open each file and concatenate all files together\n",
    "    for monthfile in monthfiles:\n",
    "        dx = open_ERA5(monthfile,atmpath)\n",
    "        dxs.append(dx)\n",
    "\n",
    "    df = pd.concat(dxs, ignore_index=True)\n",
    "    \n",
    "    # Choose specific days and hours to thin the data - we chose Day 1 00h, Day 7 12h, Day 15 6h, Day 23 18h\n",
    "    df = df[(df.time.dt.day==1)&(df.time.dt.hour==0)|(df.time.dt.day==7)&(df.time.dt.hour==12)|(df.time.dt.day==15)&(df.time.dt.hour==6)|(df.time.dt.day==23)&(df.time.dt.hour==18)]\n",
    "    # print(f'Atm number: {df[df.level==1].shape[0]}')\n",
    "    \n",
    "    # Prep SSTs\n",
    "    # Find all files for a month\n",
    "    monthfiles = [file for file in sst_list if file.endswith(f'{i}.nc')]\n",
    "    print (monthfiles)\n",
    "    \n",
    "    dxs = []\n",
    "    \n",
    "    # Open each file and concatenate all files together\n",
    "    for monthfile in monthfiles:\n",
    "        dx = open_ERA5(monthfile,sstpath)\n",
    "        dxs.append(dx)\n",
    "\n",
    "    sstf = pd.concat(dxs, ignore_index=True)\n",
    "    \n",
    "    # Choose specific days and hours to thin the data - Day 1 00h, Day 7 12h, Day 15 6h, Day 23 18h\n",
    "    sstf = sstf[(sstf.time.dt.day==1)&(sstf.time.dt.hour==0)|(sstf.time.dt.day==7)&(sstf.time.dt.hour==12)|(sstf.time.dt.day==15)&(sstf.time.dt.hour==6)|(sstf.time.dt.day==23)&(sstf.time.dt.hour==18)]\n",
    "    # print(f'SST number: {sstf.shape[0]}')\n",
    "    \n",
    "    # Remove measurements not over the ocean\n",
    "    is_sst = sstf[sstf['sst'].notna()]\n",
    "    \n",
    "    # Convert to Celsius and remove measurements over frozen ocean\n",
    "    is_sst.loc[:,'sst'] = is_sst['sst'] - 273.15\n",
    "    is_sst = is_sst[is_sst['sst'] > -1.9]\n",
    "    \n",
    "    # Trim atmospheric profiles based on remaining SST measurements and merge the dataframes\n",
    "    good_df = pd.merge(df,is_sst,how='inner',left_on=['longitude','latitude','time'],right_on=['longitude','latitude','time'])\n",
    "\n",
    "    print(f'Max SST: {good_df.sst.max()}, Min: {good_df.sst.min()}, Size: {is_sst.shape[0]}')\n",
    "    \n",
    "    # Ensure the merge produced the expected output\n",
    "    if is_sst.shape[0] != good_df[good_df['level']==1].shape[0]: \n",
    "        print ('SST and atm profiles do not match') \n",
    "        continue\n",
    "\n",
    "    # Add units to level *slow\n",
    "    good_df ['hPa'] = good_df['level'].apply(lambda x: x*units.hectopascal)\n",
    "\n",
    "    # Convert hPa to height in km *takes a long time\n",
    "    good_df['ht[km]'] = good_df['hPa'].apply(metpy.calc.pressure_to_height_std)\n",
    "    good_df['ht'] = good_df['ht[km]'].apply(lambda x: x.magnitude)\n",
    "    \n",
    "    # Prep for output and save to file\n",
    "    good_df = good_df.reset_index(drop=True)\n",
    "    good_df['ht'] = np.around(good_df['ht'],1)\n",
    "    good_df['t'] = np.around(good_df['t'],1)\n",
    "    good_df['q'] = np.around(good_df['q'],7)\n",
    "     \n",
    "    # Save atmopheric profiles and SSTs\n",
    "    outFile = basepath / f'Data/AtmCorrection/modtran_atmprofiles_{i}_20230823.txt'\n",
    "    good_df[['ht','level','t','q']].to_csv(outFile,sep='\\t',index=False,header=False, encoding='ascii')\n",
    "    \n",
    "    outFile = basepath / f'Data/AtmCorrection/modtran_sstprofiles_{i}_20230823.txt'\n",
    "    is_sst[['sst']].to_csv(outFile,sep='\\t',index=False,header=False, encoding='ascii')"
   ]
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}