Create RGB Natural Color Composite From Landsat 8.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create RGB Natural Color Composite from Landsat 8 Data\n",
    "Author: Cole Krehbiel  \n",
    "Last Updated: 10-01-017"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import glob\n",
    "import numpy as np\n",
    "from osgeo import gdal\n",
    "import scipy.misc as sm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set input directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "in_dir = '/Users/cole/Desktop/Website/Landsat/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Search directory for desired bands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "b2_file = glob.glob(in_dir + '**B2.TIF') # blue band\n",
    "b3_file = glob.glob(in_dir + '**B3.TIF') # green band\n",
    "b4_file = glob.glob(in_dir + '**B4.TIF') # red band"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define a function to normalize each band array by the min and max values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def norm(band):\n",
    "    band_min, band_max = band.min(), band.max()\n",
    "    return ((band - band_min)/(band_max - band_min))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loop through however many Landsat 8 obs are in the input directory "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for i in range(len(b2_file)):   \n",
    "    \n",
    "    # Open each band using gdal\n",
    "    b2_link = gdal.Open(b2_file[i])\n",
    "    b3_link = gdal.Open(b3_file[i])\n",
    "    b4_link = gdal.Open(b4_file[i])\n",
    "    \n",
    "    # call the norm function on each band as array converted to float\n",
    "    b2 = norm(b2_link.ReadAsArray().astype(np.float))\n",
    "    b3 = norm(b3_link.ReadAsArray().astype(np.float))\n",
    "    b4 = norm(b4_link.ReadAsArray().astype(np.float))\n",
    "    \n",
    "    # Create RGB\n",
    "    rgb = np.dstack((b4,b3,b2))\n",
    "    del b2, b3, b4\n",
    "    \n",
    "    # Visualize RGB\n",
    "    #import matplotlib.pyplot as plt\n",
    "    #plt.imshow(rgb)\n",
    "    \n",
    "    # Export RGB as TIFF file\n",
    "    # Important: Here is where you can set the custom stretch\n",
    "    # I use min as 2nd percentile and max as 98th percentile\n",
    "    sm.toimage(rgb,cmin=np.percentile(rgb,2),\n",
    "               cmax=np.percentile(rgb,98)).save(b2_file[i].split('_01_')[0]+'_RGB.tif')"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}