June 14 2021

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)

Questions or comments about this guide? Join the conversation at Planet Community.

June 14 2021

Normalized Difference Vegetation Index (NDVI), developed by a NASA scientist named Compton Tucker in 1977, is commonly used to assess whether an area contains live green vegetation or not. It can show the difference between water and plants, bare soil and grass, whether plants are under stress, and what lifecycle stage a crop is in.

It compares how much more near-infrared light is reflected by chlorophyll vs visible red light:

In this guide we will perform a basic NDVI calculation in python, using Planet's 4-Band imagery.

You will need,

To export your Planet API Key to your environment:

export PLANET_API_KEY=a3a64774d30c4749826b6be445489d3b # (not a real key)

You can find Numpy installation instructions here.

To install the Python libraries, do:

pip install rasterio
pip install matplotlib
pip install retrying
pip install requests

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 only analytic and basic_analytic asset types of the following ItemTypes:

  • PSOrthoTile
  • REOrthoTile
  • PSScene4Band

Download a 4-Band image

First, download a 4-band Planetscope image of agricultural land in Toledo, Spain (id: 20161218_101700_0e0d). You can do this using the Planet API, or with Planet Explorer, by filtering for '4 Band PlanetScope scene' (PSScene4Band) or 'Planetscope ortho tile' (PSOrthoTile), and downloading an 'analytic' asset.

Using the Planet CLI:

planet data download --item-type PSScene4Band --asset-type analytic,analytic_xml --string-in id 20161218_101700_0e0d

You now have a file called '20161228_101647_0e26_3B_AnalyticMS.tif' in your directory.


Thumbnail preview of the 4-band scene

Extract the Visible Red and NIR bands

In a python script, open the image and extract just the Red and Near-infrared bands, numbers 3 and 4.

import rasterio
import numpy

image_file = "20161228_101647_0e26_3B_AnalyticMS.tif"

# Load red and NIR bands - note all PlanetScope 4-band images have band order BGRN
with rasterio.open(image_file) as src:
    band_red = src.read(3)

with rasterio.open(image_file) as src:
    band_nir = src.read(4)

Normalize to Top of Atmosphere Reflectance

Converting the pixel values to TOA Reflectance makes the analysis more accurate, and comparable with other scenes. Load the TOA Reflectance coefficients from the metadata XML asset.

from xml.dom import minidom

xmldoc = minidom.parse("20161218_101700_0e0d_3B_AnalyticMS_metadata.xml")
nodes = xmldoc.getElementsByTagName("ps:bandSpecificMetadata")

# XML parser refers to bands by numbers 1-4
coeffs = {}
for node in nodes:
    bn = node.getElementsByTagName("ps:bandNumber")[0].firstChild.data
    if bn in ['1', '2', '3', '4']:
        i = int(bn)
        value = node.getElementsByTagName("ps:reflectanceCoefficient")[0].firstChild.data
        coeffs[i] = float(value)

Multiply the band values by the TOA Reflectance coefficients.

# Multiply by corresponding coefficients
band_red = band_red * coeffs[3]
band_nir = band_nir * coeffs[4]

Perform the NDVI calculation

Next we perform the NDVI calculation through subtraction and division of the pixel values.

# Allow division by zero
numpy.seterr(divide='ignore', invalid='ignore')

# Calculate NDVI
ndvi = (band_nir.astype(float) - band_red.astype(float)) / (band_nir + band_red)

Save the NDVI image

Finally we output these new pixel values to a new image file, making sure we mirror the GeoTIFF spatial metadata:

# Set spatial characteristics of the output object to mirror the input
kwargs = src.meta
kwargs.update(
    dtype=rasterio.float32,
    count = 1)

# Create the file
with rasterio.open('ndvi.tif', 'w', **kwargs) as dst:
        dst.write_band(1, ndvi.astype(rasterio.float32))

Apply a color map

Applying a color map can help to visually distinguish vegetation. For more color map options check out the docs for matplotlib.

import matplotlib.pyplot as plt
plt.imsave("ndvi_cmap.png", ndvi, cmap=plt.cm.summer)

Questions or comments about this guide? Join the conversation at Planet Community.

June 14 2021

THIS FUNCTIONALITY HAS BEEN DEPRECATED

the following guide is for informational purposes only, and does not reflect the current state of Planet's platform.

Overview

This guide explains how to download a subarea of a Planet image, saving time, bandwidth and storage. We will use an exciting new feature of the GDAL library, VSI Curl, which queries a part of a GeoTIFF file using an HTTP Range request.

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.

To install the Python libraries, do:

pip install retrying
pip install requests

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.

Activate a scene containing your Subarea

First, we will activate a scene of a stretch of agricultural land in Yuma, AZ, with the Planet API client.

planet data download --activate-only --item-type PSScene3Band --asset-type visual --string-in id 20161109_173041_0e0e --quiet

Activation will take a few minutes. Run the command above again until you see a URL in the 'location' element of the response.

A full-size Planet Scene of Agricultural land in Yuma, AZ

Create a GeoJSON file of your Subarea

Define the geometry of the subarea you'd like to download in GeoJSON format, by drawing it at geojson.io:

Save the GeoJSON to a file e.g. subarea.geojson,

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type":"Polygon",
        "coordinates":[
          [
             [
                -114.176423930378903,
                32.691933554346988
             ],
             [
                -114.176423930378903,
                32.685249521402788
             ],
             [
                -114.168019741897595,
                32.685249521402788
             ],
             [
                -114.168019741897595,
                32.691933554346988
             ],
             [
                -114.176423930378903,
                32.691933554346988
             ]
          ]
        ]
      }
    }
  ]
}

Download the Subarea

Next let's request a download URL for this image from the Planet Data API. Note, you must do this every time you wish to download your subarea, as these download URLs expire after 1 hour.

from osgeo import gdal
from requests.auth import HTTPBasicAuth
import os
import requests

item_id = "20161109_173041_0e0e"
item_type = "PSScene3Band"
asset_type = "visual"
item_url = 'https://api.planet.com/data/v1/item-types/{}/items/{}/assets'.format(item_type, item_id)

# Request a new download URL
result = requests.get(item_url, auth=HTTPBasicAuth(os.environ['PLANET_API_KEY'], ''))
download_url = result.json()[asset_type]['location']

Prefix the download URL with '/vsicurl/' to use the VSI Curl driver. Then call GDAL Warp, passing it an output file name, this URL, our subarea GeoJSON, a projection (in this case, the one used by GeoJSON.io), and 'cropToCutline = True'.

vsicurl_url = '/vsicurl/' + download_url
output_file = item_id + '_subarea.tif'

# GDAL Warp crops the image by our AOI, and saves it
gdal.Warp(output_file, vsicurl_url, dstSRS = 'EPSG:4326', cutlineDSName = 'subarea.geojson', cropToCutline = True)

Within a couple of seconds, you'll have an image of just the subarea in your directory.


Questions or comments about this guide? Join the conversation at Planet Community.