Message Boards Message Boards


Computing change in water level via processing of satellite images

Posted 1 year ago
1 Reply
13 Total Likes

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*)
"class dummy:
    download = 'downloading images'"]
ExternalEvaluate[session, "args = dummy()"]
(*Load full download script, downloaded from *)
ExternalEvaluate[session,File["/Users/peterro/Downloads/" ]]

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]];

 (*check if the assets are active*)
 If[ImportString[Last@StringSplit[ExternalEvaluate[session,processActivation[convertToJSON /@ checkArgs]],"\n"]] !=  "[True, True]" , 
         Print["Assets active, downloading"];
         (*download the assets*)
         Print[ExternalEvaluate[session,processDownload[convertToJSON /@ downloadArgs]]];
         Print["Download Complete"];
Print["Waiting for download: " ~~ DateString[]]

This will asynchronously download the images from the Planet servers when they are prepared for download.

Process the Images


img820 = ImageResize[First@Import["20160820_180756_0e30_analytic.tif"], Scaled[.4]];
img804 = ImageResize[First@Import["20160804_180723_0e26_analytic.tif"], Scaled[.4]];

Initial Import

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}]]]];

IR chan adjusted

Next, we'll binarize the image using thresholds derived from the image histogram:

img 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}]}

bin with lines

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]

rotated and cropped 820

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:

bounding region

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.

enter image description here - Congratulations! This post is now a Staff Pick as distinguished by a badge on your profile! Thank you, keep it coming!

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract