raster_stats.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get Point WKT from a point shapefile using ogr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "POINT (624501.568181818 5628082.59695052)\n",
      "POINT (656596.878308401 5630027.76726122)\n",
      "POINT (694388.75863061 5628499.41915995)\n",
      "POINT (718425.506041427 5639336.79660529)\n",
      "POINT (725233.602128884 5664901.89211738)\n"
     ]
    }
   ],
   "source": [
    "from osgeo import ogr\n",
    "import os\n",
    "\n",
    "shapefile = \"../points.shp\"\n",
    "driver = ogr.GetDriverByName(\"ESRI Shapefile\")\n",
    "dataSource = driver.Open(shapefile, 0)\n",
    "layer = dataSource.GetLayer()\n",
    "\n",
    "for feature in layer:\n",
    "    geom = feature.GetGeometryRef()\n",
    "    print geom.Centroid().ExportToWkt()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use rasterstats (pip install rasterstats) to get the value at a point location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.6670739601440762]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from rasterstats import point_query\n",
    "\n",
    "point = \"POINT (725233.602128884 5664901.89211738)\"\n",
    "point_query(point, \"../NDVI_example_1.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Loop over a point shapefile and extract values on a single band raster (but a multi band one could be used)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.40849066211442864]\n",
      "[0.652933816283167]\n",
      "[-0.4690519447620269]\n",
      "[0.3553161198672801]\n",
      "[0.6670739601440762]\n"
     ]
    }
   ],
   "source": [
    "#from osgeo import ogr\n",
    "#import os\n",
    "\n",
    "shapefile = \"D:/rasterStats/points.shp\"\n",
    "driver = ogr.GetDriverByName(\"ESRI Shapefile\")\n",
    "dataSource = driver.Open(shapefile, 0)\n",
    "layer = dataSource.GetLayer()\n",
    "\n",
    "for feature in layer:\n",
    "    geom = feature.GetGeometryRef()\n",
    "    point = geom.Centroid().ExportToWkt()\n",
    "    print point_query(point, \"..NDVI_example_1.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get mean and max values from a series of polygons from a shapefile with multiple features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['count', 'max', 'mean', 'min']\n",
      "[0.17415108794477352, 0.6178150208015752, 0.2820635024507548, 0.4951216243519921]\n",
      "[0.3866455554962158, 0.683485209941864, 0.45673757791519165, 0.6485721468925476]\n"
     ]
    }
   ],
   "source": [
    "from rasterstats import zonal_stats\n",
    "stats = zonal_stats(\"../polys.shp\",  \"../NDVI_example_1.tif\")\n",
    "print stats[1].keys()\n",
    "#['count', 'min', 'max', 'mean']\n",
    "print [f['mean'] for f in stats]\n",
    "print [f['max'] for f in stats]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Walk through all your .tif images in a directory and extract values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "filename:  NDVI_example_1.tif\n",
      "[0.17415108794477352, 0.6178150208015752, 0.2820635024507548, 0.4951216243519921]\n",
      "filename:  NDVI_example_2.tif\n",
      "[0.7579913160863959, 0.6720901117092226, 0.7985192637380338, 0.6334189456158745]\n"
     ]
    }
   ],
   "source": [
    "for root, dirs, files in os.walk(\"../\"):\n",
    "    for file in files:\n",
    "        if file.endswith(\".tif\"):\n",
    "            raster = (os.path.join(root, file))\n",
    "            stats = zonal_stats(\"../polys.shp\",  raster)\n",
    "            print \"filename: \", file\n",
    "            print [f['mean'] for f in stats]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "my_project",
   "language": "python",
   "name": "my_project"
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}