notebooks/unit05/unit05_09.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook was prepared by Cayetano Benavent, 2016."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reading and writing raster data with Rasterio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Importing library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import rasterio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Translate datasets in Rasterio\n",
    "\n",
    "Translate example: NetCDF datset to Geotiff."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "input_file = \"../../data/mslp_gfs/mslp_europe.nc\"\n",
    "out_tif_file = \"/tmp/mslp_europe_rst.tif\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run translate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with rasterio.drivers():\n",
    "    rasterio.copy(input_file, out_tif_file, driver='GTiff')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reading raster data and printing basic metadata:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "- Driver: netCDF\n",
      "\n",
      "- Width 320, Height 180\n",
      "\n",
      "- Shape: (180, 320)\n",
      "\n",
      "- CRS (Proj4): \n",
      "{'a': 6367470, 'proj': 'longlat', 'no_defs': True, 'b': 6367470}\n",
      "\n",
      "- CRS (WKT): \n",
      "GEOGCS[\"unnamed ellipse\",DATUM[\"unknown\",SPHEROID[\"unnamed\",6367470,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]]\n",
      "\n",
      "- Affine transf. matrix: \n",
      "| 0.25, 0.00,-30.12|\n",
      "| 0.00,-0.25, 75.12|\n",
      "| 0.00, 0.00, 1.00|\n",
      "\n",
      "- Band count: 1\n",
      "\n",
      "- Band indexes: [1]\n",
      "\n",
      "- Metadata summary: \n",
      "{'height': 180, 'dtype': 'float64', 'count': 1, 'crs': {'a': 6367470, 'proj': 'longlat', 'no_defs': True, 'b': 6367470}, 'affine': Affine(0.25, 0.0, -30.125,\n",
      "       0.0, -0.25, 75.125), 'width': 320, 'transform': (-30.125, 0.25, 0.0, 75.125, 0.0, -0.25), 'driver': 'netCDF', 'nodata': 9.969209968386869e+36}\n",
      "\n",
      "- GDAL tags: \n",
      "{'crs#semi_major_axis': '6367470', 'crs#GeoTransform': '-30.125 0.25 0 75.125 0 -0.25 ', 'Band1#GRIB_SHORT_NAME': '0-SFC', 'Band1#grid_mapping': 'crs', 'Band1#GRIB_FORECAST_SECONDS': '0 sec', 'crs#grid_mapping_name': 'latitude_longitude', 'crs#inverse_flattening': '0', 'Band1#GRIB_REF_TIME': '   283996800 sec UTC', 'crs#longitude_of_prime_meridian': '0', 'lon#long_name': 'longitude', 'crs#spatial_ref': 'GEOGCS[\"unnamed ellipse\",DATUM[\"unknown\",SPHEROID[\"unnamed\",6367470,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]]', 'Band1#_FillValue': '9.969209968386869e+36', 'lat#long_name': 'latitude', 'Band1#GRIB_UNIT': '[Pa]', 'NC_GLOBAL#history': 'Sat Apr 02 22:54:59 2016: GDAL CreateCopy( /tmp/mslp_europe.nc, ... )', 'Band1#GRIB_ELEMENT': 'MSL', 'Band1#GRIB_VALID_TIME': '   283996800 sec UTC', 'NC_GLOBAL#Conventions': 'CF-1.5', 'lat#units': 'degrees_north', 'Band1#long_name': 'GDAL Band Number 1', 'lon#units': 'degrees_east', 'lat#standard_name': 'latitude', 'Band1#GRIB_COMMENT': 'Mean sea level pressure [Pa]', 'lon#standard_name': 'longitude', 'NC_GLOBAL#GDAL': 'GDAL 1.11.2, released 2015/02/10'}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "with rasterio.drivers():\n",
    "    with rasterio.open(input_file) as src:\n",
    "        print(\"- Driver: {}\\n\".format(src.driver))\n",
    "        print(\"- Width {}, Height {}\\n\".format(src.width, src.height))\n",
    "        print(\"- Shape: {}\\n\".format(src.shape))\n",
    "        print(\"- CRS (Proj4): \\n{}\\n\".format(src.crs))\n",
    "        print(\"- CRS (WKT): \\n{}\\n\".format(src.crs_wkt))\n",
    "        print(\"- Affine transf. matrix: \\n{}\\n\".format(src.affine))\n",
    "        print(\"- Band count: {}\\n\".format(src.count))\n",
    "        print(\"- Band indexes: {}\\n\".format(src.indexes))\n",
    "        print(\"- Metadata summary: \\n{}\\n\".format(src.meta))\n",
    "        print(\"- GDAL tags: \\n{}\\n\".format(src.tags()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Affine transformation in Rasterio\n",
    "In the above result we can see affine transformation matrix which is described below:\n",
    "```\n",
    "a: Scale X \n",
    "e: Scale Y\n",
    "d: Rotation (shear in Y direction)\n",
    "b: Rotation (shear in X direction)\n",
    "c: Translation X\n",
    "f: Translation Y\n",
    "\n",
    "| x' |   | a b c | | x |\n",
    "| y' | = | d e f | | y |\n",
    "| 1  |   | 0 0 1 | | 1 |\n",
    "\n",
    "\n",
    "GDAL coefficients order in GeoTransform:\n",
    "c, a, b, f, d, e\n",
    "```\n",
    "More info:\n",
    "\n",
    "https://github.com/sgillies/affine\n",
    "\n",
    "https://en.wikipedia.org/wiki/File:2D_affine_transformation_matrix.svg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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
}