Getting Area from GeoGraphics

Question

I want to find the area covered by the reservoir of a dam using a satellite image.

This is what I have so far.

DamPosition = DamData[Entity["Dam", "TehriDam::q2zsw"], "Position"]

Which returns the coordinates. And then, I use some image manipulation to isolate the dam from geoImage.

geoImage = 
 DeleteSmallComponents[
  Dilation[
   Graphics[
    Select[Cases[
      GeoGraphics[
       GeoBoundingBox[GeoPosition[{30.377778`, 78.480556`}]], 
       GeoBackground -> 
        "VectorMinimal"], {Directive[{___, 
         RGBColor[0.6, 0.807843137254902`, 1.], ___}], ___}, 
      Infinity], Not@*FreeQ[Polygon]]], 0]]

Now, all I want to do is find the area of the blue - shaded region.

I know I can do this manually using the getcoordinates tool, but that seems too unreliable and tedious to do. I'm looking for an automated way to do it.

PS : Using GetCoordinates, we can use the following code

    damContour = {{36.14216500804076`, -98.66601000992658`},
   {36.15131161335255`, -98.6549477815357`},
   {36.1501032100147`, -98.62214348441108`},
   {36.15000037203102`, -98.60908533742007`},
   {36.13193732747872`, -98.58352140011719`},
   {36.10128112030604`, -98.57004960135653`},
   {36.08009879008038`, -98.6043835918035`},
   {36.1172565230945`, -98.61315564049649`},
   {36.12831708030723`, -98.62663775701455`},
   {36.13202116684048`, -98.64507269213401`},
   {36.13870030579935`, -98.6454368865287`},
   {36.13708633726377`, -98.6568745786338`}};
GeoGraphics[{White, Thick, Line[GeoPosition[damContour]]}, 
 GeoRange -> {{36.05, 36.17}, {-98.7, -98.52}}, 
 GeoRangePadding -> Scaled[0.1], GeoBackground -> "Satellite"];
GeoArea[Polygon[GeoPosition[damContour]]]

Answer 1

First I get just the image of the area with a 10 km scale bar:

DamPosition = DamData[Entity["Dam", "TehriDam::q2zsw"], "Position"];

l = 10;
geoImage = 
 GeoGraphics[GeoBoundingBox[DamPosition], 
  GeoBackground -> "StreetMapNoLabels", 
  GeoScaleBar -> Placed[Quantity[{0, l}, "Kilometers"], {Right, Top}]]

Then I find the areas that have the color of water, and Binarize:

waterColor = RGBColor[0.6, 0.807843137254902, 1.];
lake = ColorDetect[geoImage, ColorsNear[waterColor, 0]] // Binarize

I get just the part that looks like the reservoir by taking the largest component of lake:

res = SelectComponents[lake, "Count", -1]

Now I get the scale bar by finding the black parts of geoImage:

bar = Binarize@ColorDetect[geoImage, RGBColor[0, 0, 0]]

To get the scale bar without the labels, I select the longest length component of bar:

bar = SelectComponents[bar, "Length", -1]

Then I get the pixel length of the scale bar by grabbing the largest side length of its BoundingBox, and calculate the scale by dividing l by the pixel length of the scale bar

v = First@Values@ComponentMeasurements[bar, "BoundingBox"];
len = Max[Abs[First[v] - Last[v]]];
scale = Quantity[l, "km"]/len
(*Quantity[0.04329, "Kilometers"] *)

The area of each pixel is just scale^2:

pixelArea=scale^2;

The total area of res is now just the total of the image (since it's binarized) times `pixelArea

Total@Flatten[ImageData[res]]*pixelArea

(*Quantity[20.5019, ("Kilometers")^2]*)

So the area of the reservoir is $\approx 20.5 \mathrm{km}^2$

Answer 2

Alternate answer that I have received

DamPosition = DamData[Entity["Dam", "TehriDam::q2zsw"], "Position"];

l = 10;
geoImage = 
 GeoGraphics[GeoBoundingBox[DamPosition], 
  GeoBackground -> "StreetMapNoLabels", 
  GeoScaleBar -> Placed[Quantity[{0, l}, "Kilometers"], {Right, Top}]]
waterColor = RGBColor[0.6, 0.807843137254902, 1.];
geoImage = ColorDetect[geoImage, ColorsNear[waterColor, 0]] // 
Binarize
 ImageResize[geoImage, {128,128}]
 Area[ImageMesh[geoImage]]

This becomes useful when we resize all images to a standard size and then compare the area of different dams.