Log inSuggest an edit

Calculate change in water level

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

You will need:

  • Python 2.7 or 3
  • GDAL (Installation Instructions)
  • retrying - pip install retrying
  • requests - pip install requests
  • SimpleCV - pip install https://github.com/ingenuitas/SimpleCV (install using this fork to avoid an error in Step 8). Also requires NumPY and OpenCV.
  • BokehJS for plotting graphs - pip install bokeh
  • Your Planet API Key in your environment:
export PLANET_API_KEY=a3a64774d30c4749826b6be445489d3b # (not a real key)

Planet Imagery Product Compatability: all ItemTypes & Asset Types.

The approach:

  1. Download images at two time points
  2. Select common window to compare
  3. Find the water
  4. Threshold the Hue Transform
  5. Clean up the Masks
  6. Find regions we're not sure about
  7. Create the watershed mask
  8. Run the watershed algorithm
  9. Calculate the change

Step 1. Download images at two time points:

First, let's download images of the same Reservoir in California, 2 weeks apart. Get the download script from the API Quickstart Step 3, download.py, and create a text file that contains the following two IDs e.g. items.txt:

  • 20160707_195146_1057917_RapidEye-1
  • 20160722_194930_1057917_RapidEye-2

Activate these two items:

python download.py --idlist items.txt --activate REOrthoTile visual

If they're not already 'active', check every minute or two:

python download.py --idlist items.txt --check REOrthoTile visual

Once they become active, you can download them:

python download.py --idlist items.txt --download images/ REOrthoTile visual

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

NB 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.

Step 2: 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

Step 3: 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")

Step 4: 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")

Step 5: 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")

Step 6: 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")

Step 7: 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")

Step 8: 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")

Step 9: 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 support@planet.com.