load_process_imagery.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "97f70bdd-945f-4fb6-83ea-7acc89c475e9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from glob import glob\n",
    "import urllib\n",
    "import earthaccess\n",
    "from datatree import open_datatree\n",
    "from PIL import Image\n",
    "import xarray as xr\n",
    "from satpy.scene import Scene\n",
    "import os\n",
    "import re\n",
    "from pyresample import create_area_def\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from datetime import datetime, timedelta\n",
    "from satpy.writers import get_enhanced_image\n",
    "from trollimage.xrimage import XRImage\n",
    "from math import isnan\n",
    "from PIL import Image\n",
    "from pyproj import Transformer\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a51a125-6f94-4479-b4d8-21ecdae930cb",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Access and download netCDF files\n",
    "Data are downloaded from NASA's [EarthData server](https://search.earthdata.nasa.gov/search) to the local directory. This is accomplished using the [earthaccess](https://nsidc.github.io/earthaccess/) package. More information on using EarthData on the cloud is available in the [NASA Earthdata Cloud Cookbook](https://nasa-openscapes.github.io/earthdata-cloud-cookbook/).\n",
    "\n",
    "1. Need to create and account and login to EarthData using earthaccess.login()\n",
    "1. Specify which data products to download using NASA's \"short names\"\n",
    "   * Products with short names beginning \"VNP\" are from Suomi-NPP\n",
    "   * Products with short names beginning \"VJ\" are from NOAA-20\n",
    "1. For the Level 1b products, geolocation files accompanying the imagery files are also required\n",
    "   * Data files have short names containing \"02\"\n",
    "   * Geolocation files have short names containing \"03\"\n",
    "   \n",
    "Note that the bounding_box (geographic coordinates of interest) and temporal_limits (time range of interest) should be definted before the download_files() function is called. If a data file contains any data from within the bounding_box, the entire data file will be downloaded.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "4341009c-6984-4ef7-9b80-ef43e78d0a27",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#Function to query the NASA EarthData archive and download netCDF files into the directories ~/data/level1b/ and ~/data/level2 (if using L2 cloud masks)\n",
    "def download_files(data_path, bounding_box, temporal_limits):\n",
    "    '''Download satellite data within a specified geographic and temporal range from NASA EarthData'''\n",
    "    \n",
    "    #Requires user to set up an account\n",
    "    earthaccess.login()  \n",
    "\n",
    "    #Includes moderate resolution, imagery resolution, and day-night-band data and geolocation files for NOAA-20 and Suomi-NPP\n",
    "    short_names = [\"VNP02MOD\", \"VNP02IMG\", \"VNP02DNB\", \"VNP03MOD\", \"VNP03IMG\", \"VNP03DNB\",\n",
    "                      \"VJ102MOD\", \"VJ102IMG\", \"VJ102DNB\", \"VJ103MOD\", \"VJ103IMG\", \"VJ103DNB\"]\n",
    "      \n",
    "    #Download data for each file type\n",
    "    for short_name in short_names:\n",
    "        if short_name[0:2] == 'VN':\n",
    "            version = 2\n",
    "        elif short_name[0:2] == 'VJ':\n",
    "            version = 2.1\n",
    "                \n",
    "        results = earthaccess.search_data(\n",
    "                    short_name=short_name, \n",
    "                    cloud_hosted=True, \n",
    "                    bounding_box=bounding_box,\n",
    "                    temporal=temporal_limits,\n",
    "                    version=version\n",
    "        )\n",
    "        \n",
    "        #Update to desired directory path\n",
    "        earthaccess.download(results, data_path)           \n",
    "            "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8312e64-8c83-4a3b-8eca-66e1c633888b",
   "metadata": {},
   "source": [
    "## Define the projection\n",
    "All images will be resampled to a common projection. This is required because the tracking algorithms will operate on a pixel-by-pixel basis. This function uses the [create_area_def](https://pyresample.readthedocs.io/en/latest/howtos/geometry_utils.html#areadefinition-creation) function from [pyresample package](https://pyresample.readthedocs.io/en/latest/index.html) (which underlies resampling operations in SatPy, but can be used as a standalone package as well).\n",
    "\n",
    "Here, data are reprojected to the NSIDC Sea Ice Polar Stereographic North [(EPSG:3413)](https://nsidc.org/data/user-resources/help-center/guide-nsidcs-polar-stereographic-projection) projection. A horizonal resolution of 750 m is chosen. \n",
    "\n",
    "I ran into a little trouble implementing this projection which a SatPy engineer helped sort out. Briefly, the \"corners\" of the projection have to be input from the bounding box. This is likely because SatPy's create_area_def function is getting confused by the polar-centric projection. See this [github discussion](https://github.com/pytroll/tutorial-satpy-half-day/issues/13) for more information about this fix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "78fcb850-a5d6-4811-9f72-e02a8616a5b5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "def setup_projection(bounding_box):\n",
    "    '''Set up a custom area and projection for the images. Input bounding box is (minlon, minlat, maxlon, maxlat)'''\n",
    "\n",
    "    #Directly defining the extents from the bounding_box doesn't work with this projection. \n",
    "    #Instead, use transformer.transform to define corners \n",
    "    #See https://github.com/pytroll/tutorial-satpy-half-day/issues/13\n",
    "    \n",
    "    #Instead, need to manually determine corners\n",
    "    transformer = Transformer.from_crs(4326, 3413, always_xy=True) #3413 is the EPSG code of the desired projection\n",
    "\n",
    "    #Define 'corners'\n",
    "    (x1, y1) = transformer.transform(bounding_box[0], bounding_box[1])\n",
    "    (x2, y2) = transformer.transform(bounding_box[0], bounding_box[3])\n",
    "    (x3, y3) = transformer.transform(bounding_box[2], bounding_box[1])\n",
    "    (x4, y4) = transformer.transform(bounding_box[2], bounding_box[3])\n",
    "\n",
    "    x = (x1, x2, x3, x4)\n",
    "    y = (y1, y2, y3, y4)\n",
    "    extents = (min(x), min(y), max(x), max(y))\n",
    "           \n",
    "    area_id = 'nsidc_polar_north'\n",
    "    description = 'NSIDC Sea Ice Polar Stereographic North'\n",
    "    projection = 'EPSG:3413' \n",
    "    resolution = 750 #in meters\n",
    "\n",
    "    my_area = create_area_def(area_id=area_id, description=description, projection=projection, \n",
    "                    resolution=resolution, area_extent=extents)    \n",
    "         \n",
    "    print(my_area)\n",
    "    print(\"Pixel size: x = \", my_area.pixel_size_x, \"m, y = \", my_area.pixel_size_y, \"m\") #horizontal resolution in meters\n",
    "\n",
    "    return my_area"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2bf9ac6-9d99-4fc2-9335-e7852ffe8b4f",
   "metadata": {},
   "source": [
    "### Extract the date and time of each set of images from file names\n",
    "String manipulation function that determines the date and time of each file. The time format is YYYYDDD.HHMM. Returns the list of unique file times (there may multiple files at each file time)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "366e4ebb-9db3-4bbf-ab04-18ac2dfa178d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_file_times(filenames): \n",
    "    '''Generate a list of unique, sorted day and time strings for all files. Result can include times for both VNP and VJ files'''\n",
    "    \n",
    "    VNP_pattern = r'\\.A(\\d{7}).*(\\d{4})\\.002\\.'  #.002 is version\n",
    "    VJ_pattern = r'\\.A(\\d{7}).*(\\d{4})\\.021\\.'   #.021 is version \n",
    "    \n",
    "    time_strings = [None] * len(filenames)  # Preallocate with None\n",
    "    for i, filename in enumerate(filenames):\n",
    "                \n",
    "        if 'VN' in filename:\n",
    "            match = re.search(VNP_pattern, filename)\n",
    "        elif 'VJ' in filename:\n",
    "            match = re.search(VJ_pattern, filename)\n",
    "        \n",
    "        if match: \n",
    "            day = match.group(1)\n",
    "            hour = match.group(2)\n",
    "            time_strings[i] = day + '.' + hour\n",
    "        else:\n",
    "            print('Cannot determine time for files: ', filenames)\n",
    "    \n",
    "    unique_sorted_time_strings = sorted(set(time_strings))\n",
    "    return(unique_sorted_time_strings) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3919245-f611-4b00-a69b-f32fee3af5f7",
   "metadata": {},
   "source": [
    "## Return the names of all files at a given timestep \n",
    "Processed images will be made sequentially by iterating through all *times* for which there is imagery. Multiple files (i.e., imagery-resolution, moderate-resolution, and DNB files) are generated at each image time and need to be grouped together into a single Satpy Scene. Sometimes we also want to group together images that are taken in sequence along the satellite swath. These files are adjacent images taken six minutes apart. This function will check when the previous/next images were taken and group files together if they are on the same swath. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "3daeee2c-d741-4fbf-afa8-5113ff0b2b04",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_filenames_at_current_time(data_path, i, filetimes):\n",
    "    \n",
    "    '''Generate a list of all filenames at the current timestep. If two sets of files\n",
    "        are taken six minutes apart, assume they are part of the same swath and group them together, using the earlier time.'''\n",
    "    \n",
    "    filetime = filetimes[i]\n",
    "    filenames = glob(data_path+'*'+filetime+'*.nc') #Files at this timestep\n",
    "    filetime_dt = datetime.strptime(filetime, \"%Y%j.%H%M\")\n",
    "\n",
    "    #Check if the next batch of files is within six minutes and should be combined\n",
    "    if i < len(filetimes)-1:\n",
    "        nextfiletime = filetimes[i+1]\n",
    "        \n",
    "        # Convert timestamp strings to datetime objects\n",
    "        nextfiletime_dt = datetime.strptime(nextfiletime, \"%Y%j.%H%M\")\n",
    "\n",
    "        # Calculate time difference\n",
    "        time_difference = nextfiletime_dt - filetime_dt\n",
    "        combine_files = time_difference <= timedelta(minutes=6) \n",
    "        \n",
    "        if combine_files: \n",
    "            print(\"Combining files at \" + filetime[-4:] + ' and ' + nextfiletime[-4:])\n",
    "            nextfilenames = glob(data_path+'L1b/*'+nextfiletime+'*nc')  \n",
    "            filenames.extend(nextfilenames)\n",
    "    \n",
    "    #Check if this batch of files was already combined with the previous batch of files\n",
    "    if i > 0:\n",
    "        prevfiletime = filetimes[i-1]\n",
    "        prevfiletime_dt = datetime.strptime(prevfiletime, \"%Y%j.%H%M\")\n",
    "\n",
    "        # Calculate time difference\n",
    "        time_difference = filetime_dt - prevfiletime_dt\n",
    "        already_combined = time_difference <= timedelta(minutes=6)\n",
    "        if already_combined:\n",
    "            print(\"Files at \" + filetime[-4:] + ' were aleady combined with ' + prevfiletime[-4:]+ '. Will not reprocess the files')\n",
    "            filenames = None\n",
    "            \n",
    "    #Use version 2 data for VNP - necessary to have this at the end of the function in case files from the next timestep were added\n",
    "    if filenames is not None:\n",
    "        filenames = [filename for filename in filenames if '.001.' not in filename]\n",
    "    \n",
    "    return filenames"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "932fb67e-923a-409d-9ad5-5ee18f62efcf",
   "metadata": {},
   "source": [
    "### Determine if the image within the region of interest is from day or night\n",
    "\n",
    "Determining if the image is daylit is necessary to decide if we want to use Adaptive DNB (for low light conditions) or true and false color imagery (for daylit conditions). \n",
    "\n",
    "While the image has a DayNightFlag attribute, this is for the entire large image and can be \"both\". We want to know if it is day or night within our smaller subregion of interest. From some experimentation, it seems using a maximum solar zenith angle is an effective way to discriminate, with large angles indicating nighttime values (somewhat unintuitive!)\n",
    "\n",
    "*Note: At the moment I resample the scene here, which seems inefficient as it is later resampled again. However, I think this is unavoidable as the scene must be resampled AFTER the composite is selected, and day or night must first be determined to select which composite to make.* \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "4f6dc625-7995-4bff-b6e6-653ec42eb96e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def contains_night_pixels(scn, my_area):\n",
    "    '''Determine if the resampled scene region contains any nighttime pixels''' \n",
    "    \n",
    "    solar_zenith_angles = ['solar_zenith_angle', 'dnb_solar_zenith_angle']\n",
    "    available = scn.available_dataset_names()\n",
    "    available_zenith_angles = list(filter(lambda item: item in available, solar_zenith_angles))\n",
    "    \n",
    "    # Define a threshold solar zenith angle to distinguish day and night\n",
    "    max_zenith_angle_threshold  = 88.5 #This is from experimenting\n",
    "    \n",
    "    scn.load([available_zenith_angles[0]], generate=False)\n",
    "    resampled_scene = scn.resample(my_area)\n",
    "    \n",
    "    #From experimenting, it seems like the max angle is important in this case. \n",
    "    max_zenith_angle = np.nanmax(resampled_scene[available_zenith_angles[0]].values)\n",
    "       \n",
    "    #print('Maximum solar angle in resampled scene: ', max_zenith_angle)\n",
    "\n",
    "    if max_zenith_angle >= max_zenith_angle_threshold or isnan(max_zenith_angle): \n",
    "        print('Resampled scene contains some night time pixels')\n",
    "        return True\n",
    "    else: \n",
    "        print('Resampled scene contains only daytime pixels')\n",
    "        return False\n",
    "   "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "214c0cca-0084-4bc5-bd4d-8dbc2ec53682",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Download data files, create and save images\n",
    "Download data from the VIIRS instruments on the S-NPP and NOAA-20 satellites. Make desired composites (i.e., true color images) and save as GeoTIFF files. \n",
    "\n",
    "#### Defining a time range and region of interest\n",
    "The geographic region included in the image can be adjusted by the bounding_box parameter. The time range to make images for can be adjusted by the temporal_limits parameter \n",
    "\n",
    "#### Projection\n",
    "The default projection is the NSIDC Sea Ice Polar Stereographic North projection with horizonal resolution of 750 m. This can potentially be changed within the setup_projection function (defined above). \n",
    "\n",
    "#### Image composites\n",
    "The code determines if the region in bounding_box was daylit/nighttime when the image was taken. If daytime, true color and false color images are made. If nighttime, adaptive DNB images are made. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "87c63632-f4fb-4af8-ba51-1e79fb4d238b",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:root:shape found from radius and resolution does not contain only integers: (1225.4444314487023, 1161.3192266658189)\n",
      "Rounding shape to (1226, 1162) and resolution from (750.0, 750.0) meters to (749.5606024090914, 749.6601334310984) meters\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Area ID: nsidc_polar_north\n",
      "Description: NSIDC Sea Ice Polar Stereographic North\n",
      "Projection: {'datum': 'WGS84', 'lat_0': '90', 'lat_ts': '70', 'lon_0': '-45', 'no_defs': 'None', 'proj': 'stere', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}\n",
      "Number of columns: 1162\n",
      "Number of rows: 1226\n",
      "Area extent: (-1891420.7102, -265821.8454, -1020431.2902, 653261.4782)\n",
      "Pixel size: x =  749.5606024090915 m, y =  749.6601334310984 m\n",
      "./data_test/\n",
      "Beginning processing for images at 2019187.0030\n",
      "Resampled scene contains only daytime pixels\n",
      "true_color image at 2019187.0030 already exists. Will not remake image.\n",
      "false_color image at 2019187.0030 already exists. Will not remake image.\n",
      "Beginning processing for images at 2019187.0118\n",
      "Combining files at 0118 and 0124\n",
      "Resampled scene contains only daytime pixels\n",
      "Saving true_color at 2019187.0118\n",
      "Saving false_color at 2019187.0118\n",
      "Beginning processing for images at 2019187.0124\n",
      "Files at 0124 were aleady combined with 0118. Will not reprocess the files\n",
      "Beginning processing for images at 2019187.0812\n",
      "Resampled scene contains only daytime pixels\n",
      "Saving true_color at 2019187.0812\n",
      "Saving false_color at 2019187.0812\n",
      "Beginning processing for images at 2019187.0900\n",
      "Combining files at 0900 and 0906\n",
      "Resampled scene contains only daytime pixels\n",
      "Saving true_color at 2019187.0900\n",
      "Saving false_color at 2019187.0900\n",
      "Beginning processing for images at 2019187.0906\n",
      "Files at 0906 were aleady combined with 0900. Will not reprocess the files\n",
      "Beginning processing for images at 2019187.0954\n",
      "Resampled scene contains only daytime pixels\n",
      "Saving true_color at 2019187.0954\n",
      "Saving false_color at 2019187.0954\n",
      "Beginning processing for images at 2019187.1042\n",
      "Combining files at 1042 and 1048\n",
      "Resampled scene contains only daytime pixels\n",
      "Saving true_color at 2019187.1042\n",
      "Saving false_color at 2019187.1042\n",
      "Beginning processing for images at 2019187.1048\n",
      "Files at 1048 were aleady combined with 1042. Will not reprocess the files\n",
      "Beginning processing for images at 2019187.1136\n",
      "Resampled scene contains only daytime pixels\n",
      "Saving true_color at 2019187.1136\n",
      "Saving false_color at 2019187.1136\n"
     ]
    }
   ],
   "source": [
    "#Name the run and directories for data, images \n",
    "run = 'test'\n",
    "data_path = './data_' + run + '/'\n",
    "image_save_path = './images_' + run + '/'\n",
    "\n",
    "#Set up time and geographic range\n",
    "bounding_box=(-155, 72.5, -127, 80) #Central Beaufort Sea\n",
    "temporal_limits= ('2019-07-06 00:00', '2019-07-06 12:00') #Example time\n",
    "\n",
    "#Download files (set destination path within download_files)\n",
    "#download_files(data_path=data_path, bounding_box=bounding_box, temporal_limits=temporal_limits)\n",
    "\n",
    "#Setup projection for the region of interest for the images\n",
    "my_area = setup_projection(bounding_box)  \n",
    "\n",
    "#Create a list of times at which images were taken. These are directly from the file names (the format is YYYYDDD.HHMM)\n",
    "all_VNP_filenames = glob(data_path+'VNP*.002.*nc') #S-NPP data. 002 indicates Version 2\n",
    "all_VJ_filenames = glob(data_path+'VJ*.021.*nc') #NOAA-20 data. 021 indicates Version 2.1\n",
    "filetimes = get_file_times(all_VNP_filenames + all_VJ_filenames) #If you just want images from one instrument, you can use e.g., filetimes = get_file_times(all_VJ_filenames)\n",
    "\n",
    "print(data_path)\n",
    "\n",
    "#Make images for each satellite pass\n",
    "for i, filetime in enumerate(filetimes):\n",
    "     \n",
    "    print() # This will print a blank line to distinguish this time's information\n",
    "    print('Beginning processing images taken at ' + filetime)\n",
    "    \n",
    "    #All files at the current time - combine files on same swath but not the same image (i.e., combine pairs of images taken six minutes apart)\n",
    "    filenames = get_filenames_at_current_time(data_path, i, filetimes)\n",
    "    \n",
    "    #filenames will be None if the timestep was part of the same swath as the previous image and was already plotted at previous timestep\n",
    "    if filenames is None: \n",
    "        continue  \n",
    "    \n",
    "    # Try to build a Satpy Scene using all files at this timestep. \n",
    "    # Skip this timestep if this fails (I have gotten \"possibly corrupted file\" errors, I think due to downloading errors)\n",
    "    try:\n",
    "        scn = Scene(filenames=filenames, reader='viirs_l1b')           \n",
    "    except Exception as e:\n",
    "        print(f\"An error occurred while processing files at '{filetime}': {e}\")\n",
    "        continue\n",
    "    \n",
    "    # Decide what composites to plot\n",
    "    if contains_night_pixels(scn, my_area):\n",
    "        channels = ['adaptive_dnb']\n",
    "    else:\n",
    "        channels = ['true_color', 'false_color'] \n",
    "        \n",
    "    #Remove any datasets composites that are not available for this scene\n",
    "    available = scn.available_composite_names() + scn.available_dataset_names()\n",
    "    channels = list(filter(lambda item: item in available, channels))\n",
    "      \n",
    "    #Iterate through the desired composites and make images\n",
    "    for channel in channels:\n",
    "         \n",
    "        #Generate composite and resample to the common grid my_area defined in setup_projection\n",
    "        try:\n",
    " \n",
    "            #Data will have to be resampled before the composite can actually be generated. \n",
    "            #This is because different required bands have different resolution. \n",
    "            #generate=False allows us to \"load\" the ccomposite and then resample\n",
    "            # See also https://satpy.readthedocs.io/en/latest/faq.html#how-can-i-speed-up-creation-of-composites-that-need-resampling\n",
    "            scn.load([channel], generate=False)          \n",
    "            \n",
    "            #Resample to the common grid defined in my_area\n",
    "            resampled_scene = scn.resample(my_area)          \n",
    "            \n",
    "            #Make and save the image          \n",
    "            savename = f\"{image_save_path}{filetime}_{channel}.tiff\"    \n",
    "            if os.path.exists(savename): #Dont remake this figure\n",
    "                print(channel+' image at ' + filetime + ' already exists. Will not remake image.')\n",
    "            \n",
    "            else:\n",
    "                print('Saving ' + channel + ' at ' + filetime)# + 'as' + savename)\n",
    "                \n",
    "                #Save the image. fill_value=0 excludes the alpha channel and makes an RGB image (techinically it means \"set invalid values to black\")\n",
    "                resampled_scene.save_dataset(channel,  writer='geotiff', filename=savename, fill_value=0)\n",
    "\n",
    "                    \n",
    "        except Exception as e:\n",
    "            print(f\"An error occurred while processing '{channel}' for '{filetime}': {e}\")\n",
    "            "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Sea Ice Oscillations",
   "language": "python",
   "name": "sea_ice_oscillations"
  },
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}