area1_classification.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee\n",
    "ee.Initialize()\n",
    "from IPython.display import display, Image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scale = 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roi = ee.FeatureCollection(\"users/nkeikon/myanmar_sr/area1\")\n",
    "op = ee.FeatureCollection(\"users/nkeikon/myanmar_sr/OP_young\")\n",
    "PA = ee.FeatureCollection(\"WCMC/WDPA/current/polygons\")\n",
    "s2TNI = ee.Image(\"users/nkeikon/myanmar_sr/s2image\").clip(roi)\n",
    "palm = ee.FeatureCollection(\"users/nkeikon/myanmar_sr/palm_area1\")\n",
    "rubber = ee.FeatureCollection(\"users/nkeikon/myanmar_sr/rubber_area1\")\n",
    "other = ee.FeatureCollection(\"users/nkeikon/myanmar_sr/other_area1\")\n",
    "bare = ee.FeatureCollection(\"users/nkeikon/myanmar_sr/bare_area1\")\n",
    "water = ee.FeatureCollection(\"users/nkeikon/myanmar_sr/water_area1\")\n",
    "shrub = ee.FeatureCollection(\"users/nkeikon/myanmar_sr/shrub_area1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image(\n",
    "    url=s2TNI.getThumbUrl(\n",
    "        {\"bands\": \"red_mean,green_mean,blue_mean\", \"min\": 0, \"max\": 0.4, \"size\": \"400\"}\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add Sentinel-1 data (only ascending for Jan to Mar 2018 due to noise)\n",
    "s1_JanMar18 = (\n",
    "    ee.ImageCollection(\"COPERNICUS/S1_GRD\")\n",
    "    .filter(ee.Filter.eq(\"instrumentMode\", \"IW\"))\n",
    "    .filter(ee.Filter.listContains(\"transmitterReceiverPolarisation\", \"VH\"))\n",
    "    .filter(ee.Filter.listContains(\"transmitterReceiverPolarisation\", \"VV\"))\n",
    "    .filter(ee.Filter.eq(\"orbitProperties_pass\", \"ASCENDING\"))\n",
    "    .filterDate(\"2018-01-01\", \"2018-03-31\")\n",
    "    .filterBounds(roi)\n",
    ")\n",
    "\n",
    "s1_Apr18Jan19 = (\n",
    "    ee.ImageCollection(\"COPERNICUS/S1_GRD\")\n",
    "    .filter(ee.Filter.eq(\"instrumentMode\", \"IW\"))\n",
    "    .filter(ee.Filter.listContains(\"transmitterReceiverPolarisation\", \"VH\"))\n",
    "    .filter(ee.Filter.listContains(\"transmitterReceiverPolarisation\", \"VV\"))\n",
    "    .filterDate(\"2018-04-01\", \"2019-01-31\")\n",
    "    .filterBounds(roi)\n",
    ")\n",
    "\n",
    "s1 = s1_JanMar18.merge(s1_Apr18Jan19)\n",
    "\n",
    "#s1mean = s1.mean().select([\"VV\", \"VH\"], [\"VV_mean\", \"VH_mean\"]).clip(roi)\n",
    "#s1stdev = s1.reduce(ee.Reducer.stdDev()).clip(roi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Update: Add Sentinel-1 data (only ascending for Jan to Mar 2018 due to noise) using natural values\n",
    "\n",
    "def natural_s1(img):\n",
    "    #= s1.map(function(img)\n",
    "    return ee.Image(10.0).pow(img.divide(10.0)).copyProperties(img, ['system:time_start'])\n",
    "\n",
    "natural_s1 = s1.map(natural_s1)\n",
    "s1mean_n = natural_s1.mean().select(['VV','VH'],['VV_mean','VH_mean'])\n",
    "s1stdev_n = natural_s1.reduce(ee.Reducer.stdDev())\n",
    "\n",
    "# convert back to dB\n",
    "s1mean_dB = s1mean_n.log10().multiply(10.0)\n",
    "\n",
    "print('S1 data', s1.size().getInfo())\n",
    "\n",
    "s1mean = s1mean_dB\n",
    "s1stdev = s1stdev_n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the Normalized Difference Vegetation Index (NDVI)\n",
    "red = s2TNI.select(\"red_mean\")\n",
    "nir = s2TNI.select(\"nir_mean\")\n",
    "ndvi = nir.subtract(red).divide(nir.add(red)).rename(\"NDVI\")\n",
    "\n",
    "# Compute standard deviation (stDev) as texture of the NDVI\n",
    "texture = ndvi.reduceNeighborhood(\n",
    "    reducer=ee.Reducer.stdDev(), kernel=ee.Kernel.square(5)\n",
    ")\n",
    "\n",
    "# Compute elevation, slope and aspect\n",
    "SRTM = ee.Image(\"USGS/SRTMGL1_003\").clip(roi)\n",
    "slope = ee.Terrain.slope(SRTM)\n",
    "\n",
    "# Add them to bands\n",
    "s2final = (\n",
    "    ee.Image(s2TNI)\n",
    "    .addBands(ndvi)\n",
    "    .addBands(texture)\n",
    "    .addBands(slope)\n",
    "    .addBands(s1mean)\n",
    "    .addBands(s1stdev)\n",
    "    .float()\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bands = [\n",
    "    \"blue_mean\",\n",
    "    \"green_mean\",\n",
    "    \"red_mean\",\n",
    "    \"red1_mean\",\n",
    "    \"red2_mean\",\n",
    "    \"red3_mean\",\n",
    "    \"nir_mean\",\n",
    "    \"red4_mean\",\n",
    "    \"swir1_mean\",\n",
    "    \"swir2_mean\",\n",
    "    \"NDVI\",\n",
    "    \"NDVI_stdDev\",\n",
    "    \"slope\",\n",
    "    \"VH_mean\",\n",
    "    \"VV_mean\",\n",
    "    \"VH_stdDev\",\n",
    "    \"VV_stdDev\",\n",
    "]\n",
    "s2finalWbands = s2final.select(bands)\n",
    "\n",
    "# Scale the image for classification\n",
    "s2classification = s2finalWbands.reproject(\n",
    "    ee.Projection(\"EPSG:4326\").atScale(scale)\n",
    ").reduceResolution(reducer=ee.Reducer.mean(), maxPixels=65536)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "randomSeed = 0\n",
    "n = randomSeed\n",
    "split = 0.5\n",
    "\n",
    "# Get the values for all pixels in each polygon in the training.\n",
    "palmSample = s2classification.sampleRegions(\n",
    "    collection=palm, properties=[\"class\"], scale=scale\n",
    ")\n",
    "rubberSample = s2classification.sampleRegions(\n",
    "    collection=rubber, properties=[\"class\"], scale=scale\n",
    ")\n",
    "otherSample = s2classification.sampleRegions(\n",
    "    collection=other, properties=[\"class\"], scale=scale\n",
    ")\n",
    "shrubSample = s2classification.sampleRegions(\n",
    "    collection=shrub, properties=[\"class\"], scale=scale\n",
    ")\n",
    "bareSample = s2classification.sampleRegions(\n",
    "    collection=bare, properties=[\"class\"], scale=scale\n",
    ")\n",
    "waterSample = s2classification.sampleRegions(\n",
    "    collection=water, properties=[\"class\"], scale=scale\n",
    ")\n",
    "\n",
    "randomPalm = palmSample.randomColumn(\"random\", n)\n",
    "randomRubber = rubberSample.randomColumn(\"random\", n)\n",
    "randomOther = otherSample.randomColumn(\"random\", n)\n",
    "randomShrub = shrubSample.randomColumn(\"random\", n)\n",
    "randomBare = bareSample.randomColumn(\"random\", n)\n",
    "randomWater = waterSample.randomColumn(\"random\", n)\n",
    "\n",
    "trainingSample = (\n",
    "    randomPalm.filter(ee.Filter.lt(\"random\", split))\n",
    "    .merge(randomRubber.filter(ee.Filter.lt(\"random\", split)))\n",
    "    .merge(randomOther.filter(ee.Filter.lt(\"random\", split)))\n",
    "    .merge(randomShrub.filter(ee.Filter.lt(\"random\", split)))\n",
    "    .merge(randomBare.filter(ee.Filter.lt(\"random\", split)))\n",
    "    .merge(randomWater.filter(ee.Filter.lt(\"random\", split)))\n",
    ")\n",
    "\n",
    "testingSample = (\n",
    "    randomPalm.filter(ee.Filter.gte(\"random\", split))\n",
    "    .merge(randomRubber.filter(ee.Filter.gte(\"random\", split)))\n",
    "    .merge(randomOther.filter(ee.Filter.gte(\"random\", split)))\n",
    "    .merge(randomShrub.filter(ee.Filter.gte(\"random\", split)))\n",
    "    .merge(randomBare.filter(ee.Filter.gte(\"random\", split)))\n",
    "    .merge(randomWater.filter(ee.Filter.gte(\"random\", split)))\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"palm\", randomPalm.size().getInfo())\n",
    "print(\"rubber\", randomRubber.size().getInfo())\n",
    "print(\"other\", randomOther.size().getInfo())\n",
    "print(\"shrub\", randomShrub.size().getInfo())\n",
    "print(\"bare\", randomBare.size().getInfo())\n",
    "print(\"water\", randomWater.size().getInfo())\n",
    "# print('training', trainingSample.size().getInfo())\n",
    "# print('testing', testingSample.size().getInfo())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "classifier = ee.Classifier.randomForest(numberOfTrees=100, variablesPerSplit=4).train(\n",
    "    trainingSample, \"class\", bands\n",
    ")\n",
    "\n",
    "classified = s2classification.classify(classifier)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get a confusion matrix representing resubstitution accuracy\n",
    "# trainAccuracy = classifier.confusionMatrix()\n",
    "# print('Resubstitution error matrix: ', trainAccuracy.getInfo())\n",
    "# print('Training overall accuracy: ', trainAccuracy.accuracy().getInfo())\n",
    "\n",
    "# Classify the validation data\n",
    "validated = testingSample.classify(classifier)\n",
    "\n",
    "# Get a confusion matrix representing expected accuracy\n",
    "testAccuracy = validated.errorMatrix(\"class\", \"classification\",[1,2,3,4,5,6])\n",
    "errorMatrix = testAccuracy.array()\n",
    "\n",
    "overallAccuracy = testAccuracy.accuracy()\n",
    "producersAccuracy = testAccuracy.producersAccuracy()\n",
    "consumersAccuracy = testAccuracy.consumersAccuracy()\n",
    "\n",
    "feature = ee.Feature(None)\n",
    "feature = feature.set(\"matrix\", errorMatrix)\n",
    "feature = feature.set(\"results\", overallAccuracy)\n",
    "feature = feature.set(\"producer\", producersAccuracy)\n",
    "feature = feature.set(\"consumer\", consumersAccuracy)\n",
    "accuracy_results = ee.FeatureCollection(feature)\n",
    "\n",
    "# Export as csv\n",
    "export1 = ee.batch.Export.table.toDrive(\n",
    "    collection=accuracy_results,\n",
    "    description=\"export_area1_accuracy\",\n",
    "    fileNamePrefix=\"area1_accuracy\",\n",
    ")\n",
    "export1.start()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# all reference data\n",
    "reference = testingSample.merge(trainingSample)\n",
    "\n",
    "# Classify the reference data\n",
    "validated_r = reference.classify(classifier)\n",
    "\n",
    "# Get a confusion matrix representing expected accuracy\n",
    "testAccuracy = validated_r.errorMatrix(\"class\", \"classification\",[1,2,3,4,5,6])\n",
    "errorMatrix = testAccuracy.array()\n",
    "\n",
    "overallAccuracy = testAccuracy.accuracy()\n",
    "producersAccuracy = testAccuracy.producersAccuracy()\n",
    "consumersAccuracy = testAccuracy.consumersAccuracy()\n",
    "\n",
    "feature = ee.Feature(None)\n",
    "feature = feature.set(\"matrix\", errorMatrix)\n",
    "feature = feature.set(\"results\", overallAccuracy)\n",
    "feature = feature.set(\"producer\", producersAccuracy)\n",
    "feature = feature.set(\"consumer\", consumersAccuracy)\n",
    "accuracy_results = ee.FeatureCollection(feature)\n",
    "\n",
    "# Export as csv\n",
    "export2 = ee.batch.Export.table.toDrive(\n",
    "    collection=accuracy_results,\n",
    "    description=\"export_area1_accuracy_reference\",\n",
    "    fileNamePrefix=\"area1_accuracy_reference\",\n",
    ")\n",
    "export2.start()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run a majority filter\n",
    "# radius (Float, default: 1.5) kernelType (String, default: \"circle\"):units (String, default: \"pixels\")\n",
    "filtered = classified.focal_mode()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export the image to estimate area\n",
    "# to drive\n",
    "export3 = ee.batch.Export.image.toDrive(\n",
    "    image=filtered,\n",
    "    description=\"export_area1_classified_drive\",\n",
    "    fileNamePrefix=\"area1_classified\",\n",
    "    scale=scale,\n",
    "    maxPixels=1e13,\n",
    ")\n",
    "export3.start()\n",
    "\n",
    "# to asset (replace with your username)\n",
    "export4 = ee.batch.Export.image.toAsset(\n",
    "    image=filtered,\n",
    "    description=\"export_area1_classified_asset\",\n",
    "    assetId=\"users/nkeikon/area1_classified\",\n",
    "    scale=20,\n",
    "    maxPixels=1e13,\n",
    ")\n",
    "# export4.start()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Produce final map by manually classifying one young oil palm plantation\n",
    "area1_finalmap = filtered.where(\n",
    "    (classified.clip(op).eq(2).Or(classified.clip(op).eq(3))), 1\n",
    ")\n",
    "\n",
    "export5 = ee.batch.Export.image.toDrive(\n",
    "    image=area1_finalmap,\n",
    "    description=\"export_area1_map_drive\",\n",
    "    fileNamePrefix=\"area1_map\",\n",
    "    scale=scale,\n",
    "    maxPixels=1e13,\n",
    ")\n",
    "export5.start()\n",
    "\n",
    "# to asset (replace with your username)\n",
    "export6 = ee.batch.Export.image.toAsset(\n",
    "    image=area1_finalmap,\n",
    "    description=\"export_area1_map_asset\",\n",
    "    assetId=\"users/nkeikon/area1_map\",\n",
    "    scale=20,\n",
    "    maxPixels=1e13,\n",
    ")\n",
    "# export6.start()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# query current tasks\n",
    "tasks = ee.batch.Task.list()\n",
    "print(tasks[0])\n",
    "print(tasks[1])\n",
    "print(tasks[2])\n",
    "print(tasks[3])\n",
    "print(tasks[4])\n",
    "print(tasks[5])"
   ]
  },
  {
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}