Generate the complement of the Earth's land mass polygons

Question

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

where i need to include seas as well

Answer 1

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]

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

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]

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

Answer 2

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

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

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.

Answer 3

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]

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