01_Working with vectors and rasters.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vector Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p> Vector data are data in which geographic features are represented as discrete geometries—specifically, points, lines, and polygons. However, continuous data, such as elevation, cannot be considered as vector data </p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### OGR Library"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The OGR Simple Features Library is part of the Geospatial Data Abstraction Library (GDAL), an extremely popular open source library for reading and writing spatial data. The OGR portion of GDAL is the part that provides the ability to read and write many different vector data formats. OGR also allows you to create and manipulate geometries; edit attribute values; filter vector data based on attribute values or spatial location; and it also offers data analysis capabilities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/NEONScience/NEON-Data-Skills/dev-aten/graphics/vector-general/Attribute_Table.png\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import convention for OGR :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import ogr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Drivers\n",
    "Get a driver (ESRI Shapefile):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ESRI Shapefile'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "driver = ogr.GetDriverByName('ESRI shapefile')\n",
    "driver.name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Open and iterate over a shapefile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<osgeo.ogr.DataSource; proxy of <Swig Object of type 'OGRDataSourceShadow *' at 0x7f9db5755b70> >"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "shapefile = 'CAWa_CropType_samples.shp' # found in kaggle\n",
    "ds = ogr.Open(shapefile, 0)\n",
    "if ds is None:\n",
    "    sys.exit('Error opening file')\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Get Layers of a shapefie"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'CAWa_CropType_samples'"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "layer = ds.GetLayer(0)\n",
    "layerName = layer.GetName()\n",
    "layerName"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8435"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "layer.GetFeatureCount()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Layer Extent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(60.201329730000005, 72.35394193117526, 37.42410181284845, 41.82521504999999)\n",
      "Upper-Left x:60.201329730000005, y:41.82521504999999\n",
      "Lower-Right x:72.35394193117526, y:37.42410181284845\n"
     ]
    }
   ],
   "source": [
    "ext = layer.GetExtent()\n",
    "print(ext)\n",
    "print('Upper-Left x:{}, y:{}'.format(ext[0], ext[3]))\n",
    "print('Lower-Right x:{}, y:{}'.format(ext[1], ext[2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Attributes of layer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('sampler', 'String', 80),\n",
       " ('country', 'String', 80),\n",
       " ('region', 'String', 80),\n",
       " ('date', 'String', 80),\n",
       " ('year', 'String', 80),\n",
       " ('label_1', 'String', 80),\n",
       " ('label_2', 'String', 80),\n",
       " ('area', 'Real', 24)]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[(fld.GetName(), ogr.GetFieldTypeName(fld.GetType()), fld.GetWidth()) for fld in layer.schema]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Iterate over the features of layer and work with geometry and values of fields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.696569129837691e-05 POLYGON ((71.5306558053016 40.4970238287256,71.5330325978266 40.4959484127077,71.5345893969305 40.4952977324651,71.531594638349 40.4927491740512,71.5267809646491 40.4946972078007,71.5295730212381 40.4963736786199,71.5300811981627 40.4967867847025,71.5306558053016 40.4970238287256)) cotton\n"
     ]
    }
   ],
   "source": [
    "chosenField = 'label_1'\n",
    "for feature in layer:\n",
    "    value = feature.GetField(chosenField)\n",
    "    geom = feature.geometry()\n",
    "    area = geom.Area()\n",
    "    wkt = geom.ExportToWkt()\n",
    "    print(area,wkt,value)\n",
    "    break\n",
    "layer.ResetReading()\n",
    "ds = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding values to a shapefile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adding new features to existing layers is exactly the same as adding them to empty layers. Create an empty feature based on the layer definition, populate it, and insert it into the layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "layer.CreateField(ogr.FieldDefn('Confidence', ogr.OFTInteger))\n",
    "confidenceLevel = 90\n",
    "for feature in layer:\n",
    "    feature.SetField('Confidence', confidenceLevel)\n",
    "    layer.SetFeature(feature)\n",
    "    break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Raster Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Raster data is a two- or threedimensional array of data values such as the pixels in a photograph. Raster datasets have the potential to include any type of data. Nonetheless, objects that can be thought of as points, lines, or polygons are usually better left as vectors. For example, country boundaries lend themselves perfectly to a polygon vector dataset. This same data could be stored in a raster, but it would take up more disk space and the boundaries wouldn’t be nice, smooth lines. You also couldn’t use vector data analysis functions such as buffering and intersecting. Raster is a perfect choice when values change continuously instead of at sharply defined boundaries. This includes common datasets such as elevation, slope, aspect, precipitation, temperature, and satellite data, but it can include many other things, too."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GDAL Library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import gdal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Translate netCDF to GeoTiff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_filename = \"test.nc\"\n",
    "src_ds = gdal.Open(input_filename)\n",
    "out_format_r = \"GTiff\"\n",
    "driver_r = gdal.GetDriverByName(out_format_r)\n",
    "out_rasterfile = \"test.tif\"\n",
    "dst_ds = driver_r.CreateCopy(out_rasterfile, src_ds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read Raster Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x7f9db5755cf0> >"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "image = 'August012020_B04.tif'\n",
    "ds = gdal.Open(image)\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get image extent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3601552.551268, 4170287.4574389, 3627982.551268, 4147477.4574389)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gt = ds.GetGeoTransform()\n",
    "xmin,xmax,ymin,ymax = (gt[0], gt[3], gt[0] + gt[1] * ds.RasterXSize,gt[3] + gt[5] * ds.RasterYSize)\n",
    "xmin,xmax,ymin,ymax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate Dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(51950, 56874)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import math\n",
    "rows = math.ceil((ymax - ymin) / -gt[5])\n",
    "columns = math.ceil((xmax - xmin) / gt[1])\n",
    "rows,columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read Data of Raster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1384,  985,  745, ...,  719,  411,  411],\n",
       "       [1650, 1252, 1146, ...,  255,  143,  143],\n",
       "       [1650, 1252, 1134, ...,  130,  176,  176],\n",
       "       ...,\n",
       "       [4420, 3760, 3632, ..., 1158, 1532, 1560],\n",
       "       [4420, 3760, 3632, ..., 1420, 1526, 1456],\n",
       "       [3416, 4152, 4204, ..., 1490, 1524, 1396]], dtype=uint16)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = ds.GetRasterBand(1).ReadAsArray()\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "del data\n",
    "del ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use Case: NDVI Calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.4052428  0.52994514 0.6310968  ... 0.6194761  0.75053114 0.75053114]\n",
      " [0.26796806 0.37742418 0.45141217 ... 0.85531914 0.8879749  0.8879749 ]\n",
      " [0.26796806 0.37742418 0.40814197 ... 0.89930284 0.8500852  0.8500852 ]\n",
      " ...\n",
      " [0.11030596 0.1790393  0.18198198 ... 0.4427334  0.2887651  0.26622766]\n",
      " [0.11030596 0.1790393  0.18198198 ... 0.33050448 0.26599327 0.25371605]\n",
      " [0.13780919 0.12992457 0.13604604 ... 0.303087   0.24103586 0.24824986]]\n",
      "Total number of vegetation pixels: 4493164\n"
     ]
    }
   ],
   "source": [
    "red_ds = gdal.Open('August012020_B04.tif')\n",
    "nir_ds = gdal.Open('August012020_B08.tif')\n",
    "red = red_ds.GetRasterBand(1).ReadAsArray().astype(np.float32)\n",
    "nir = nir_ds.GetRasterBand(1).ReadAsArray().astype(np.float32)\n",
    "shape = red.shape\n",
    "\n",
    "check = np.logical_or(red > 0, nir > 0)\n",
    "ndvi = np.where(check, (nir - red) / (nir + red), -9999)\n",
    "print(ndvi)\n",
    "print(\"Total number of vegetation pixels:\",np.count_nonzero(ndvi>0.2))\n",
    "\n",
    "del check\n",
    "del red\n",
    "del nir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save NDVI data as a raster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = red_ds.GetGeoTransform()\n",
    "proj = red_ds.GetProjection()\n",
    "driver = gdal.GetDriverByName(\"GTiff\")\n",
    "dst_ds = driver.Create('ndvi.tif', shape[1],shape[0], 1, gdal.GDT_Float32)\n",
    "dst_ds.SetGeoTransform(geo)\n",
    "dst_ds.SetProjection(proj)\n",
    "dst_ds.GetRasterBand(1).WriteArray(ndvi)\n",
    "\n",
    "del red_ds\n",
    "del nir_ds\n",
    "del dst_ds\n",
    "del ndvi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}