abrazeau/gdal_warp_Tiffs_Batch.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Projecting the National Elevation Dataset to North America Albers Equal Area Conic using GDAL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Required:\n",
    " - Unprojected NED Cloud Optimized GeoTIFFS (COGs)\n",
    "     - A globebd list of COGs that covers the whole country (1/3rd arc, 1/9th arc, 1m)\n",
    " - GDAL Warp command line function\n",
    " \n",
    "## Workflow\n",
    "\n",
    " 1. Set paths for input rasters, output porjected rasters and the .cmd file that will do all the hard work for us\n",
    " 2.  Glob for a list of all unprojected COGs\n",
    " 3.  Define the function that will write the .cmd file\n",
    "     - Easiest to write it to the same folder that the gdal tool is located in\n",
    " 4.  Write the .cmd file\n",
    " 5. Open cmd prompt, cd to the location of the folder containing both the gdal tool and the .cmd file, run the .cmd file\n",
    "     - Better than just opening the .cmd so that you can see the cmd prompt log and identify errors when it is finished/crashes\n",
    " 6. Check to see if any COGs were missed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the necessary modules\n",
    "import os\n",
    "from glob import glob"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Workflow step 1.\n",
    "#### Set paths for input rasters, output porjected rasters and the .cmd file that will do all the hard work for us"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set paths\n",
    "inputs_tif_dir = r'T:\\CCSI\\TECH\\FEMA\\2017_NFIP_RiskRating_Milliman\\GIS\\Harddrive\\NED\\NED_COGTIFFS'\n",
    "outputs_tif_dir = r'T:\\CCSI\\TECH\\FEMA\\2017_NFIP_RiskRating_Milliman\\GIS\\Harddrive\\NED\\NED_COGTIFFS_AEA'\n",
    "cmd_file = r'C:\\Users\\abrazeau\\Desktop\\test_warp_processing\\RunBatch.cmd'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Workflow step 2.\n",
    "#### Glob for a list of all unprojected COGs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13199 Cloud Optimized GeoTIFFS found in input directory\n"
     ]
    }
   ],
   "source": [
    "# File types, globbed list of COGs\n",
    "file_type = '.tif'\n",
    "tifs = glob(os.path.join(inputs_tif_dir, '*{}'.format(file_type)))\n",
    "\n",
    "# SANITY CHECK\n",
    "print('{} Cloud Optimized GeoTIFFS found in input directory'.format(len(tifs)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Workflow step 3.\n",
    "#### Define the function that will write the .cmd file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MAIN FUNCTION\n",
    "def WriteCMD(cmd_file, raster_paths, file_type, gdal_function, new_file_dir):\n",
    "    '''This function opens a .cmd file in write mode, goes through a list of rasters and writes a command to use the gdal warp\n",
    "    tool on each raster. This tool projects a raster to a specified coordinate system, then the function places the projected COG\n",
    "    in 'new_file_dir'. The new GeoTIFFs are saved using their same file name from before projection.'''\n",
    "    quotes = '\"'\n",
    "    with open(cmd_file, 'w') as f:\n",
    "        for tif in raster_paths:\n",
    "            cmd = '{} '.format(gdal_function)\n",
    "            crs_options = ' -t_srs \"+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs' #NA Albers Eq. Area PROJ4 crs\n",
    "            gdal_options = ' -co TILED=YES -co COPY_SRC_OVERVIEWS=YES -co COMPRESS=LZW' # Maintain Cloud Optimized GeoTIFF properties\n",
    "            if '.tif' in tif:\n",
    "                raw_file_name = os.path.basename(tif).strip(file_type)\n",
    "            else:\n",
    "                raw_file_name = tif.split('\\\\')[-2]\n",
    "            new_file_name = os.path.join(new_file_dir, raw_file_name + '.tif')\n",
    "            terminal_input =cmd + crs_options + quotes+ ' ' + quotes + tif +quotes+ ' '+ quotes + new_file_name + quotes + gdal_options + '\\n'\n",
    "            f.write(terminal_input)\n",
    "    return None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Workflow step 4.\n",
    "#### Write the .cmd file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The .cmd file to project all NED rasters can be found at C:\\Users\\abrazeau\\Desktop\\test_warp_processing\\RunBatch.cmd\n"
     ]
    }
   ],
   "source": [
    "# GDAL COMMAND\n",
    "gdalwarp = 'gdalwarp'\n",
    "\n",
    "WriteCMD(cmd_file, tifs, file_type, gdalwarp, outputs_tif_dir)\n",
    "\n",
    "print ('The .cmd file to project all NED rasters can be found at {}'.format(cmd_file))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Workflow step 5.\n",
    "#### Open cmd prompt, cd to the location of the folder containing both the gdal tool and the .cmd file, run the .cmd file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Workflow step 6.\n",
    "#### Check to see if any COGs were missed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are  0  COGs missing...\n"
     ]
    }
   ],
   "source": [
    "itifs = glob(os.path.join(inputs_tif_dir, '*{}'.format(file_type)))\n",
    "otifs = glob(os.path.join(outputs_tif_dir, '*{}'.format(file_type)))\n",
    "\n",
    "print ('There are ', len(itifs) - len(otifs), ' COGs missing...')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run this cell to find which COGs are missing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If the cell above prints any number but 0, run this cell to find which files are missing\n",
    "tif_files = []\n",
    "for itif in itifs:\n",
    "    tif_files.append(tif.split('NED_COGTIFFS\\\\')[1])\n",
    "\n",
    "otif_files = []\n",
    "for otif in otifs:\n",
    "    otif_files.append(otif.split('AEA\\\\')[1])\n",
    "    \n",
    "otifs_items = set(otif_files)\n",
    "missing = [item for item in tif_files if item not in otifs_items]\n",
    "print( 'The following COGs were missed: ', missing)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "##### CRS projection = North America Albers Equal Area Conic, ESRI:102008, info at https://epsg.io/102008\n",
    "##### Adding 'Start ' to the beginning of each line in the .cmd file begins all wraps at once.. would crash the machine\n",
    "##### Cloud Optimized GeoTIFF information at: https://trac.osgeo.org/gdal/wiki/CloudOptimizedGeoTIFF#HowtoreaditwithGDAL\n",
    "##### gdalwarp info at: http://www.gdal.org/gdalwarp.html"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:geos]",
   "language": "python",
   "name": "conda-env-geos-py"
  },
  "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}