notebooks/unit06/unit06_05.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook was prepared by Cayetano Benavent, 2016."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Map Algebra with GDAL and NumPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from osgeo import gdal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "src_file01 = '../../data/temp_gfs/tmp_2m_k_20150903.tif'\n",
    "src_file02 = '../../data/temp_gfs/tmp_80m_k_20150903.tif'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read datasets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "src_ds01 = gdal.Open(src_file01)\n",
    "src_ds02 = gdal.Open(src_file02)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "src01_band = src_ds01.GetRasterBand(1)\n",
    "src02_band = src_ds02.GetRasterBand(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read bands as NumPy arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "src01_ar = src01_band.ReadAsArray()\n",
    "src02_ar = src02_band.ReadAsArray()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "src_ds01 = None\n",
    "src_ds02 = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(src01_ar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(src02_ar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 270.55001831,  270.55001831,  270.55001831, ...,  270.55001831,\n",
       "         270.55001831,  270.55001831],\n",
       "       [ 270.55001831,  270.55001831,  270.55001831, ...,  270.55001831,\n",
       "         270.55001831,  270.55001831],\n",
       "       [ 270.55001831,  270.55001831,  270.55001831, ...,  270.55001831,\n",
       "         270.55001831,  270.55001831],\n",
       "       ..., \n",
       "       [ 218.55000305,  218.55000305,  218.55000305, ...,  218.55000305,\n",
       "         218.55000305,  218.55000305],\n",
       "       [ 219.3500061 ,  219.3500061 ,  219.3500061 , ...,  219.3500061 ,\n",
       "         219.3500061 ,  219.3500061 ],\n",
       "       [ 219.75001526,  219.75001526,  219.75001526, ...,  219.75001526,\n",
       "         219.75001526,  219.75001526]], dtype=float32)"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "src01_ar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 269.8500061 ,  269.8500061 ,  269.8500061 , ...,  269.8500061 ,\n",
       "         269.8500061 ,  269.8500061 ],\n",
       "       [ 269.95001221,  269.95001221,  269.95001221, ...,  269.95001221,\n",
       "         269.95001221,  269.95001221],\n",
       "       [ 269.95001221,  269.95001221,  270.05001831, ...,  269.95001221,\n",
       "         269.95001221,  269.95001221],\n",
       "       ..., \n",
       "       [ 218.25001526,  218.25001526,  218.25001526, ...,  218.25001526,\n",
       "         218.25001526,  218.25001526],\n",
       "       [ 218.95001221,  218.95001221,  218.95001221, ...,  218.95001221,\n",
       "         218.95001221,  218.95001221],\n",
       "       [ 219.45001221,  219.45001221,  219.45001221, ...,  219.45001221,\n",
       "         219.45001221,  219.45001221]], dtype=float32)"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "src02_ar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Operating with bands as NumPy arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "res = (src01_ar -273.15) - (src02_ar -273.15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.70001221,  0.70001221,  0.70001221, ...,  0.70001221,\n",
       "         0.70001221,  0.70001221],\n",
       "       [ 0.6000061 ,  0.6000061 ,  0.6000061 , ...,  0.6000061 ,\n",
       "         0.6000061 ,  0.6000061 ],\n",
       "       [ 0.6000061 ,  0.6000061 ,  0.5       , ...,  0.6000061 ,\n",
       "         0.6000061 ,  0.6000061 ],\n",
       "       ..., \n",
       "       [ 0.29998779,  0.29998779,  0.29998779, ...,  0.29998779,\n",
       "         0.29998779,  0.29998779],\n",
       "       [ 0.3999939 ,  0.3999939 ,  0.3999939 , ...,  0.3999939 ,\n",
       "         0.3999939 ,  0.3999939 ],\n",
       "       [ 0.30000305,  0.30000305,  0.30000305, ...,  0.30000305,\n",
       "         0.30000305,  0.30000305]], dtype=float32)"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Store NumPy array result on a new raster dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dst_file = '/tmp/out_mapalgebra_gdal.tif'\n",
    "\n",
    "drv_gtif = gdal.GetDriverByName(\"GTiff\")\n",
    "src_ds = gdal.Open(src_file01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-180.0, 0.25, 0.0, 90.0, 0.0, -0.2498266296809986)\n",
      "1440 721\n"
     ]
    }
   ],
   "source": [
    "src_geot = src_ds.GetGeoTransform()\n",
    "\n",
    "x_sz = src_ds.RasterXSize\n",
    "y_sz = src_ds.RasterYSize\n",
    "\n",
    "print(src_geot)\n",
    "print(x_sz, y_sz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "datatype = gdal.GDT_Float32\n",
    "n_band = 1\n",
    "nodata_val = -99999\n",
    "\n",
    "dst_raster = drv_gtif.Create(dst_file, x_sz, y_sz, n_band, datatype)\n",
    "dst_raster.SetGeoTransform(src_geot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write NumPy results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "out_band = dst_raster.GetRasterBand(n_band)\n",
    "out_band.SetNoDataValue(nodata_val)\n",
    "out_band.WriteArray(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Close datasets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "src_ds = None\n",
    "dst_raster = None"
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}