examples/warp radar.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: Warp radar image\n",
    "\n",
    "The [National Weather Service radar mosaic](http://radar.weather.gov/Conus/full.php) shows the 48 Continental states of the US. Unfortunately, it's in unprojected lat/lon coordindates. In this example, we'll use [rasterio](https://github.com/mapbox/rasterio) to reproject the image to a more standard equal-area projection.\n",
    "\n",
    "Based on [a blog post using GDALwarp]((http://egb13.net/2009/03/using-gdalwarp) and the [rasterio reproject example](https://github.com/mapbox/rasterio/blob/master/examples/reproject.py)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "import os\n",
    "from IPython.display import Image\n",
    "image_location = 'http://radar.weather.gov/Conus/RadarImg/latest.gif'\n",
    "# uncomment the following to use local file (note: absolute paths don't work for Image())\n",
    "#image_location = os.path.join(os.path.abspath('..'), 'data', 'latest.gif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"http://radar.weather.gov/Conus/RadarImg/latest.gif\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Image(url='http://radar.weather.gov/Conus/RadarImg/latest.gif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import rasterio\n",
    "from rasterio.warp import reproject, RESAMPLING\n",
    "from pyproj import Proj"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately, the GIF image does not have geographic metadata, so we will need to set the input projection. We get the transform parameters from the world file at http://radar.weather.gov/Conus/RadarImg/latest_radaronly.gfw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "src_crs = {'init': 'EPSG:4326'}\n",
    "west = -127.620375523875420\n",
    "north = 50.406626367301044\n",
    "dx = 0.017971305190311\n",
    "dy = -dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Desired width of the output image (height will be calculated):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "width = 1600"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to set the output projection. We'll use a Lambert Equal-Area projection, specifically [EPSG 2163](http://epsg.io/2163), which is used by the [National Map](http://nationalmap.gov/) to display the continental US.\n",
    "\n",
    "Here's an example that displays average annual precipication for the CONUS:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"http://nationalmap.gov/small_scale/printable/images/preview/precip/pageprecip_us3.gif\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Image(url='http://nationalmap.gov/small_scale/printable/images/preview/precip/pageprecip_us3.gif')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use pyproj to transform the corners of the radar image to get the bounds to use for our output grid. First we'll compute the south and east boundaries of the map in lat/lon space.  To do that, we'll open the remote file with rasterio to get the image size.\n",
    "\n",
    "Next, we'll compute the corners in projected coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with rasterio.drivers():\n",
    "    with rasterio.open(image_location) as src:\n",
    "        south = north + src.height * dy\n",
    "        east = west + src.width * dx\n",
    "        src_transform = rasterio.transform.from_bounds(west, south, east, north, src.width, src.height)\n",
    "\n",
    "dst_crs = {'init': 'EPSG:2163'}\n",
    "us_equal_area = Proj(**dst_crs)\n",
    "left, bottom = us_equal_area(west, south)\n",
    "right, _ = us_equal_area(east, south)\n",
    "_, top = us_equal_area(east, north)\n",
    "height = width * (top - bottom) / (right - left)\n",
    "dst_transform = rasterio.transform.from_bounds(left, bottom, right, top, width, height)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll initialize a NumPy array of bytes to transform the image data into."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dst_shape = (height, width)\n",
    "destination = np.zeros(dst_shape, np.uint8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now for the actual transformation. We open the radar file with rasterio, and use `rasterio.warp.reproject` to transform it to our new coordinate system. It's important to use `resampling=RESAMPLING.nearest` in this case, because we don't want to interpolate values from the GIF image. For continuous data, other resampling methods may be appropriate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with rasterio.drivers():\n",
    "    with rasterio.open(image_location) as src:\n",
    "        data = src.read(1)\n",
    "        cmap = src.colormap(1)\n",
    "        reproject(data, destination,\n",
    "                  src_transform=src_transform, src_crs=src_crs,\n",
    "                  dst_crs=dst_crs, dst_transform=dst_transform,\n",
    "                  resampling=RESAMPLING.nearest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the result is in a NumPy array. The values in the array are the pixel values from the GIF image, which will not make sense without that file's colormap."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use rasterio to write the output to a new local GIF file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with rasterio.open('warped.gif', 'w', driver='GIF',\n",
    "                   width=width, height=height,\n",
    "                   count=1, dtype=np.uint8) as dst:\n",
    "    dst.write_band(1, destination)\n",
    "    dst.write_colormap(1, cmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"warped.gif\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Image(url='warped.gif')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we've cut off the bottom of the image because we calculated the bounding box from the corners. We could fit it in by extending the lower boundary of the output box. Alternatively, we might want to crop the image more tightly to the 48 states. Try adjusting the values of `left`, `right`, `top`, and `bottom` above that go into the output projection.\n",
    "\n",
    "## Further exploration\n",
    "\n",
    "- Try warping the map to a different map projection. One common example (ubiquitous on the web, though unloved by cartographers) is [EPSG 3857](https://epsg.io/3857), web Mercator.\n",
    "\n",
    "- The National Weather Service also provides a [radar_only image](http://radar.weather.gov/Conus/RadarImg/latest_radaronly.gif) as a transparent GIF without the basemap. Try plotting your own map as a base layer and using the radar image as an overlay."
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}