Take HLS data and write to GeoTiff.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "May 6, 2020\n",
    "\n",
    "This notebook will loop through a given system folder and look for hdf files\n",
    "from the HLS dataset and convert them to geotiffs with masking applied.\n",
    "\n",
    "Notes:\n",
    "- Water is not masked, but the code to do so is merely commented out (easily reinstated) \n",
    "- For L30 data, the hdr file is REQUIRED\n",
    "- There is a threshold option for excluding images with too high percentage masked\n",
    "- A list of images excluded from converstion is generated "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, re, sys, datetime\n",
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "from pyhdf.SD import SD, SDC\n",
    "import gdal\n",
    "import rasterio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Location of HLS data to be converted \n",
    "system_folder = 'filepath to HLS data folder'\n",
    "\n",
    "#Percent of image masked for which an image in the folder \n",
    "# is skipped and not convered to geotiff\n",
    "masking_exclusion_threshold = .5  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_hdf_data(file_obj, band_name, scaled=True):\n",
    "    sds_obj = file_obj.select(band_name)\n",
    "    data = sds_obj.get()\n",
    "    attributes = sds_obj.attributes()\n",
    "    if scaled:\n",
    "        add_offset = float(attributes['add_offset'])\n",
    "        scale_factor = float(attributes['scale_factor'])\n",
    "        data = (data-add_offset)*scale_factor\n",
    "    return data\n",
    "\n",
    "def QA_mask(file, product_id):\n",
    "    if product_id == 'L30':\n",
    "        blue = read_hdf_data(file, 'band02')\n",
    "        green = read_hdf_data(file, 'band03')\n",
    "        red = read_hdf_data(file, 'band04')\n",
    "        nir = read_hdf_data(file, 'band05')\n",
    "    else:\n",
    "        blue = read_hdf_data(file, 'B02')\n",
    "        green = read_hdf_data(file, 'B03')\n",
    "        red = read_hdf_data(file, 'B04')\n",
    "        nir = read_hdf_data(file, 'B8A')\n",
    "    \n",
    "    # Remove pixels covered by cirrus, clouds, cloud shadows, snow/ice, water and thick aerosol.\n",
    "    qa = read_hdf_data(file, 'QA', scaled=False)\n",
    "    cirrus = np.right_shift(np.bitwise_and(qa, int('00000001', 2)), 0)\n",
    "    cloud = np.right_shift(np.bitwise_and(qa, int('00000010', 2)), 1)\n",
    "    cloud_shadow = np.right_shift(np.bitwise_and(qa, int('00001000', 2)), 3)\n",
    "    snow_ice = np.right_shift(np.bitwise_and(qa, int('00010000', 2)), 4)\n",
    "    # water = np.right_shift(np.bitwise_and(qa, int('00100000', 2)), 5)  Commented this out because I use water\n",
    "    aerosol = np.right_shift(np.bitwise_and(qa, int('11000000', 2)), 6)\n",
    "    mask1 = (cirrus==1)|(cloud==1)|(cloud_shadow==1)|(snow_ice==1)|(aerosol==3)# |(water==1)\n",
    "\n",
    "    # Remove pixels with a reflectance out of the range 0.0-1.0.\n",
    "    mask2 = (blue<=0.0)|(blue>=1.0)|(green<=0.0)|(green>=1.0)|(red<=0.0)|(red>=1.0)|(nir<=0.0)|(nir>=1.0)\n",
    "\n",
    "    # Make the final mask.\n",
    "    mask = mask1|mask2\n",
    "    mask = ~mask1\n",
    "\n",
    "    del mask1, mask2, qa, cirrus, cloud, cloud_shadow, snow_ice, aerosol, # water    \n",
    "   \n",
    "    return mask\n",
    "\n",
    "def read_hdf_to_arr(hdf_path, band, datatype=np.int16):\n",
    "    \"\"\"\n",
    "    Functionalizing the process of reading HDF files into arrays\n",
    "    read a single band out of the hdf and load it into a numpy array\n",
    "    \"\"\"\n",
    "    if os.path.isfile(hdf_path):\n",
    "        src = gdal.Open(hdf_path)\n",
    "        band_ds = gdal.Open(src.GetSubDatasets()[band][0], gdal.GA_ReadOnly)\n",
    "        band_array = band_ds.ReadAsArray().astype(datatype)\n",
    "        print(band_array)\n",
    "        del src\n",
    "        return band_array\n",
    "    else:\n",
    "        print(\"That file does not exist\")\n",
    "        return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading the metadata of a hdf file.\n",
      "This is an S30 scene.\n",
      "    File name:HLS.S30.T12SUC.2020148.v1.4.hdf\n",
      "    Cloud coverage: 1\n",
      "    Upper-left corner map coordinates (X, Y in m): 300000.0, 3800040.0\n",
      "    Pixel size (in m): 30.0\n",
      "    Image rows and columns: 3660, 3660\n",
      "    UTM zone: 12\n",
      "    EPSG code is: 32612\n",
      "    Percent of image masked: 0.01225058974588671\n",
      "[[ 706  563  593 ...  655  656  651]\n",
      " [ 762  611  592 ...  645  655  645]\n",
      " [ 673  607  601 ...  658  658  659]\n",
      " ...\n",
      " [1091 1116 1131 ... 1972 1740 1548]\n",
      " [1128 1139 1127 ... 1826 1724 1715]\n",
      " [1039  999  984 ... 1438 1832 1833]]\n",
      "[[ 971  797  841 ...  899  882  895]\n",
      " [1041  846  822 ...  877  888  879]\n",
      " [ 942  858  841 ...  901  892  909]\n",
      " ...\n",
      " [1510 1532 1558 ... 2395 2174 2006]\n",
      " [1570 1576 1545 ... 2283 2107 2113]\n",
      " [1437 1382 1373 ... 1746 2252 2326]]\n",
      "[[1258 1016 1111 ... 1246 1246 1259]\n",
      " [1342 1093 1058 ... 1218 1243 1242]\n",
      " [1245 1126 1093 ... 1224 1239 1268]\n",
      " ...\n",
      " [2089 2126 2141 ... 2859 2697 2435]\n",
      " [2163 2177 2152 ... 2728 2609 2599]\n",
      " [2014 1933 1933 ... 2110 2824 2885]]\n",
      "[[2069 1862 2000 ... 1905 1890 1980]\n",
      " [2185 1860 1801 ... 1856 1852 1903]\n",
      " [2083 1837 1787 ... 1844 1834 1909]\n",
      " ...\n",
      " [2976 2998 3006 ... 3234 3145 3239]\n",
      " [3031 3013 2938 ... 3258 3093 3043]\n",
      " [2858 2832 2808 ... 2752 3201 3281]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/sarahwegmueller/anaconda3/envs/Image_prep/lib/python3.7/site-packages/rasterio/__init__.py:219: NotGeoreferencedWarning: Dataset has no geotransform set. The identity matrix may be returned.\n",
      "  s = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)\n"
     ]
    }
   ],
   "source": [
    "#Prepare a list to populate with images that are 50% masked or greater\n",
    "high_masking = []\n",
    "\n",
    "for dirpath, dirnames, filenames in os.walk(system_folder):\n",
    "    if not filenames:\n",
    "        continue\n",
    "    for f in filenames:\n",
    "        if f.endswith(\".hdf\"):\n",
    "            pathname = os.path.join(dirpath, f)\n",
    "            outname = pathname[:-4]+\".tif\"\n",
    "            \n",
    "            print(\"Reading the metadata of a hdf file.\")\n",
    "    \n",
    "            product_id = os.path.basename(pathname).split(\".\")[1]\n",
    "            print(\"This is an \" + product_id + ' scene.')\n",
    "            print('    File name:' + f)\n",
    "\n",
    "            try:\n",
    "                file = SD(pathname, SDC.READ)\n",
    "            except:\n",
    "                raise IOError(\"Cannot open the hdf file: %s\" %pathname)\n",
    "\n",
    "            metadata = file.attributes()\n",
    "\n",
    "            # acquisition time\n",
    "            # if product_id == 'L30':\n",
    "            #    tm = datetime.datetime.strptime(metadata['SENSING_TIME'], '%Y-%m-%dT%H:%M:%S.%f0Z')\n",
    "            #    year = str(tm.timetuple().tm_year)\n",
    "            #    doy = tm.timetuple().tm_yday\n",
    "            #    print(\"    Acqusition time (YYYY-MM-DD hh:mm:ss): %s\" %tm)\n",
    "\n",
    "            # else:\n",
    "            #    tm = datetime.datetime.strptime(metadata['SENSING_TIME'], '%Y-%m-%dT%H:%M:%S.%fZ')\n",
    "            #    year = str(tm.timetuple().tm_year)\n",
    "            #    doy = tm.timetuple().tm_yday\n",
    "            #    print(\"    Acqusition time (YYYY-MM-DD hh:mm:ss): %s\" %tm)\n",
    "                \n",
    "            # get cloud cover percent\n",
    "            print('    Cloud coverage: ' + str(metadata['cloud_coverage']))\n",
    "\n",
    "            # upper-left corner map coordinates\n",
    "            ulx      = float(metadata['ULX'])\n",
    "            uly      = float(metadata['ULY'])\n",
    "            print(\"    Upper-left corner map coordinates (X, Y in m): %s, %s\" %(ulx, uly))\n",
    "\n",
    "            # pixel size\n",
    "            p_size   = float(metadata['SPATIAL_RESOLUTION'])\n",
    "            print(\"    Pixel size (in m): %s\" %p_size)\n",
    "\n",
    "            # image rows and columns\n",
    "            nrows    = int(metadata['NROWS'])\n",
    "            ncols    = int(metadata['NCOLS'])\n",
    "            print(\"    Image rows and columns: %s, %s\" %(nrows, ncols))\n",
    "\n",
    "            # utm zone\n",
    "            utm_zone = int(re.search('ZONE \\d+', metadata['HORIZONTAL_CS_NAME'].upper()).group()[-2:])\n",
    "            print(\"    UTM zone: %s\" %utm_zone)\n",
    "\n",
    "            # Get EPSG code from S30 image; L30 does not have this option\n",
    "            if product_id == 'S30':\n",
    "                epsg = metadata['HORIZONTAL_CS_CODE'].split(':')\n",
    "                epsg = (int(epsg[1]))\n",
    "                print('    EPSG code is: ' + str(epsg))\n",
    "            \n",
    "            if product_id == 'L30':\n",
    "                bands = [1,2,3,4]   \n",
    "            else:\n",
    "                bands = [1,2,3,8] # b g r n (index 8 is the 8A band using zero-indexing)\n",
    "                \n",
    "            # Get Mask\n",
    "            qa_mask = QA_mask(file, product_id)  # turn the array into a binary mask\n",
    "                     \n",
    "            # What percent of the image is masked?\n",
    "            mask_vals = np.unique(qa_mask, return_counts=True)\n",
    "            percent_masked = (mask_vals[1][0])/(sum(mask_vals[1]))\n",
    "            print('    Percent of image masked: ' + str(percent_masked))\n",
    "                                    \n",
    "            # If image is over a given threshold percent of masking, put it in a list\n",
    "            if percent_masked >= masking_exclusion_threshold:\n",
    "                high_masking.append(f)\n",
    "            \n",
    "            else:\n",
    "                qa_mask = np.invert(qa_mask) #Inverstion is necessary to work with ma.masked_array\n",
    "                negatives = np.zeros_like(qa_mask, dtype=bool, subok=False)  # get the dimensions of the image from the QA mask\n",
    "                arr = np.zeros((len(bands), np.shape(qa_mask)[0], np.shape(qa_mask)[1]))  # pre-allocate\n",
    "                \n",
    "                # apply QA mask and remove all negative values\n",
    "                for count, b in enumerate(bands):\n",
    "                    tmp_band = read_hdf_to_arr(pathname, b)\n",
    "                    tmp = ma.masked_array(tmp_band, qa_mask)\n",
    "                    arr[count,:,:] = ma.filled(tmp, 0)  # fill mask using nodata value. add to 'arr' to create a 4-band image\n",
    "                    # find all the negative pixels\n",
    "                    negatives = negatives + (tmp_band<0)  # True where there are negative values. False elsewhere.\n",
    "                arr = ma.masked_array(arr, mask=np.broadcast_to(negatives[np.newaxis,:,:], arr.shape))\n",
    "                arr = ma.filled(arr, 0)\n",
    "\n",
    "                # Write data to a geotiff\n",
    "                with rasterio.open(pathname) as src:\n",
    "                    kwds = src.profile\n",
    "                    kwds['nodata'] = 0\n",
    "                    kwds['driver'] = 'GTiff'\n",
    "                    kwds['dtype'] = rasterio.int16\n",
    "                    kwds['width'] = arr.shape[2]\n",
    "                    kwds['height'] = arr.shape[1]\n",
    "                    kwds['count'] = arr.shape[0]\n",
    "                    kwds['transform'] = rasterio.transform.from_origin(ulx, uly, 30., 30.) \n",
    "                    \n",
    "                    #Get crs info\n",
    "                    if product_id == 'S30':\n",
    "                        kwds['crs'] = rasterio.crs.CRS.from_epsg(epsg) \n",
    "                    else: \n",
    "                        #for some reason, L30 data is incomplete in the hdf file.  \n",
    "                        # Need to use hdr.\n",
    "                        hdr_fp = f + '.hdr'\n",
    "                        hdr = os.path.join(dirpath, hdr_fp)\n",
    "                        info = open(hdr, 'r')\n",
    "                        guts = info.read()\n",
    "                        pieces = guts.split('\\n')\n",
    "                        crs = pieces[11].split('=')\n",
    "                        kwds['crs'] = (crs[1][2:-1])\n",
    "\n",
    "                    with rasterio.open(outname, 'w', **kwds) as dst:\n",
    "                         dst.write(np.array(arr).astype(rasterio.int16))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}