4
$\begingroup$

I would like to get polygon covering of all water area of earth. I though of doing by inverting

Polygon[EntityClass["Country", "LandMasses"]]

How do i invert this polygon, such that it covers all the oceans and seas?

Specific use case example:

GeoListDensityPlot[data_] := Module[
  {colorRange, colorRegion, color},
  colorRange = MinMax[data];
  color = ColorData["TemperatureMap"];
  colorRegion[pos_ -> value_] := { 
    color[(value - colorRange[[1]])/(colorRange[[2]] - 
        colorRange[[1]])], area[pos]};
  GeoGraphics[{
    colorRegion /@ Normal[data],
    EdgeForm[Black], GeoStyling["ReliefMap"], 
    Polygon[EntityClass["Ocean", "Oceans"]]
    } , GeoBackground -> GeoStyling["Coastlines"] ]
  ]

enter image description here

where i need to include seas as well

$\endgroup$
7
  • $\begingroup$ How are you intending to use the resulting polygon? Do you just need to for mapping or for some sort of computation? $\endgroup$
    – Edmund
    Aug 13, 2017 at 19:14
  • $\begingroup$ added a specific example, where the seas are missing in the plot $\endgroup$
    – Karolis
    Aug 13, 2017 at 19:36
  • $\begingroup$ Would using GeoBackground -> GeoStyling[LightBlue] suffice? $\endgroup$
    – Greg Hurst
    Aug 14, 2017 at 14:22
  • 1
    $\begingroup$ @ChipHurst data and the notebook is available on github.com/kmisiunas/data-sunshine-vs-population $\endgroup$
    – Karolis
    Aug 14, 2017 at 14:41
  • 1
    $\begingroup$ [This comment is unrelated to inverting a Polygon, but my guess is you're wanting to just make a nice display. Correct me if I'm wrong...]. For the purposes of display, I don't see why you can't just set the GeoBackground to a color. Also it seems you have a regular grid of data points. If you're wanting to color the background, would making an Image be a better choice? $\endgroup$
    – Greg Hurst
    Aug 14, 2017 at 15:00

3 Answers 3

4
$\begingroup$

Polygon itself doesn't allow for holes and inverting a polygon will usually result in holes.

Two ways around this are with FilledCurve or BoundaryMeshRegion.

FilledCurve

This requires some manual labor. Here's the area I'd like to invert:

polygon = CountryData["World", "SchematicPolygon"];

GeoGraphics[polygon]

enter image description here

world = Line[GeoPosition[{{-90, -180}, {-90, 180}, {90, 180}, {90, -180}}]];
holes = List /@ Line /@ GeoPosition /@ pp[[1, 1]];
reg = FilledCurve[Prepend[holes, {world}]];

GeoGraphics[{Red, reg}, GeoRange -> "World"]

enter image description here

BoundaryMeshRegion

We can start out with a filled in earth $\left([-180, 180] \times [-90, 90]\right)$ and use RegionDifference to subtract away the polygons. I'm converting everything to mesh regions because RegionDifference usually will respect open/closed-ness of a point set, but ignores this for meshes.

For whatever reason, BoundaryDiscretizeGraphics is slow for large polygons. But Nonetheless

world = BoundaryDiscretizeGraphics[Rectangle[{-180, -90}, {180, 90}]];

holes = BoundaryDiscretizeGraphics /@ EntityValue[EntityClass["GeographicRegion", "Continents"], "Polygon"];

diff = Fold[RegionDifference, world, holes]

enter image description here

This is a pretty large mesh, which might be a reason why it was so slow:

MeshCellCount[diff]
{241629, 241629}

Let's drop that stray horizontal line at the bottom:

First[MaximalBy[ConnectedMeshComponents[diff], Area]]

enter image description here

$\endgroup$
1
  • $\begingroup$ It's a pity that the generated mesh is huge, it looks fine otherwise. $\endgroup$ Aug 15, 2017 at 14:53
6
$\begingroup$

The answer is almost there in the OP. The "Ocean" Entity is used there but only the "Oceans" EntityClass is used. There are more classes of water bodies needed other than "Oceans".

If you are not particular about overlapping polygons you may use EntityClass["Ocean", {"HasPolygon" -> True}].

GeoGraphics[{
  GeoStyling[],
  EdgeForm[Black],
  Opacity[.5, Red],
  EntityValue[EntityClass["Ocean", {"HasPolygon" -> True}], "Polygon"]
  }]

Mathematica graphics

As mentioned you do get some overlaps as the both the North and South polygons of the Atlantic and Pacific Oceans cover the same area as the larger (not split) Atlantic and Pacific Ocean polygons.

The can be avoided by not including the "SevenSeas" EntityClass. Unfortunately, excluding this class also excludes the oceans. Therefore, I included all classes except for the "SevenSeas" class. One quirk of this is that the "NorthwestPassage" is not in any class so it has to be included separately.

GeoGraphics[{
  GeoStyling[],
  EdgeForm[Black],
  Opacity[.5, Red],
  OceanData["NorthwestPassages", "Polygon"],
  EntityValue[
   EntityClass["Ocean",
    {"HasPolygon" -> True,
     "Memberships" -> 
      ContainsAny[{"Oceans", "Seas", "Bays", "Basins", "Straits", "Channels", "Gulfs"}]}], 
   "Polygon"]
  }]

Mathematica graphics

There is still some overlap in the Mediterranean Sea where it is overlapped by the polygons for is Western and Eastern basins. It should not be too difficult for you to remove these to entities before graphing if this overlap is of concern.

Hope this helps.

$\endgroup$
1
  • $\begingroup$ This is a great answer. Thank you. I am still curious if anyone can produce an inverse from landmass polygon $\endgroup$
    – Karolis
    Aug 13, 2017 at 21:27
3
$\begingroup$

Just for fun. Here's a way to visualize your data and remove the small islands using image processing.

colorRange = MinMax[sunshine];
color = ColorData["Rainbow"];
cf[value_] := color[(value - colorRange[[1]])/(colorRange[[2]] - colorRange[[1]])]

im = Image[Reverse[
  Partition[If[#2, cf[#1], Black] & @@@ Values[Merge[{sunshine, landQ}, Join]], 360]]];

mask = DeleteSmallComponents[Image[Reverse[
  Partition[If[#2, White, Black] & @@@ Values[Merge[{sunshine, landQ}, Join]], 360]]], 60];

heatmap = ColorReplace[ImageMultiply[im, mask], Black -> LightGray]

enter image description here

And here's the image with geographic borders to show better contrast (though the islands show here):

GeoGraphics[{EdgeForm[Black], FaceForm[], 
  Polygon[GeoVariant[Entity["GeographicRegion", "World"], "SimplifiedArea"]]}, 
  GeoBackground -> GeoStyling[{"Image", heatmap}], GeoCenter -> {0, 0}]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.