Calculate change in water level

Use Planet satellite imagery to calculate change in reservoir water levels

analysis, data

Overview

In this tutorial we are going to use Planet imagery to calculate an approximate percentage change in reservoir water levels.

You will need:

To export your Planet API Key to your environment:

export PLANET_API_KEY=a3a64774d30c4749826b6be445489d3b # (not a real key)

You can find GDAL installation instructions here. For installation of scipy & numpy, see this page.

To install the Python libraries, do:

pip install retrying
pip install requests
pip install bokeh
pip install simplecv

# optionally:
pip install planet

Once installation is complete, you're ready to start with Step 1 below.

A note about Planet Imagery Product Compatability: this tutorial is compatable with all ItemTypes & Asset Types.

Download images at two time points

First, let's download images of the same Reservoir in California, 2 weeks apart.

You can use the Data API to search for, activate & download these images, or optionally you can use the planet CLI.

To use the CLI, do:

$ planet data download --item-type REOrthoTile --asset-type visual --string-in id 20160707_195146_1057917_RapidEye-1
$ planet data download --item-type REOrthoTile --asset-type visual --string-in id 20160722_194930_1057917_RapidEye-2

You now have the two Planet 'visual' GeoTIFF format images in your current directory.

N.B.: As this is just an example, we are using Planet's 'visual' asset type. If we wanted a more accurate measurement we would use the higher bit-depth 'analytic' product.

Select common window to compare

Our scenes don't overlap perfectly, and for the calculation to be accurate we'd prefer they did. The GDAL Warp command enables us to do this crop.

With QGIS we can find the overlapping rectangle between the two scenes. Move your mouse to where you estimate the corners might be, and take note of the numbers from the 'coordinates' box on the bottom menu bar.

We then run the following bash commands:

gdalwarp -te 547500 4511500 556702 4527000 20160707_195146_1057917_RapidEye-1_visual.tif 20160707.tif
gdalwarp -te 547500 4511500 556702 4527000 20160722_194930_1057917_RapidEye-2_visual.tif 20160722.tif

Find the Water

Finding the water is easy, it's blue. We can use HUE. The hue distance transform turns all the water white. We also shrink the images, as we don't need the extra resolution for this simple analysis.

import SimpleCV as scv
a = scv.Image("20160707.tif")
b = scv.Image("20160722.tif")

innerA =  a.hueDistance(175)
innerB =  b.hueDistance(171)
innerA.scale(0.4).rotateLeft().save("hue_a.png")
innerB.scale(0.4).rotateLeft().save("hue_b.png")

Threshold the Hue Transform

Next we apply a threshold to the image based upon hue, to create a mask. The hue distance transform gives us the following images, which are base masks:

innerA =  a.hueDistance(175).threshold(200)
innerB =  b.hueDistance(171).threshold(200)
innerA.scale(0.4).rotateLeft().save("hue_a_threshold.png")
innerB.scale(0.4).rotateLeft().save("hue_b_threshold.png")

Clean up the Masks

We then use morphological operations to clean up these masks. We will call this our "inner mask." We are pretty sure now that areas labeled in white are water.

innerA =  a.hueDistance(175).threshold(200).erode(2).dilate(3).erode(2)
innerB =  b.hueDistance(171).threshold(200).erode(2).dilate(3).erode(2)
innerA.scale(0.4).rotateLeft().save("hue_a_morph.png")
innerB.scale(0.4).rotateLeft().save("hue_b_morph.png")

Find regions we're not sure about

Now we want to find regions we are less sure about. We will try to label these regions as gray.

To do this we will grow our mask, then invert it.

# we use /2 to change white into gray
middleA =  a.hueDistance(175).threshold(200).erode(1).dilate(2).invert()/2
middleB =  b.hueDistance(171).threshold(200).erode(1).dilate(2).invert()/2
middleA.scale(0.4).rotateLeft().save("middle_a.png")
middleB.scale(0.4).rotateLeft().save("middle_b.png")

Create the watershed mask

Now we put it all together. Unlike a normal mask, a watershed mask has three parts: White -- definitely our object. Gray -- definitely not our object. * Black -- unsure / boundary.

The watershed algorithm will pick the optimal point between the gray and white boundaries. Because gray pixel values are a higher number than black ones, summing them gives us what we want, as long as the white area is smaller than the black area.

innerA =  a.hueDistance(175).threshold(190).erode(1).dilate(1).erode(1)
innerB =  b.hueDistance(171).threshold(200).erode(1).dilate(1).erode(1)
middleA =  a.hueDistance(175).threshold(190).erode(1).dilate(3).invert()/2
middleB =  b.hueDistance(171).threshold(200).erode(1).dilate(3).invert()/2
a_mask = innerA+middleA
b_mask = innerB+middleB
a_mask.scale(0.4).rotateLeft().save("mask_a.png")
b_mask.scale(0.4).rotateLeft().save("mask_b.png")

Run the watershed algorithm

aBlobs = a.findBlobsFromWatershed(mask=a_mask)
bBlobs = b.findBlobsFromWatershed(mask=b_mask)
aBlobs.draw(color=scv.Color.RED,width=-1,alpha=128)
a = a.applyLayers()
a.scale(0.4).rotateLeft().save("ws_a.png")
bBlobs.draw(color=scv.Color.HOTPINK,width=-1,alpha=128)
b = b.applyLayers()
b.scale(0.4).rotateLeft().save("ws_b.png")

Calculate the water-level change

import numpy as np

# best support is with data in a format that is table-like
aa =  np.sum(aBlobs.area())
bb =  np.sum(bBlobs.area())
print aa
print bb
print bb/float(aa)
print "Percent change {0}%".format(((bb/float(aa))-1)*100)

The output you get will look like this:

1414665.5
1357179.0
0.95936389203
Percent change -4.06361079704%

Finally, we can also plot these data in a bar chart using BokehJS:

from bokeh.charts import Bar
from bokeh.io import output_file, show
import bokeh.plotting as bk

data = {
    'date': ['2016-07-13', '2016-09-10'],
    'pixels': [aa, bb]}
bar2 = Bar(data, values='pixels', label=['date'], title="Reservoir pixels", width=800)

# Saves chart as waterlevel.html
bk.save(bar2)

If you have any trouble with this guide, let us know at devex@planet.com.