Tile_April_20th-final.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.\n",
      "  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Fixing figure output for all figures\n",
    "rcParams['axes.labelsize'] = 12\n",
    "rcParams['xtick.labelsize'] = 12\n",
    "rcParams['ytick.labelsize'] = 12\n",
    "rcParams['legend.fontsize'] = 12\n",
    "rcParams['figure.dpi'] = 160\n",
    "rcParams['savefig.dpi'] = 200"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from osgeo import gdal\n",
    "import georasters as gr\n",
    "import hashlib\n",
    "import numpy as np\n",
    "import osr\n",
    "import os\n",
    "import glob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def resample_tif(tif_file, pixel_spacing=0.5):\n",
    "    source = gdal.Open(tif_file)\n",
    "    source.GetRasterBand(1).SetNoDataValue(-32767)\n",
    "    wgs84 = osr.SpatialReference() # slopppy \n",
    "    wgs84.ImportFromEPSG(3413)     # only for polar stereographic!\n",
    "    # Get the Geotransform vector\n",
    "    geo_t = source.GetGeoTransform ()\n",
    "    x_size = source.RasterXSize # Raster xsize\n",
    "    y_size = source.RasterYSize # Raster ysize\n",
    "    # Work out the boundaries of the new dataset in the target projection\n",
    "    ulx, uly = geo_t[0], geo_t[3]\n",
    "    lrx = geo_t[0] + geo_t[1]*x_size\n",
    "    lry = geo_t[3] + geo_t[5]*y_size\n",
    "    # Now, we create an in-memory raster\n",
    "    mem_drv = gdal.GetDriverByName('MEM')\n",
    "    # The size of the raster is given the new projection and pixel spacing\n",
    "    # Using the values we calculated above. Also, setting it to store one band\n",
    "    # and to use Float32 data type.\n",
    "    dest = mem_drv.Create('', int((lrx - ulx)/pixel_spacing), \\\n",
    "            int((uly - lry)/pixel_spacing), 1, gdal.GDT_Float32)\n",
    "    # Calculate the new geotransform\n",
    "    new_geo = ( ulx, pixel_spacing, geo_t[2], \\\n",
    "                uly, geo_t[4], -pixel_spacing )\n",
    "    # Set the geotransform\n",
    "    dest.SetGeoTransform(new_geo)\n",
    "    dest.SetProjection(source.GetProjection())\n",
    "    # Perform the projection/resampling \n",
    "    res = gdal.ReprojectImage(source, dest, None, None,\n",
    "                gdal.GRA_Bilinear)\n",
    "    graster = gr.GeoRaster(dest.GetRasterBand(1).ReadAsArray(),\n",
    "                           dest.GetGeoTransform(),0,\n",
    "                           wgs84,gdal.GDT_Float32)\n",
    "    return graster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def grid_write(graster):\n",
    "    rows, cols = graster.shape\n",
    "    for i in r_[0:rows:200]:\n",
    "        for j in r_[0:cols:200]:\n",
    "            temp = graster[i:i+200,j:j+200]\n",
    "            if np.sum(temp.raster.mask) > 0:\n",
    "                pass\n",
    "            else:\n",
    "                temp = temp - temp.mean()\n",
    "                temp.to_tiff(hashlib.sha256(temp.raster).hexdigest())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "input_files = \"/Users/grigsbye/Downloads/April_20th_2012_ortho_DEM/ftp.fireballit.com/Icebridge/20120420/terrain_val/*.tif\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/grigsbye/Downloads/April_20th_2012_ortho_DEM/ftp.fireballit.com/Icebridge/20120420/1000/output\n"
     ]
    }
   ],
   "source": [
    "cd \"/Users/grigsbye/Downloads/April_20th_2012_ortho_DEM/ftp.fireballit.com/Icebridge/20120420/output/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "u'/Users/grigsbye/Downloads/April_20th_2012_ortho_DEM/ftp.fireballit.com/Icebridge/20120420/output'"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for filename in glob.glob(input_files):\n",
    "    grid_write(resample_tif(filename))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}