Inspired by Planet’s guide on calculating change in water level, I wanted to re-create their visual analysis environment from within the Wolfram Language.
Using the Planet API
The Planet API is easily leveraged from their quickstart python here.
Browsing the Planet Explorer we can find some cool assets, note their names, and store to variables:
idList = {"20160804_180723_0e26", "20160820_180756_0e30"}; (*id of image*)
itemType = "PSScene4Band"; (*Type of satellite/image *)
assetType = "analytic"; (*analytic = cool multi spectrum image, visual = visible spectrum image*)
path = "images/"; (*where you want to output the images*)
The new ExternalEvaluate will allow us to write an easy wrapper to access the APIs functionality:
session=StartExternalSession[<|"System"->"Python","Version"->"2.7.12 :: Anaconda 4.2.0 (x86_64)", "ReturnType"->"String"|>];
(*change working directory and add the api key to the environment*)
ExternalEvaluate[session,"import os; os.chdir('/Users/peter/Downloads')"]
ExternalEvaluate[session,"os.environ['PLANET_API_KEY'] = "~~PersistentValue["PLANETAPIKEY"]]
(*overload an uninitialized global reference*)
ExternalEvaluate[session,
"class dummy:
download = 'downloading images'"]
ExternalEvaluate[session, "args = dummy()"]
(*Load full download script, downloaded from *)
ExternalEvaluate[session,File["/Users/peterro/Downloads/download.py" ]]
convertToJSON[strAr_] := ExportString[strAr,"JSON", "Compact"->True];
inputArgs = <|"idList"->idList,"itemType"-> itemType, "assetType"->assetType |>;
activateArgs = Join[inputArgs,<|"activateOrCheck"->"activate"|>];
checkArgs = Join[inputArgs,<|"activateOrCheck"->"activate"|>];
downloadArgs = Join[inputArgs,<|"path"->path|>];
processActivation = StringTemplate["process_activation(activate, `idList`, `itemType`, `assetType`, `activateOrCheck`)"];
processDownload = StringTemplate["process_download(`path`, `idList`, `itemType`, `assetType`, '')"];
(*activate the assets*)
ExternalEvaluate[session,processActivation[convertToJSON /@ activateArgs]];
obj=SessionSubmit[ScheduledTask[
(*check if the assets are active*)
If[ImportString[Last@StringSplit[ExternalEvaluate[session,processActivation[convertToJSON /@ checkArgs]],"\n"]] != "[True, True]" ,
Module[{},
Print["Assets active, downloading"];
TaskRemove[$CurrentTask];
(*download the assets*)
Print[ExternalEvaluate[session,processDownload[convertToJSON /@ downloadArgs]]];
DeleteObject[session];
Print["Download Complete"];
]
,
Print["Waiting for download: " ~~ DateString[]]
],{10,10}]];
This will asynchronously download the images from the Planet servers when they are prepared for download.
Process the Images
Import:
img820 = ImageResize[First@Import["20160820_180756_0e30_analytic.tif"], Scaled[.4]];
img804 = ImageResize[First@Import["20160804_180723_0e26_analytic.tif"], Scaled[.4]];
These look weird because we’re viewing the stacked {R,G,B,IR} channels. For the purposes of calculating water level change, I wanted to use the IR channel:
infra820 = ImageAdjust[ColorCombine[ColorSeparate[img820][[{4}]]]];
infra804 = ImageAdjust[ColorCombine[ColorSeparate[img804][[{4}]]]];
Next, we’ll binarize the image using thresholds derived from the image histogram:
binarized820 = MorphologicalBinarize[infra820, {0.172, 0.896}];
binarized804 = MorphologicalBinarize[infra804, {0.172, 0.896}];
These images are hard to work with because they are rotated and not cropped so we can use some cool image rotation to correct that
First identify the lines:
lines820=ImageLines[GradientFilter[binarized820,1], MaxFeatures->5];
lines804=ImageLines[GradientFilter[binarized804,1], MaxFeatures->5];
Show[#,ImageSize->Large]&/@ {HighlightImage[binarized820,{Orange,Line/@lines}],HighlightImage[binarized804,{Orange,Line/@lines2}]}
Rotate and crop:
rotateCrop820 = ColorNegate@ImagePad[ImageCrop[ImageRotate[binarized820,Median[Select[-ArcTan[(#[[2,2]]-#[[1,2]])/(#[[2,1]]-#[[1,1]])]&/@
ImageLines[GradientFilter[binarized820,1]],#<20/360 2 \[Pi]&&#>-20/360 2 \[Pi]&]]]],-20]
Next we’ll use the overlapping key points to create a comparable area between the two images:
crops= { rotateCrop804,rotateCrop820};
matches = ImageCorrespondingPoints[Sequence@@crops, KeypointStrength->.01] ;
{matches804,matches820} = matches;
Column@MapThread[Show[#1,Graphics[{Blue,MapIndexed[Inset[#2[[1]],#1]& ,#2]}], ImageSize->Large]&,{crops,matches}]
Use the bounding regions to create the masks:
Tie it all together and measure the change in water level
waterPixels820 = Total@Flatten@ImageData@ImageTrim[rotateCrop820,Round /@ BoundingRegion[matches820] /. Cuboid->List]
waterPixels804 = Total@Flatten@ImageData@ImageTrim[rotateCrop804,Round /@ BoundingRegion[matches804] /. Cuboid->List]
We get 575304 and 597415 meaning the water level decreased by 3.7% over these days.
