Plot GMMA results.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Base imports for plotting and geospatial data management"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import os\n",
    "import numpy as np\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import cm\n",
    "from matplotlib.collections import LineCollection\n",
    "from matplotlib.patches import Polygon\n",
    "\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "\n",
    "import shapefile\n",
    "\n",
    "# OGR/GDAL imports:\n",
    "from osgeo import gdal\n",
    "from osgeo import osr\n",
    "from osgeo import ogr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from osgeo.gdalconst import GDT_Float32\n",
    "gdal.UseExceptions()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def getProjection(shape_file):\n",
    "    \"\"\"\n",
    "    Get the projection of a shape file\n",
    "    \n",
    "    :param str shape_file: Name of a valid shape file\n",
    "\n",
    "    :returns: :class:`osgeo.osr.SpatialReference` of the shape file\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    driver = ogr.GetDriverByName('ESRI Shapefile')\n",
    "    dataset = driver.Open(shape_file)\n",
    "    layer = dataset.GetLayer()\n",
    "    spatial_ref = layer.GetSpatialRef()\n",
    "    \n",
    "    return spatial_ref"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def getExtent(inputShapeFile):\n",
    "    inDriver = ogr.GetDriverByName(\"ESRI Shapefile\")\n",
    "    inDataSource = inDriver.Open(inputShapeFile, 0)\n",
    "    inLayer = inDataSource.GetLayer()\n",
    "    extent = inLayer.GetExtent()\n",
    "    return extent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def reproject(inputShapeFile, outputShapeFile, outEPSG=4326):\n",
    "    \"\"\"\n",
    "    Reproject a shape file to a known projection. \n",
    "    \n",
    "    :param str inputShapeFile: Source shape file to be reprojected.\n",
    "    :param str outputShapeFile: Destination shape file.\n",
    "    :param int outEPSG: EPSG code for the output projection. Default is \n",
    "                        4326 (WGS 1984)\n",
    "                        \n",
    "    \"\"\"\n",
    "    driver = ogr.GetDriverByName('ESRI Shapefile')\n",
    "    insref = getProjection(inputShapeFile)\n",
    "    outsref = osr.SpatialReference()\n",
    "    outsref.ImportFromEPSG(outEPSG)\n",
    "    \n",
    "    coordTransform = osr.CoordinateTransformation(insref, outsref)\n",
    "    \n",
    "    inDataSet = driver.Open(inputShapeFile)\n",
    "    inLayer = inDataSet.GetLayer()\n",
    "\n",
    "    # create the output layer\n",
    "    if os.path.exists(outputShapeFile):\n",
    "        driver.DeleteDataSource(outputShapeFile)\n",
    "    outDataSet = driver.CreateDataSource(outputShapeFile)\n",
    "    outPrjFile = outputShapeFile.replace('.shp', '.prj')\n",
    "        \n",
    "    fh = open(outPrjFile, 'w')\n",
    "    fh.write(outsref.ExportToWkt())\n",
    "    fh.close()\n",
    "    \n",
    "    outLayer = outDataSet.CreateLayer(\"\", geom_type=ogr.wkbMultiPolygon)\n",
    "\n",
    "    # add fields\n",
    "    inLayerDefn = inLayer.GetLayerDefn()\n",
    "    for i in range(0, inLayerDefn.GetFieldCount()):\n",
    "        fieldDefn = inLayerDefn.GetFieldDefn(i)\n",
    "        outLayer.CreateField(fieldDefn)\n",
    "\n",
    "    # get the output layer's feature definition\n",
    "    outLayerDefn = outLayer.GetLayerDefn()\n",
    "\n",
    "    # loop through the input features\n",
    "    inFeature = inLayer.GetNextFeature()\n",
    "    while inFeature:\n",
    "        # get the input geometry\n",
    "        geom = inFeature.GetGeometryRef()\n",
    "        # reproject the geometry\n",
    "        geom.Transform(coordTransform)\n",
    "        # create a new feature\n",
    "        outFeature = ogr.Feature(outLayerDefn)\n",
    "        # set the geometry and attribute\n",
    "        outFeature.SetGeometry(geom)\n",
    "        for i in range(0, outLayerDefn.GetFieldCount()):\n",
    "            outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))\n",
    "        # add the feature to the shapefile\n",
    "        outLayer.CreateFeature(outFeature)\n",
    "        # destroy the features and get the next input feature\n",
    "        outFeature.Destroy()\n",
    "        inFeature.Destroy()\n",
    "        inFeature = inLayer.GetNextFeature()\n",
    "\n",
    "    # close the shapefiles\n",
    "    inDataSet.Destroy()\n",
    "    outDataSet.Destroy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "inshpfile = \"C:/workspace/gmma/output/event_pop_affect.shp\"\n",
    "outshpfile = \"C:/workspace/data/Glenda/BGY_impact_gcs.shp\"\n",
    "#reproject(inshpfile, outshpfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def drawmap(kwargs):\n",
    "    fig = plt.figure(figsize=(10,10))\n",
    "    ax = plt.subplot(111)\n",
    "    m = Basemap(**kwargs)\n",
    "    m.drawcoastlines()\n",
    "    m.drawmapboundary()\n",
    "    m.drawcountries()\n",
    "    m.fillcontinents('0.75')\n",
    "    \n",
    "    return m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "extent = getExtent(outshpfile)\n",
    "llcrnrlon = int(10*extent[0])/10.\n",
    "urcrnrlon = int(10*extent[1] + 1)/10.\n",
    "llcrnrlat = int(10*extent[2])/10.\n",
    "urcrnrlat = int(10*extent[3] + 1)/10.\n",
    "mapkwargs = dict(projection='cyl', llcrnrlon=llcrnrlon,\n",
    "                  llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon,\n",
    "                  urcrnrlat=urcrnrlat, resolution='f')\n",
    "m = drawmap(mapkwargs)\n",
    "m.readshapefile(\"C:/workspace/data/Glenda/BGY_impact_gcs\", 'results')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "r = shapefile.Reader(outshpfile)\n",
    "shapes = r.shapes()\n",
    "records = r.records()\n",
    "print r.fields\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}