Message Boards Message Boards

Satellite Image Processing with Python and the Wolfram Language

GROUPS:

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
Answer
15 days ago

Group Abstract Group Abstract