{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n=================================\nParsing rasterio's geocoordinates\n=================================\n\n\nConverting a projection's cartesian coordinates into 2D longitudes and\nlatitudes.\n\nThese new coordinates might be handy for plotting and indexing, but it should\nbe kept in mind that a grid which is regular in projection coordinates will\nlikely be irregular in lon/lat. It is often recommended to work in the data's\noriginal map projection (see `recipes.rasterio_rgb`).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import os\nimport urllib.request\n\nimport cartopy.crs as ccrs\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom rasterio.warp import transform\n\nimport xarray as xr\n\n# Download the file from rasterio's repository\nurl = 'https://github.com/mapbox/rasterio/raw/master/tests/data/RGB.byte.tif'\nurllib.request.urlretrieve(url, 'RGB.byte.tif')\n\n# Read the data\nda = xr.open_rasterio('RGB.byte.tif')\n\n# Compute the lon/lat coordinates with rasterio.warp.transform\nny, nx = len(da['y']), len(da['x'])\nx, y = np.meshgrid(da['x'], da['y'])\n\n# Rasterio works with 1D arrays\nlon, lat = transform(da.crs, {'init': 'EPSG:4326'},\n                     x.flatten(), y.flatten())\nlon = np.asarray(lon).reshape((ny, nx))\nlat = np.asarray(lat).reshape((ny, nx))\nda.coords['lon'] = (('y', 'x'), lon)\nda.coords['lat'] = (('y', 'x'), lat)\n\n# Compute a greyscale out of the rgb image\ngreyscale = da.mean(dim='band')\n\n# Plot on a map\nax = plt.subplot(projection=ccrs.PlateCarree())\ngreyscale.plot(ax=ax, x='lon', y='lat', transform=ccrs.PlateCarree(),\n               cmap='Greys_r', add_colorbar=False)\nax.coastlines('10m', color='r')\nplt.show()\n\n# Delete the file\nos.remove('RGB.byte.tif')"
      ]
    }
  ],
  "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.6.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}