glm_extended.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f0935f0a-b796-4454-9080-349c98b7088a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: xarray in /srv/conda/envs/notebook/lib/python3.10/site-packages (2024.3.0)\n",
      "Requirement already satisfied: rasterio in /srv/conda/envs/notebook/lib/python3.10/site-packages (1.3.9)\n",
      "Collecting rio-cogeo\n",
      "  Downloading rio_cogeo-5.3.3-py3-none-any.whl.metadata (4.6 kB)\n",
      "Requirement already satisfied: numpy in /srv/conda/envs/notebook/lib/python3.10/site-packages (1.26.4)\n",
      "Requirement already satisfied: packaging>=22 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from xarray) (24.0)\n",
      "Requirement already satisfied: pandas>=1.5 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from xarray) (2.2.1)\n",
      "Requirement already satisfied: affine in /srv/conda/envs/notebook/lib/python3.10/site-packages (from rasterio) (2.4.0)\n",
      "Requirement already satisfied: attrs in /srv/conda/envs/notebook/lib/python3.10/site-packages (from rasterio) (23.2.0)\n",
      "Requirement already satisfied: certifi in /srv/conda/envs/notebook/lib/python3.10/site-packages (from rasterio) (2024.2.2)\n",
      "Requirement already satisfied: click>=4.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from rasterio) (8.1.7)\n",
      "Requirement already satisfied: cligj>=0.5 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from rasterio) (0.7.2)\n",
      "Requirement already satisfied: snuggs>=1.4.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from rasterio) (1.4.7)\n",
      "Requirement already satisfied: click-plugins in /srv/conda/envs/notebook/lib/python3.10/site-packages (from rasterio) (1.1.1)\n",
      "Requirement already satisfied: setuptools in /srv/conda/envs/notebook/lib/python3.10/site-packages (from rasterio) (69.2.0)\n",
      "Collecting morecantile<6.0,>=5.0 (from rio-cogeo)\n",
      "  Downloading morecantile-5.4.0-py3-none-any.whl.metadata (6.9 kB)\n",
      "Requirement already satisfied: pydantic~=2.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from rio-cogeo) (2.6.4)\n",
      "Requirement already satisfied: pyproj~=3.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from morecantile<6.0,>=5.0->rio-cogeo) (3.6.1)\n",
      "Requirement already satisfied: python-dateutil>=2.8.2 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from pandas>=1.5->xarray) (2.8.2)\n",
      "Requirement already satisfied: pytz>=2020.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from pandas>=1.5->xarray) (2024.1)\n",
      "Requirement already satisfied: tzdata>=2022.7 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from pandas>=1.5->xarray) (2024.1)\n",
      "Requirement already satisfied: annotated-types>=0.4.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from pydantic~=2.0->rio-cogeo) (0.6.0)\n",
      "Requirement already satisfied: pydantic-core==2.16.3 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from pydantic~=2.0->rio-cogeo) (2.16.3)\n",
      "Requirement already satisfied: typing-extensions>=4.6.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from pydantic~=2.0->rio-cogeo) (4.11.0)\n",
      "Requirement already satisfied: pyparsing>=2.1.6 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from snuggs>=1.4.1->rasterio) (3.1.2)\n",
      "Requirement already satisfied: six>=1.5 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from python-dateutil>=2.8.2->pandas>=1.5->xarray) (1.16.0)\n",
      "Downloading rio_cogeo-5.3.3-py3-none-any.whl (20 kB)\n",
      "Downloading morecantile-5.4.0-py3-none-any.whl (49 kB)\n",
      "\u001b[2K   \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m49.4/49.4 kB\u001b[0m \u001b[31m399.9 kB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m \u001b[36m0:00:01\u001b[0m\n",
      "\u001b[?25hInstalling collected packages: morecantile, rio-cogeo\n",
      "Successfully installed morecantile-5.4.0 rio-cogeo-5.3.3\n"
     ]
    }
   ],
   "source": [
    "!pip install xarray rasterio rio-cogeo numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "109505a2-55fc-4262-9676-ec8c5cd230a0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "copying lat, lon, data ...\n",
      "generating cog ...\n"
     ]
    },
    {
     "ename": "OSError",
     "evalue": "[Errno 39] Directory not empty: '/home/jovyan/lmatools-notebooks/data/GLM_input/'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mOSError\u001b[0m                                   Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[7], line 275\u001b[0m\n\u001b[1;32m    271\u001b[0m     preprocess_glm(file, variable, each_file_location)\u001b[38;5;241m.\u001b[39mconnector()\n\u001b[1;32m    274\u001b[0m shutil\u001b[38;5;241m.\u001b[39mrmtree(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mBASE_PATH\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124mdata/NALMA_input/\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m--> 275\u001b[0m \u001b[43mshutil\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrmtree\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43mf\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;132;43;01m{\u001b[39;49;00m\u001b[43mBASE_PATH\u001b[49m\u001b[38;5;132;43;01m}\u001b[39;49;00m\u001b[38;5;124;43mdata/GLM_input/\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/srv/conda/envs/notebook/lib/python3.10/shutil.py:731\u001b[0m, in \u001b[0;36mrmtree\u001b[0;34m(path, ignore_errors, onerror)\u001b[0m\n\u001b[1;32m    729\u001b[0m         os\u001b[38;5;241m.\u001b[39mrmdir(path)\n\u001b[1;32m    730\u001b[0m     \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mOSError\u001b[39;00m:\n\u001b[0;32m--> 731\u001b[0m         \u001b[43monerror\u001b[49m\u001b[43m(\u001b[49m\u001b[43mos\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrmdir\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpath\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msys\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mexc_info\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    732\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m    733\u001b[0m     \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m    734\u001b[0m         \u001b[38;5;66;03m# symlinks to directories are forbidden, see bug #1669\u001b[39;00m\n",
      "File \u001b[0;32m/srv/conda/envs/notebook/lib/python3.10/shutil.py:729\u001b[0m, in \u001b[0;36mrmtree\u001b[0;34m(path, ignore_errors, onerror)\u001b[0m\n\u001b[1;32m    727\u001b[0m     os\u001b[38;5;241m.\u001b[39mclose(fd)\n\u001b[1;32m    728\u001b[0m     fd_closed \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m\n\u001b[0;32m--> 729\u001b[0m     \u001b[43mos\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrmdir\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpath\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    730\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mOSError\u001b[39;00m:\n\u001b[1;32m    731\u001b[0m     onerror(os\u001b[38;5;241m.\u001b[39mrmdir, path, sys\u001b[38;5;241m.\u001b[39mexc_info())\n",
      "\u001b[0;31mOSError\u001b[0m: [Errno 39] Directory not empty: '/home/jovyan/lmatools-notebooks/data/GLM_input/'"
     ]
    }
   ],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import os\n",
    "import shutil\n",
    "import xarray as xa\n",
    "\n",
    "def read_nc_file(location):\n",
    "    return read_local_nc_file(location)\n",
    "\n",
    "def read_local_nc_file(location):\n",
    "    file = xa.open_dataset(location, engine=\"netcdf4\", decode_coords='all', decode_times=False)\n",
    "    return file\n",
    "\n",
    "def read_s3_nc_file(location):\n",
    "    pass\n",
    "\n",
    "def output_cog(file, directory_name, filename):\n",
    "    make_directory(directory_name)\n",
    "    file.rio.to_raster(f\"{directory_name}/{filename}.tif\", driver='COG')\n",
    "\n",
    "\n",
    "\n",
    "class preprocess_glm():\n",
    "    def __init__(self, file, variable, location):\n",
    "        self.file = file\n",
    "        self.variable = variable\n",
    "        self.location = location\n",
    "        self.r_eq, self.r_pol, self.H, self.lon0 = self.extract_variables()\n",
    "        \n",
    "    def connector(self):\n",
    "        lat, lon, data = self.copy_lat_lon_data()\n",
    "        file = self.generate_new_xarray(lat, lon, data)\n",
    "        self.generate_cog(file)\n",
    "    \n",
    "    def generate_cog(self, file):\n",
    "        print(\"generating cog ...\")\n",
    "        file = file[self.variable]\n",
    "        file = file.transpose('lat', 'lon')\n",
    "        file.rio.set_spatial_dims(x_dim='lon', y_dim='lat', inplace=True)\n",
    "        file.rio.crs\n",
    "        file.rio.set_crs('epsg:4326', inplace=True)\n",
    "\n",
    "        directory_name, filename = self.extract_file_directory_name(self.location)\n",
    "        output_cog(file, directory_name, filename)\n",
    "\n",
    "    def extract_file_directory_name(self, location):\n",
    "        filename = location.split('/')[-1]\n",
    "        filename = filename[:-3]\n",
    "        directory_name = \"GLM\"\n",
    "\n",
    "        return directory_name, filename\n",
    "\n",
    "    def copy_lat_lon_data(self):\n",
    "        print('copying lat, lon, data ...')\n",
    "        x_cpy = self.file.x.copy()\n",
    "        y_cpy = self.file.y.copy()\n",
    "\n",
    "        lat = np.zeros(self.file.x.size)\n",
    "        lon = np.zeros(self.file.x.size)\n",
    "        data = self.file[self.variable].data.copy()\n",
    "\n",
    "        for i in range(x_cpy.size):\n",
    "            lat[i], lon[i] = self.calc_lat_lon(x_cpy[i], y_cpy[i])\n",
    "\n",
    "        _lat = lat.copy()\n",
    "        _lon = lon.copy()\n",
    "        _data = self.file[self.variable].data.copy()\n",
    "\n",
    "        end_index = len(lon)\n",
    "        for index, value in np.ndenumerate(lon[::-1]):\n",
    "            end_index -= 1\n",
    "            if value == -999.0:\n",
    "                _data = np.delete(_data, end_index, axis=0)\n",
    "                _lon = np.delete(_lon, end_index)\n",
    "\n",
    "        end_index = len(lat)\n",
    "        for index, value in np.ndenumerate(lat[::-1]):\n",
    "            end_index -= 1\n",
    "            if value == -999.0:\n",
    "                _data = np.delete(_data, end_index, axis=1)\n",
    "                _lat = np.delete(_lat, end_index)\n",
    "\n",
    "        return _lat, _lon, _data\n",
    "    \n",
    "    def generate_new_xarray(self, lat, lon, data):\n",
    "        file = xa.Dataset(\n",
    "            data_vars={self.variable: ((\"lon\", \"lat\"), data)},\n",
    "            coords={\"lon\": lon, \"lat\": lat},\n",
    "            attrs={\"instrument_ID\": \"GLM-1\"}\n",
    "        )\n",
    "        return file\n",
    "\n",
    "    def extract_variables(self):\n",
    "        #radius of Earth at the equator in meters\n",
    "        r_eq = self.file['goes_imager_projection'].semi_major_axis\n",
    "        #radius of Earth at the pole in meters\n",
    "        r_pol = self.file['goes_imager_projection'].semi_minor_axis\n",
    "        #distance from satellite to earth center in meters\n",
    "        H = self.file['goes_imager_projection'].perspective_point_height+r_eq\n",
    "        #longitude of projection origin in units of radians\n",
    "        lon0 = self.file['goes_imager_projection'].longitude_of_projection_origin * (math.pi/180.)\n",
    "        return r_eq, r_pol, H, lon0\n",
    "\n",
    "    def calc_lat_lon(self, xx,yy):\n",
    "        #intermediate calculations\n",
    "        a = (math.sin(xx))**2 + (math.cos(xx))**2 * ((math.cos(yy))**2 + self.r_eq**2/self.r_pol**2*(math.sin(yy))**2)\n",
    "        b = -2 * self.H * math.cos(xx) * math.cos(yy)\n",
    "        c = self.H**2 - self.r_eq**2\n",
    "\n",
    "        if b**2-4*a*c < 0:\n",
    "            lat = -999.\n",
    "            lon = -999.\n",
    "        else:\n",
    "            #distance from the satellite to point P\n",
    "            r_s = (-b-math.sqrt(b**2-4*a*c))/2/a\n",
    "\n",
    "            #derived using satellite location and earth geometry s_x, s_y, s_z\n",
    "            s_x = r_s * math.cos(xx) * math.cos(yy)\n",
    "            s_y = -r_s * math.sin(xx)\n",
    "            s_z = r_s * math.cos(xx) * math.sin(yy)\n",
    "\n",
    "            #geodetic latitude and longitude for point P\n",
    "            #math.atan(x) return the arc tangent of x, in radians. The result is between -pi/2 and pi/2\n",
    "            lat = math.atan(((self.r_eq**2)/(self.r_pol**2)) * (s_z/math.sqrt((self.H-s_x)**2 + s_y**2))) #radians\n",
    "            lon = self.lon0 - math.atan(s_y/(self.H-s_x)) #radians\n",
    "\n",
    "            #convert radians to degrees\n",
    "            lat = lat * 180./math.pi #degree\n",
    "            lon = lon * 180./math.pi #degree\n",
    "\n",
    "        return lat, lon\n",
    "\n",
    "\n",
    "class preprocess_nalma():\n",
    "    def __init__(self, file, variable, location):\n",
    "        self.latitude = 'latitude'\n",
    "        self.longitude = 'longitude'\n",
    "        self.variable = variable\n",
    "        print(f\">>>>>>> Variable: {self.variable}\")\n",
    "        self.file = file\n",
    "        self.location = location\n",
    "    \n",
    "    def connector(self):\n",
    "        total_grids = len(self.file[self.variable].data)\n",
    "        print(f\"total_grids: {total_grids}\")\n",
    "        for grid_number in range(total_grids):\n",
    "            lat, lon, data = self.copy_lat_lon_data(self.file, number=grid_number)\n",
    "            if self.check_flashes(data):\n",
    "                lat, lon, data = self.delete_row_col(lat, lon, data)\n",
    "                lat, lon, data = self.flip_lat_data(lat, lon, data)\n",
    "                new_file = self.make_new_xarray(lat, lon, data)\n",
    "                self.generate_cog(new_file, grid_number)\n",
    "            else:\n",
    "                directory_name, filename = self.extract_file_directory_name(self.location, grid_number)\n",
    "                print(f\"### ### ### No flashes on file: {filename}\")\n",
    "                continue\n",
    "\n",
    "\n",
    "    def copy_lat_lon_data(self, file, number=0):\n",
    "        lat = file[self.latitude].data.copy()\n",
    "        lon = file[self.longitude].data.copy()\n",
    "        data = file[self.variable].data[number].copy()\n",
    "        return lat, lon, data\n",
    "    \n",
    "    def make_new_xarray(self, lat, lon, data):\n",
    "        file = xa.Dataset(\n",
    "            {\n",
    "                self.variable: ([self.longitude, self.latitude], data),\n",
    "            },\n",
    "            coords={\n",
    "                self.longitude: ([self.longitude], lon),\n",
    "                self.latitude: ([self.latitude], lat),\n",
    "            },\n",
    "        )\n",
    "        return file\n",
    "    \n",
    "    def extract_file_directory_name(self, location, grid_number):\n",
    "        #extract original filename from location\n",
    "        filename = location.split('/')[-1][:-3]\n",
    "        directory_name =  filename[0:21]\n",
    "        return directory_name, f\"{filename}_{grid_number}\"\n",
    "    \n",
    "    def generate_cog(self, file, grid_number):\n",
    "        file = file[self.variable]\n",
    "        file = file.transpose(self.latitude, self.longitude)\n",
    "        file.rio.set_crs('epsg:4326')\n",
    "        file.rio.set_spatial_dims(x_dim=self.longitude, y_dim=self.latitude, inplace=True)\n",
    "        directory_name, filename = self.extract_file_directory_name(self.location, grid_number)\n",
    "        output_cog(file, directory_name, filename)\n",
    "\n",
    "    def delete_row_col(self, lat, lon, data):\n",
    "        #finding row with all zeroes    \n",
    "        zero_row = np.where(np.all(data == 0, axis=1))[0]\n",
    "        #removing lon with all zeroes\n",
    "        lon = np.delete(lon, zero_row)\n",
    "        # Find rows with non-zero elements\n",
    "        non_zero_rows = np.any(data != 0, axis=1)\n",
    "        # Filter the array based on non-zero rows\n",
    "        data = data[non_zero_rows]\n",
    "\n",
    "        #finding columns with all zeroes\n",
    "        zero_columns = np.where(np.all(data == 0, axis=0))[0]\n",
    "        #removing lat with all zeroes\n",
    "        lat = np.delete(lat, zero_columns)\n",
    "        # Find rows with non-zero elements\n",
    "        non_zero_rows = np.any(data != 0, axis=1)\n",
    "        # Find columns with non-zero elements\n",
    "        non_zero_columns = np.any(data != 0, axis=0)\n",
    "        # Filter the array based on non-zero rows and columns\n",
    "        data = data[non_zero_rows][:, non_zero_columns]\n",
    "        return lat, lon, data\n",
    "    \n",
    "    def flip_lat_data(self, lat, lon, data):\n",
    "        lat = lat[::-1]\n",
    "        data = np.flip(data, axis=1)\n",
    "        return lat, lon, data\n",
    "\n",
    "    def check_flashes(self, data):\n",
    "        count = 0\n",
    "        for row in data:\n",
    "            for value in row:\n",
    "                if value != 0:\n",
    "                    count += 1\n",
    "        if count == 0:\n",
    "            return False\n",
    "        return True\n",
    "\n",
    "\n",
    "\n",
    "def make_directory(path):\n",
    "    # path = os.path.join(base_path, directory_name)\n",
    "    if os.path.exists(path) == False:\n",
    "        os.mkdir(path)\n",
    "        return True\n",
    "    return False\n",
    "\n",
    "def get_all_nalma_files_in_a_directory(directory_path):\n",
    "    files = []\n",
    "    for f in os.listdir(directory_path):\n",
    "        if os.path.isfile(os.path.join(directory_path, f)) and f[0:5] == 'NALMA':\n",
    "            files.append(os.path.join(directory_path, f))\n",
    "    return files\n",
    "\n",
    "def get_all_files_in_a_directory(directory_path):\n",
    "    files = []\n",
    "    for f in os.listdir(directory_path):\n",
    "        if os.path.isfile(os.path.join(directory_path, f)):\n",
    "            files.append(os.path.join(directory_path, f))\n",
    "    return files\n",
    "\n",
    "\n",
    "\n",
    "BASE_PATH = '/home/jovyan/lmatools-notebooks/'\n",
    "\n",
    "\n",
    "#NALMA\n",
    "# all_local_nalma_files = get_all_nalma_files_in_a_directory(f\"{BASE_PATH}data/NALMA_input/\")\n",
    "# for each_file_location in all_local_nalma_files:\n",
    "#     #read NALMA netcdf4 file\n",
    "#     file = read_nc_file(each_file_location)\n",
    "#     #define variable that needs to be converted to geotiff \n",
    "#     variable = \"flash_extent\"\n",
    "#     #preprocess it for COG conversion\n",
    "#     preprocess_nalma(file, variable, each_file_location).connector()\n",
    "\n",
    "#GLM\n",
    "all_local_glm_files = get_all_files_in_a_directory(f\"{BASE_PATH}data/GLM_input/\")\n",
    "for each_file_location in all_local_glm_files:\n",
    "    file = read_nc_file(each_file_location)\n",
    "    variable = 'Flash_extent_density'\n",
    "    preprocess_glm(file, variable, each_file_location).connector()\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf173cf9-2487-424f-9d2d-de01c12d0165",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}