Overview¶
The purpose of this tutorial is to use Planet imagery to calculate an approximate percentage change in reservoir water levels.
How To Get Started¶
Before we start our analysis we need to ensure we have the required packages and Planet satellite images.
Requirements¶
Software & Python Libraries¶
You will need:
- Python 3
- GDAL
- scipy & numpy
- skimage
- requests
- OpenCV
- imutils
- The Planet CLI (optional)
- Your Planet API Key
How to Export Planet API Key and Install Packages¶
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:
$ python -m pip install requests
$ python -m pip install opencv-python
$ python -m pip install scikit-image
$ python -m pip install imutils
# optionally:
$ python -m 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.
How to Download Images with 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.
Use QGIS to 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.
Crop Images With GDAL Warp Command¶
We then run the following bash commands:
gdalwarp -te 547500 4511500 556702 4527000 1057917_2016-07-07_RE1_3A_Visual.tif 20160707.tif
gdalwarp -te 547500 4511500 556702 4527000 1057917_2016-07-22_RE2_3A_Visual.tif 20160722.tif
Find the Water¶
In order to find the water we want to extract the blue hues within the image. Using OpenCV, we convert the BGR colorspace to HSV and create a threshold to extract the blue water. As you can see below, the mask we created differentiates the water from the land.
Create a Mask Using Threshold¶
import cv2 as cv
a = cv.imread('20160707.tif')
b = cv.imread('20160722.tif')
hsv_a = cv.cvtColor(a, cv.COLOR_BGR2HSV)
hsv_b = cv.cvtColor(b, cv.COLOR_BGR2HSV)
low = np.array([55, 0, 0])
high = np.array([118, 255, 255])
inner_a = cv.inRange(hsv_a, low, high)
inner_b = cv.inRange(hsv_b, low, high
cv.imwrite('inner_a.png', inner_a)
cv.imwrite('inner_b.png', inner_b)
Mask A¶
Mask B¶
Cleaning up the Mask Using Erosion and Dilation¶
We currently have some noise in our mask that we can eliminate using two morphological operations, erosion and dilation. The purpose of these operations to ensure we have a clear separation between the background and foreground, or land and water.
kernel = np.ones((5,5),np.uint8)
erosion_a = cv.erode(inner_a,kernel,iterations = 2)
erosion_b = cv.erode(inner_b,kernel,iterations = 2)
innerA = cv.dilate(erosion_a,kernel,iterations = 1)
innerB = cv.dilate(erosion_b,kernel,iterations = 1)
cv.imwrite('innerA.png', innerA)
cv.imwrite('innerB.png', innerB)
Mask A¶
Mask B¶
Watershed Algorithm¶
In order to calculate the changes in water levels we need to know the area of the water in the image. We will use a segmentation technique so we only focus on the water and ignore everything else. The watershed algorithm returns a numpy array of labels with unique values corresponding to the pixel value.
from scipy import ndimage
from skimage.feature import peak_local_max
from skimage.segmentation import watershed
import numpy as np
eucl_a = ndimage.distance_transform_edt(innerA)
eucl_b = ndimage.distance_transform_edt(innerB)
localMaxA = peak_local_max(eucl_a, indices=False, labels=innerA)
localMaxB = peak_local_max(eucl_b, indices=False, labels=innerB)
markers_a = ndimage.label(localMaxA, structure=np.ones((3, 3)))[0]
markers_b = ndimage.label(localMaxB, structure=np.ones((3, 3)))[0]
labels_a = watershed(-eucl_a, markers_a, mask=innerA)
labels_b = watershed(-eucl_b, markers_b, mask=innerB)
Calculating Area Using Contours and Overlaying Mask¶
Once we have our labels, we loop over each unique label and look for the values that correspond with the foreground, or water. After applying those values to our mask, we then grab the contour of each object in our mask in order to perform our calculations on it. By adding the area of each individual contour to the total area, we are able to approximate the area of the water in our image.
Note: I added lines to fill the contours and save the image for visualization purposes but this isn't necessary for our calculation.
Python Function to Find Area¶
def get_area(labels, inner_mask, img):
area = 0
for label in np.unique(labels):
if label== 0:
continue
mask = np.zeros(inner_mask.shape, dtype="uint8")
mask[labels == label] = 255
contours = cv.findContours(mask, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE)
contours = imutils.grab_contours(contours)
cv.fillPoly(img, pts=contours, color=(0,255,255))
a = cv.contourArea(contours[0])
area+=a
cv.imwrite('mask.png', img)
return area
area_a = get_area(labels_a, innerA, a)
area_b = get_area(labels_b, innerB, b)
Image A¶
Image B¶
Apply Function to get Area and Water Change¶
After we apply our function to both planet images, we now know the approximate change in water reservoir water levels between the times both these images were captured.
area_a = mask_img(labels_a, innerA, a)
area_b = mask_img(labels_b, innerB, b)
water_level_diff = area_b/float(area_a)
Output¶
Area A = 1164765.0
Area B = 1120738.5
Water Level Difference = 0.9622013882628685
Percent change = -3.7798611737131504%
Data Visualization Section¶
Our last step is to plot a bar chart to represent the difference in water levels.
from bokeh.io import output_file, show
from bokeh.plotting import figure
import bokeh.plotting as bk
dates = ['2016-07-13', '2016-09-10']
pixels = [area_a, area_b]
plt = figure(x_range=dates, plot_height=275, title="Reservoir Pixels",
toolbar_location=None, tools="")
plt.vbar(x=dates, top=pixels, width=0.3, fill_color="#cfe31e")
plt.xgrid.grid_line_color = None
plt.y_range.start = 0
plt.xaxis.axis_label = "Date"
plt.yaxis.axis_label = "Sum(Pixels)"
#saves as a html file
bk.save(plt)
show(plt)
Water Changes¶
Questions or comments about this guide? Join the conversation at Planet Community.
Rate this guide: