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$
  • $\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 '17 at 19:14
  • $\begingroup$ added a specific example, where the seas are missing in the plot $\endgroup$ – Karolis Aug 13 '17 at 19:36
  • $\begingroup$ Would using GeoBackground -> GeoStyling[LightBlue] suffice? $\endgroup$ – Chip Hurst Aug 14 '17 at 14:22
  • 1
    $\begingroup$ @ChipHurst data and the notebook is available on github.com/kmisiunas/data-sunshine-vs-population $\endgroup$ – Karolis Aug 14 '17 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$ – Chip Hurst Aug 14 '17 at 15:00
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

| improve this answer | |
$\endgroup$
  • $\begingroup$ It's a pity that the generated mesh is huge, it looks fine otherwise. $\endgroup$ – J. M.'s discontentment Aug 15 '17 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.

| improve this answer | |
$\endgroup$
  • $\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 '17 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

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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