Working_with_geotiffs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some examples for Working with geotiffs\n",
    "\n",
    "Using IOOS3 conda environment\n",
    "\n",
    "A good source for info is http://www.gis.usu.edu/~chrisg/python/2009/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bands:  1 rows:  12265 cols:  13052\n",
      "xOrigin, yOrigin:  376206.5224272934 4625437.967243103\n",
      "Shape of array:  (12265, 13052)\n",
      "upper left x, y:  376206.5224272934 4625437.967243103\n",
      "lower right x, y:  377511.72242729337 4624211.467243103\n"
     ]
    }
   ],
   "source": [
    "from osgeo import gdal, gdalnumeric, ogr, osr\n",
    "import sys\n",
    "gdal.UseExceptions()\n",
    "\n",
    "fn = \"D:/crs/proj/2015_Sandwich/Sandwich_all_surveys_preliminary/2016-03-30_DEM_10cm_cropped.tif\"\n",
    "ds=gdal.Open(fn)\n",
    "rows = ds.RasterYSize\n",
    "cols = ds.RasterXSize\n",
    "bands = ds.RasterCount\n",
    "print(\"bands: \",bands, \"rows: \", rows, \"cols: \",cols)\n",
    "\n",
    "# get georeference info\n",
    "transform = ds.GetGeoTransform()\n",
    "xOrigin = transform[0]\n",
    "yOrigin = transform[3]\n",
    "pixelWidth = transform[1]\n",
    "pixelHeight = transform[5]\n",
    "print(\"xOrigin, yOrigin: \",xOrigin,yOrigin)\n",
    "\n",
    "# read an n-dimensional array into an ndArray Appropriate when RasterCount > 1.\n",
    "nda=ds.ReadAsArray()\n",
    "print(\"Shape of array: \",nda.shape)\n",
    "\n",
    "# different way of accessing transformation info\n",
    "# coordinates are located at the top left corner of the pixels\n",
    "ulx, xres, xskew, uly, yskew, yres  = ds.GetGeoTransform()\n",
    "lrx = ulx + (ds.RasterXSize * xres)\n",
    "lry = uly + (ds.RasterYSize * yres)\n",
    "print(\"upper left x, y: \", ulx,uly)\n",
    "print(\"lower right x, y: \", lrx,lry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pixx, pixy:  7934 14379\n"
     ]
    }
   ],
   "source": [
    "# get the pixel coordinates for some point x, y\n",
    "x = 377000.\n",
    "y = 4624000.\n",
    "xOffset = int((x-xOrigin) / pixelWidth)\n",
    "yOffset = int((y-yOrigin) / pixelHeight)\n",
    "print(\"pixx, pixy: \",xOffset,yOffset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get the projection and print in nice format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PROJCS[\"NAD83 / UTM zone 19N\",\n",
      "    GEOGCS[\"NAD83\",\n",
      "        DATUM[\"North_American_Datum_1983\",\n",
      "            SPHEROID[\"GRS 1980\",6378137,298.2572221010042,\n",
      "                AUTHORITY[\"EPSG\",\"7019\"]],\n",
      "            AUTHORITY[\"EPSG\",\"6269\"]],\n",
      "        PRIMEM[\"Greenwich\",0],\n",
      "        UNIT[\"degree\",0.0174532925199433],\n",
      "        AUTHORITY[\"EPSG\",\"4269\"]],\n",
      "    PROJECTION[\"Transverse_Mercator\"],\n",
      "    PARAMETER[\"latitude_of_origin\",0],\n",
      "    PARAMETER[\"central_meridian\",-69],\n",
      "    PARAMETER[\"scale_factor\",0.9996],\n",
      "    PARAMETER[\"false_easting\",500000],\n",
      "    PARAMETER[\"false_northing\",0],\n",
      "    UNIT[\"metre\",1,\n",
      "        AUTHORITY[\"EPSG\",\"9001\"]],\n",
      "    AUTHORITY[\"EPSG\",\"26919\"]]\n"
     ]
    }
   ],
   "source": [
    "raster_wkt = ds.GetProjection()\n",
    "spatial_ref = osr.SpatialReference()\n",
    "spatial_ref.ImportFromWkt(raster_wkt)\n",
    "print(spatial_ref.ExportToPrettyWkt())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### or in Proj4 format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+proj=utm +zone=19 +datum=NAD83 +units=m +no_defs \n"
     ]
    }
   ],
   "source": [
    "print(spatial_ref.ExportToProj4())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(376206.5224272934, 0.09999999999999822, 0.0, 4625437.967243103, 0.0, -0.10000000000003037)\n"
     ]
    }
   ],
   "source": [
    "# Load the source data as a gdalnumeric array\n",
    "srcArray = gdalnumeric.LoadFile(fn)\n",
    "\n",
    "# Also load as a gdal image to get geotransform\n",
    "# (world file) info\n",
    "srcImage = gdal.Open(fn)\n",
    "geoTrans = srcImage.GetGeoTransform()\n",
    "print(geoTrans)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get some info on a band number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ds_info( band_num, input_file ):\n",
    "    src_ds = gdal.Open( input_file )\n",
    "    if src_ds is None:\n",
    "        print('Unable to open {%s}'.format(input_file))\n",
    "\n",
    "    try:\n",
    "        srcband = src_ds.GetRasterBand(band_num)\n",
    "    except RuntimeError as e:\n",
    "        print('No band {%i} found'.format(band_num))\n",
    "        print(e)\n",
    "\n",
    "    print(\"[ NO DATA VALUE ] = {}\".format(srcband.GetNoDataValue()))\n",
    "    print(\"[ MIN ] = {}\".format(srcband.GetMinimum()))\n",
    "    print(\"[ MAX ] = {}\".format(srcband.GetMaximum()))\n",
    "    print (\"[ SCALE ] = {}\".format(srcband.GetScale()))\n",
    "    print (\"[ UNIT TYPE ] = {}\".format(srcband.GetUnitType()))\n",
    "    ctable = srcband.GetColorTable()\n",
    "\n",
    "    if ctable is None:\n",
    "        print('No ColorTable found')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ NO DATA VALUE ] = -32767.0\n",
      "[ MIN ] = None\n",
      "[ MAX ] = None\n",
      "[ SCALE ] = 1.0\n",
      "[ UNIT TYPE ] = metre\n",
      "No ColorTable found\n"
     ]
    }
   ],
   "source": [
    "ds_info(1,fn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here is an example of a command line that\n",
    "# gdalbuildvrt -separate -te xmin ymin xmax ymax -input_file_list my_filenames.txt output_file.vrt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CompletedProcess(args=['gdalinfo', 'D:/crs/proj/2015_Sandwich/Sandwich_all_surveys_preliminary/2016-03-30_DEM_10cm_cropped.tif'], returncode=0, stdout=b'Driver: GTiff/GeoTIFF\\r\\nFiles: D:/crs/proj/2015_Sandwich/Sandwich_all_surveys_preliminary/2016-03-30_DEM_10cm_cropped.tif\\r\\nSize is 13052, 12265\\r\\nCoordinate System is:\\r\\nPROJCS[\"NAD83 / UTM zone 19N\",\\r\\n    GEOGCS[\"NAD83\",\\r\\n        DATUM[\"North_American_Datum_1983\",\\r\\n            SPHEROID[\"GRS 1980\",6378137,298.2572221010042,\\r\\n                AUTHORITY[\"EPSG\",\"7019\"]],\\r\\n            AUTHORITY[\"EPSG\",\"6269\"]],\\r\\n        PRIMEM[\"Greenwich\",0],\\r\\n        UNIT[\"degree\",0.0174532925199433],\\r\\n        AUTHORITY[\"EPSG\",\"4269\"]],\\r\\n    PROJECTION[\"Transverse_Mercator\"],\\r\\n    PARAMETER[\"latitude_of_origin\",0],\\r\\n    PARAMETER[\"central_meridian\",-69],\\r\\n    PARAMETER[\"scale_factor\",0.9996],\\r\\n    PARAMETER[\"false_easting\",500000],\\r\\n    PARAMETER[\"false_northing\",0],\\r\\n    UNIT[\"metre\",1,\\r\\n        AUTHORITY[\"EPSG\",\"9001\"]],\\r\\n    AUTHORITY[\"EPSG\",\"26919\"]]\\r\\nOrigin = (376206.522427293413784,4625437.967243103310466)\\r\\nPixel Size = (0.099999999999998,-0.100000000000030)\\r\\nMetadata:\\r\\n  AREA_OR_POINT=Point\\r\\nImage Structure Metadata:\\r\\n  INTERLEAVE=BAND\\r\\nCorner Coordinates:\\r\\nUpper Left  (  376206.522, 4625437.967) ( 70d29\\'21.94\"W, 41d46\\'16.08\"N)\\r\\nLower Left  (  376206.522, 4624211.467) ( 70d29\\'21.02\"W, 41d45\\'36.33\"N)\\r\\nUpper Right (  377511.722, 4625437.967) ( 70d28\\'25.43\"W, 41d46\\'16.81\"N)\\r\\nLower Right (  377511.722, 4624211.467) ( 70d28\\'24.52\"W, 41d45\\'37.05\"N)\\r\\nCenter      (  376859.122, 4624824.717) ( 70d28\\'53.23\"W, 41d45\\'56.57\"N)\\r\\nBand 1 Block=13052x1 Type=Float32, ColorInterp=Gray\\r\\n  NoData Value=-32767\\r\\n  Unit Type: metre\\r\\n')\n"
     ]
    }
   ],
   "source": [
    "import subprocess\n",
    "p=subprocess.run(['gdalinfo',fn], stdout=subprocess.PIPE)\n",
    "print(p)"
   ]
  },
  {
   "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}