Ch_1 Raster_Vector_datasets_exploration.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import ogr\n",
    "import gdal\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open a GDAL raster dataset\n",
    "raster = gdal.Open(\"F:/works/practical/S2A/cliped/stack.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'AREA_OR_POINT': 'Area'}\n"
     ]
    }
   ],
   "source": [
    "# Does the raster have metadata?\n",
    "print (raster.GetMetadata())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of bands in image: 4\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# How many bands in the image have?\n",
    "num_bands = raster.RasterCount\n",
    "print('Number of bands in image: {n}\\n'.format(n=num_bands))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image size is: 1770 rows x 2607 columns\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# How many rows and columns?\n",
    "rows = raster.RasterYSize\n",
    "cols = raster.RasterXSize\n",
    "print('Image size is: {r} rows x {c} columns\\n'.format(r=rows, c=cols))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Raster description: F:/works/practical/S2A/cliped/stack.tif\n",
      "Raster metadata:\n",
      "{'AREA_OR_POINT': 'Area'}\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Does the raster have a description or metadata?\n",
    "desc = raster.GetDescription()\n",
    "metadata = raster.GetMetadata()\n",
    "\n",
    "print('Raster description: {desc}'.format(desc=desc))\n",
    "print('Raster metadata:')\n",
    "print(metadata)\n",
    "print('\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Raster driver: GTiff\n",
      "\n"
     ]
    }
   ],
   "source": [
    "driver = raster.GetDriver()\n",
    "print('Raster driver: {d}\\n'.format(d=driver.ShortName))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image projection:\n",
      "PROJCS[\"WGS_1984_UTM_Zone_43N\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",75],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AUTHORITY[\"EPSG\",\"32643\"]]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# What is the raster's projection?\n",
    "proj = raster.GetProjection()\n",
    "print('Image projection:')\n",
    "print(proj + '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image geo-transform: (379903.94753758144, 10.0, 0.0, 3509518.713358881, 0.0, -10.0)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# What is the raster's \"geo-transform\"\n",
    "gt = raster.GetGeoTransform()\n",
    "print('Image geo-transform: {gt}\\n'.format(gt=gt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<osgeo.gdal.Band; proxy of <Swig Object of type 'GDALRasterBandShadow *' at 0x000001A49D08E300> >\n"
     ]
    }
   ],
   "source": [
    "# Access the single (blue) band in the  stacked image\n",
    "blue = raster.GetRasterBand(1)\n",
    "\n",
    "print(blue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Band datatype: 6\n"
     ]
    }
   ],
   "source": [
    "# What is the band's datatype?\n",
    "datatype = blue.DataType\n",
    "print('Band datatype: {dt}'.format(dt=blue.DataType))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Band datatype: UInt32\n"
     ]
    }
   ],
   "source": [
    "datatype_name = gdal.GetDataTypeName(blue.DataType)\n",
    "print('Band datatype: {dt}'.format(dt=datatype_name))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Band datatype size: 32 bytes\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# We can also ask how much space does this datatype take up\n",
    "bytes = gdal.GetDataTypeSize(blue.DataType)\n",
    "print('Band datatype size: {b} bytes\\n'.format(b=bytes))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Band range: 0.33660000562668 - 0.015299990773201\n",
      "Band mean, stddev: 0.049093857786703, 0.014083102996686\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# blue band statistics i.e. mean , median, min, max values?\n",
    "band_max, band_min, band_mean, band_stddev = blue.GetStatistics(0, 1)\n",
    "print('Band range: {minimum} - {maximum}'.format(maximum=band_max,minimum=band_min))\n",
    "print('Band mean, stddev: {m}, {s}\\n'.format(m=band_mean, s=band_stddev))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-999. -999. -999. ... -999. -999. -999.]\n",
      " [-999. -999. -999. ... -999. -999. -999.]\n",
      " [-999. -999. -999. ... -999. -999. -999.]\n",
      " ...\n",
      " [-999. -999. -999. ... -999. -999. -999.]\n",
      " [-999. -999. -999. ... -999. -999. -999.]\n",
      " [-999. -999. -999. ... -999. -999. -999.]]\n",
      "\n",
      "Blue band mean is: -27.97496795654297\n",
      "Size is: (1770, 2607)\n"
     ]
    }
   ],
   "source": [
    "# read the band in the form of numpyarray\n",
    "blue_data = blue.ReadAsArray()\n",
    "\n",
    "print(blue_data)\n",
    "print()\n",
    "print('Blue band mean is: {m}'.format(m=blue_data.mean()))\n",
    "print('Size is: {sz}'.format(sz=blue_data.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  ...\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]]\n",
      "\n",
      " [[-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  ...\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]]\n",
      "\n",
      " [[-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  ...\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]]\n",
      "\n",
      " ...\n",
      "\n",
      " [[-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  ...\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]]\n",
      "\n",
      " [[-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  ...\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]]\n",
      "\n",
      " [[-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  ...\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]\n",
      "  [-999. -999. -999. -999.]]]\n",
      "float64\n"
     ]
    }
   ],
   "source": [
    "# Initialize a 3d array -- use the size properties of our image for portability!\n",
    "image = np.zeros((raster.RasterYSize, raster.RasterXSize, raster.RasterCount))\n",
    "\n",
    "# Loop over all bands in dataset\n",
    "for b in range(raster.RasterCount):\n",
    "    # Remember, GDAL index is on 1, but Python is on 0 -- so we add 1 for our GDAL calls\n",
    "    band = raster.GetRasterBand(b + 1)\n",
    "    \n",
    "    # Read in the band's data into the third dimension of our array\n",
    "    image[:, :, b] = band.ReadAsArray()\n",
    "\n",
    "print(image)\n",
    "print(image.dtype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raster.GetRasterBand(1).DataType\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compare datatypes: \n",
      "    when unspecified: float64\n",
      "    when specified: float32\n",
      "Red band mean: -27.942620814044673\n",
      "NIR band mean: -27.87640348750272\n"
     ]
    }
   ],
   "source": [
    "\n",
    "from osgeo import gdal_array\n",
    "\n",
    "# DataType is a property of the individual raster bands\n",
    "image_datatype = raster.GetRasterBand(1).DataType\n",
    "\n",
    "# Allocate our array, but in a more efficient way\n",
    "image_correct = np.zeros((raster.RasterYSize, raster.RasterXSize, raster.RasterCount),\n",
    "                 dtype=gdal_array.GDALTypeCodeToNumericTypeCode(image_datatype))\n",
    "\n",
    "# Loop over all bands in dataset\n",
    "for b in range(raster.RasterCount):\n",
    "    # Remember, GDAL index is on 1, but Python is on 0 -- so we add 1 for our GDAL calls\n",
    "    band = raster.GetRasterBand(b + 1)\n",
    "    \n",
    "    # Read in the band's data into the third dimension of our array\n",
    "    image_correct[:, :, b] = band.ReadAsArray()\n",
    "\n",
    "print(\"Compare datatypes: \")\n",
    "print(\"    when unspecified: {dt}\".format(dt=image.dtype))\n",
    "print(\"    when specified: {dt}\".format(dt=image_correct.dtype))\n",
    "print('Red band mean: {r}'.format(r=image[:, :, 2].mean()))\n",
    "print('NIR band mean: {nir}'.format(nir=image[:, :, 3].mean()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0. -0. -0. ... -0. -0. -0.]\n",
      " [-0. -0. -0. ... -0. -0. -0.]\n",
      " [-0. -0. -0. ... -0. -0. -0.]\n",
      " ...\n",
      " [-0. -0. -0. ... -0. -0. -0.]\n",
      " [-0. -0. -0. ... -0. -0. -0.]\n",
      " [-0. -0. -0. ... -0. -0. -0.]]\n",
      "Maximum value of NDVI = 0.908464587950874 \n"
     ]
    }
   ],
   "source": [
    "b_red = 2\n",
    "b_nir = 3\n",
    "\n",
    "ndvi = (image[:, :, b_nir] - image[:, :, b_red]) / (image[:, :, b_red] + image[:, :, b_nir])\n",
    "\n",
    "print(ndvi)\n",
    "print('Maximum value of NDVI = {ndvi} '.format(ndvi=ndvi.max()) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NDVI matrix: \n",
      "[[-0. -0. -0. ... -0. -0. -0.]\n",
      " [-0. -0. -0. ... -0. -0. -0.]\n",
      " [-0. -0. -0. ... -0. -0. -0.]\n",
      " ...\n",
      " [-0. -0. -0. ... -0. -0. -0.]\n",
      " [-0. -0. -0. ... -0. -0. -0.]\n",
      " [-0. -0. -0. ... -0. -0. -0.]]\n",
      "\n",
      "Max NDVI: 0.908464587950874\n",
      "Mean NDVI: 0.2568754514932775\n",
      "Median NDVI: 0.19104636375367037\n",
      "Min NDVI: -0.4062500429942998\n"
     ]
    }
   ],
   "source": [
    "ndvi = (image[:, :, b_nir] - image[:, :, b_red]) / \\\n",
    "        (image[:, :, b_nir] + image[:, :, b_red]).astype(np.float64)\n",
    "\n",
    "print('NDVI matrix: ')\n",
    "print(ndvi)\n",
    "\n",
    "print('\\nMax NDVI: {m}'.format(m=ndvi.max()))\n",
    "print('Mean NDVI: {m}'.format(m=ndvi.mean()))\n",
    "print('Median NDVI: {m}'.format(m=np.median(ndvi)))\n",
    "print('Min NDVI: {m}'.format(m=ndvi.min()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dataset driver is: ESRI Shapefile\n",
      "\n",
      "Number of layers = 1\n",
      "\n",
      "The layer is named: Punjab Districts\n",
      "\n",
      "The layer's geometry is: Polygon\n",
      "\n",
      "Layer projection is: +proj=longlat +datum=WGS84 +no_defs \n",
      "\n",
      "Layer has 36 features\n",
      "\n",
      "Layer has 2 fields\n",
      "Their names are: \n",
      "\tOBJECTID - Integer64\n",
      "\tDistt_name - String\n"
     ]
    }
   ],
   "source": [
    "### Let's get starting with the vector dataset\n",
    "dataset = ogr.Open('F:/works/practical/Dist-Boundary/Punjab Districts.shp')\n",
    "\n",
    "# the driver of the shapefile \n",
    "driver = dataset.GetDriver()\n",
    "print('Dataset driver is: {n}\\n'.format(n=driver.name))\n",
    "\n",
    "# How many layers contain in shapefile\n",
    "layers = dataset.GetLayerCount()\n",
    "print('Number of layers = {n}\\n'.format(n=layers))\n",
    "\n",
    "### What is the name of the 1 layer?\n",
    "layer = dataset.GetLayerByIndex(0)\n",
    "print('The layer is named: {n}\\n'.format(n=layer.GetName()))\n",
    "\n",
    "geometry = layer.GetGeomType()\n",
    "\n",
    "# So we need to translate it to the name of the enum\n",
    "geometry_name = ogr.GeometryTypeToName(geometry)\n",
    "print(\"The layer's geometry is: {geom}\\n\".format(geom=geometry_name))\n",
    "\n",
    "# Get the spatial reference\n",
    "spatial_ref = layer.GetSpatialRef()\n",
    "\n",
    "# Export this spatial reference to something we can read... like the Proj4\n",
    "proj4 = spatial_ref.ExportToProj4()\n",
    "print('Layer projection is: {proj4}\\n'.format(proj4=proj4))\n",
    "\n",
    "### How many features are in the layer?\n",
    "feature_count = layer.GetFeatureCount()\n",
    "print('Layer has {n} features\\n'.format(n=feature_count))\n",
    "\n",
    "### How many fields are in the shapefile, and what are their names?\n",
    "# First we need to capture the layer definition\n",
    "defn = layer.GetLayerDefn()\n",
    "\n",
    "# How many fields\n",
    "field_count = defn.GetFieldCount()\n",
    "print('Layer has {n} fields'.format(n=field_count))\n",
    "\n",
    "# What are their names?\n",
    "print('Their names are: ')\n",
    "for i in range(field_count):\n",
    "    field_defn = defn.GetFieldDefn(i)\n",
    "    print('\\t{name} - {datatype}'.format(name=field_defn.GetName(),\n",
    "                                         datatype=field_defn.GetTypeName()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "     ### End chapter 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.6"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}