notebooks/unit06/unit06_06.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "This notebook was prepared by Cayetano Benavent, 2016."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Map algebra with Rasterio and NumPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import rasterio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "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'\n",
    "dst_file = '/tmp/out_mapalgebra_rst.tif'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Explore input datasets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 204.15 314.55\n",
      "1 216.35 312.85\n"
     ]
    }
   ],
   "source": [
    "with rasterio.open(src_file01) as src:\n",
    "    for i in src.indexes:\n",
    "        band = src.read(i)\n",
    "        print(i, band.min(), band.max())\n",
    "\n",
    "with rasterio.open(src_file02) as src:\n",
    "    for i in src.indexes:\n",
    "        band = src.read(i)\n",
    "        print(i, band.min(), band.max())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read bands as NumPy arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'numpy.ndarray'>\n",
      "(721, 1440) 2 1038240\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",
      " [ 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   219.3500061\n",
      "   219.3500061 ]\n",
      " [ 219.75001526  219.75001526  219.75001526 ...,  219.75001526\n",
      "   219.75001526  219.75001526]]\n"
     ]
    }
   ],
   "source": [
    "with rasterio.open(src_file01) as src:\n",
    "    for i in src.indexes:\n",
    "        band = src.read(i)\n",
    "        print(type(band))\n",
    "        print(band.shape, band.ndim, band.size)\n",
    "        print(band)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 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   219.3500061\n",
      "   219.3500061 ]\n",
      " [ 219.75001526  219.75001526  219.75001526 ...,  219.75001526\n",
      "   219.75001526  219.75001526]]\n",
      "<class 'numpy.ndarray'>\n",
      "[[ 269.8500061   269.8500061   269.8500061  ...,  269.8500061   269.8500061\n",
      "   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]]\n",
      "<class 'numpy.ndarray'>\n"
     ]
    }
   ],
   "source": [
    "with rasterio.open(src_file01) as src01:\n",
    "    with rasterio.open(src_file02) as src02:\n",
    "        for i in src01.indexes:\n",
    "            band_src01 = src01.read(i)\n",
    "            print(band_src01)\n",
    "            print(type(band_src01))\n",
    "        for i in src02.indexes:\n",
    "            band_src02 = src02.read(i)\n",
    "            print(band_src02)\n",
    "            print(type(band_src01))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Operating with bands as NumPy arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.70001221  0.70001221  0.70001221 ...,  0.70001221  0.70001221\n",
      "   0.70001221]\n",
      " [ 0.6000061   0.6000061   0.6000061  ...,  0.6000061   0.6000061\n",
      "   0.6000061 ]\n",
      " [ 0.6000061   0.6000061   0.5        ...,  0.6000061   0.6000061\n",
      "   0.6000061 ]\n",
      " ..., \n",
      " [ 0.29998779  0.29998779  0.29998779 ...,  0.29998779  0.29998779\n",
      "   0.29998779]\n",
      " [ 0.3999939   0.3999939   0.3999939  ...,  0.3999939   0.3999939\n",
      "   0.3999939 ]\n",
      " [ 0.30000305  0.30000305  0.30000305 ...,  0.30000305  0.30000305\n",
      "   0.30000305]]\n"
     ]
    }
   ],
   "source": [
    "with rasterio.open(src_file01) as src01:\n",
    "    with rasterio.open(src_file02) as src02:\n",
    "        for i in src01.indexes:\n",
    "            band_src01 = src01.read(i)\n",
    "        for i in src02.indexes:\n",
    "            band_src02 = src02.read(i)\n",
    "        res = (band_src01 -273.15) - (band_src02 -273.15)\n",
    "        print(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computing with NumPy and storing results on a new raster dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.4/dist-packages/rasterio/__init__.py:99: FutureWarning: GDAL-style transforms are deprecated and will not be supported in Rasterio 1.0.\n",
      "  transform = guard_transform(transform)\n"
     ]
    }
   ],
   "source": [
    "with rasterio.open(src_file01) as src01:\n",
    "    with rasterio.open(src_file02) as src02:\n",
    "        for i in src01.indexes:\n",
    "            band_src01 = src01.read(i)\n",
    "        for i in src02.indexes:\n",
    "            band_src02 = src02.read(i)\n",
    "        res = (band_src01 -273.15) - (band_src02 -273.15)\n",
    "        kwargs = src.meta\n",
    "        with rasterio.open(dst_file, 'w', **kwargs) as dst:\n",
    "            dst.write_band(1, res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Key 'transform' of kwargs dict is a GDAL-style transform (raise a warning):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'transform': (-180.0, 0.25, 0.0, 90.0, 0.0, -0.2498266296809986), 'count': 1, 'crs': {'b': 6371229, 'a': 6371229, 'no_defs': True, 'proj': 'longlat', 'wktext': True}, 'affine': Affine(0.25, 0.0, -180.0,\n",
      "       0.0, -0.2498266296809986, 90.0), 'dtype': 'float32', 'driver': 'GTiff', 'width': 1440, 'nodata': -3.4028234663852886e+38, 'height': 721}\n"
     ]
    }
   ],
   "source": [
    "print(kwargs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we remove 'transform' key from kwargs dict, the warning disappear. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "with rasterio.open(src_file01) as src01:\n",
    "    with rasterio.open(src_file02) as src02:\n",
    "        for i in src01.indexes:\n",
    "            band_src01 = src01.read(i)\n",
    "        for i in src02.indexes:\n",
    "            band_src02 = src02.read(i)\n",
    "        res = (band_src01 -273.15) - (band_src02 -273.15)\n",
    "        \n",
    "        kwargs = src.meta\n",
    "        if kwargs.get('transform'):\n",
    "            kwargs.pop('transform')\n",
    "        \n",
    "        with rasterio.open(dst_file, 'w', **kwargs) as dst:\n",
    "            dst.write_band(1, res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'count': 1, 'crs': {'b': 6371229, 'a': 6371229, 'no_defs': True, 'proj': 'longlat', 'wktext': True}, 'affine': Affine(0.25, 0.0, -180.0,\n",
      "       0.0, -0.2498266296809986, 90.0), 'dtype': 'float32', 'driver': 'GTiff', 'width': 1440, 'nodata': -3.4028234663852886e+38, 'height': 721}\n"
     ]
    }
   ],
   "source": [
    "print(kwargs)"
   ]
  }
 ],
 "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
}