Message Boards Message Boards

Computing change in water level via processing of satellite images

Posted 6 years ago

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

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

binarized

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

labels

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.

POSTED BY: Peter Roberge
5 Replies
Posted 6 months ago

Hello peter, have nice day. The method of receiving the satellite image was interesting. Is it possible to get an image from the USGS site?

POSTED BY: Erfan Abdi

I am currently working towards this. Will update you as soon as the code is done

POSTED BY: Navvye Anand
Posted 6 months ago

Thank you very much Navvye.

POSTED BY: Erfan Abdi
Posted 2 years ago

H Peter,

Great work. Congrats!

Is it possible to do the same compuations for a lake e.g. Urmia and Caspian ?

Thank you so much.

POSTED BY: Alex Teymouri

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!

POSTED BY: Moderation Team
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract