How can we find the longitude that cuts all the land on Earth in half?

Question

Mathematica has already a specific function, RegionCentroid, that can be utilized to calculate the geographic center of any region. As it has been discussed in here and maybe some other questions.

There is also a nice post by Christopher Wolfram which discusses some issues on calculating the centroid of a large geographical region. I have used the ideas in there to calculate the geographical center of all lands on Earth. The result perfectly matched with this post on GIS SE. The main function was this:

geoRegion3D[meshI : (_MeshRegion | _BoundaryMeshRegion)] := 
 Block[{mesh = Quiet@TriangulateMesh@DiscretizeRegion@meshI}, 
  If[Head[mesh] =!= MeshRegion, Missing[], MeshRegion[GeoPositionXYZ[
    GeoPosition[Reverse/@ MeshCoordinates[mesh]]][[1]], MeshCells[mesh, 2]]]]

which divides the region into some 3D mesh and then, for example:

r = DiscretizeGraphics[Polygon /@ Map[Reverse,
  Entity["Country", "France"]["Polygon"][[1, 1]], {2}]];

GeoPosition@GeoPositionXYZ@RegionCentroid@geoRegion3D@r

calculates the centroid of France. It suffices to replace the polygon of France with the world's polygon, which would take a looong time for calculation.

Now a different question popped into my mind which I wasn't able to find any ideas about. So after an hour of fruitless googling, I decided to try my luck in here.

There is a longitude on Earth which cuts all the land in half. i.e. there would be equal amounts of land on both hemispheres. Can we find this longitude using Mathematica?

Answer 1

Edit: Thanks to @polfosol for pointing out that using "Country" in GeoIdentify ignores Antarctica. That issue can be fixed by using "Ocean" instead of "Country". If one also wants to exclude lakes as being different from land, then I suppose that could be incorporated into the code. As with everything else, the answer depends on what exactly you want to know and exactly as to what information is available that approximates the conditions you desire.

Here is an alternative approach that uses random sampling.

(* Generate uniformly random points on the surface of a sphere with the Earth's radius *)
SeedRandom[12345];
n = 1000;
pts = RandomPoint[Sphere[{0, 0, 0}, 6.378*^6], n];

(* Get corresponding latitudes and longitudes *)
latlong = GeoPosition[GeoPositionXYZ[#]] & /@ pts;

(* Get the "status" for each point:  ocean (0) or not ocean (1) *)
status = If[GeoIdentify["Ocean", GeoPosition[#]] != {}, 0, 1] & /@ latlong;

(* Combine status and longitudes and sort by longitude *)
latlong = {#[[1]], #[[2]]} & /@ Flatten[latlong /. GeoPosition -> List, 1];
info = SortBy[Select[Transpose[Join[Transpose[latlong], {status}]], #[[3]] == 1 &][[All, {2, 3}]], #[[2]] &];

Note that the non-ocean areas on Earth comprise approximately 29% of the surface area and that the sampling has resulted in 28.8% of the points in non-ocean areas.

nPoints = Length[info]
(* 288 *)

Now figure out the proportion of non-ocean points between x and x+180 degrees longitude.

h = Join[{info[[nPoints]]}, info, info];
h[[1, 1]] = h[[1, 1]] - 360;
h[[nPoints + 2 ;; 2 nPoints + 1, 1]] = info[[All, 1]] + 360;
acc = N[Accumulate[h[[All, 2]]]/nPoints];
g = Interpolation[Transpose[{h[[All, 1]], acc}], InterpolationOrder -> 0];
t = Table[{long, g[long + 180] - g[long]}, {long, -180, 180, 0.1}];
ListPlot[{{{-150, 0}, {-150, 0.5}, {30, 0.5}, {30, 0}}, t}, Joined -> True]

We see 0.5 occurs at just two places 180 degrees of longitude apart: around 30 degrees East and 150 degrees West.

Answer 2

I accidentally did this for latitide first, so bear with me. This solution also isn't even close to scientifically rigorous (rectangular geographic projections aren't safe to use in this kind of calculation, but I did it anyway).

We're going to do this with some image processing. First, let's create a black-and-white image of the world.

ri = Binarize@
  GeoGraphics[GeoRange -> "World", 
   GeoBackground -> 
    GeoStyling[{"Coastlines", "Land" -> White, "Ocean" -> Black, 
      "Border" -> Black}]]

And now let's get the area of white pixels (components):

ComponentMeasurements[ri, "Area"] // Values // Total

46254.4

This output makes sense, since the total area of the image is 156800 (Times@@ImageDimensions@ri).

Now we're simply going to black out parts of the image and measure the area, and aim for half of the original area. We can do this with ReplacePixelValue, like so:

ReplacePixelValue[ri, {;; , 40 ;;} -> 0]

Our goal is 23127.2. A bit of trial and error later, we can get very close:

Total@Values@
  ComponentMeasurements[ReplacePixelValue[ri, {;; , 167 ;;} -> 0], 
   "Area"]

23120.9

So, 167 pixels up is our halfway area. We'll just convert this by eye to latitudes, by messing with GeoRange and getting an image that looks the same.

GeoGraphics[
 GeoBackground -> 
  GeoStyling[{"Coastlines", "Land" -> White, "Ocean" -> Black, 
    "Border" -> None}], GeoRange -> {{-90, 17}, All}]

So, 17 degrees north is the latitude that cuts the area of the land on Earth in half. More or less.

We can do the exact same thing for longitude.

goal = Total@Values@ComponentMeasurements[ri, "Area"]/2
f[x_] := goal - 
  Total@Values@
    ComponentMeasurements[ReplacePixelValue[ri, {x ;;, ;;} -> 0], 
     "Area"]

A bit of trial and error and we can get that the minimum difference is at column 324. We can then do the same thing, just eyeballing it with GeoGraphics:

ReplacePixelValue[ri, {324 ;;, ;;} -> 0]

GeoGraphics[
 GeoBackground -> 
  GeoStyling[{"Coastlines", "Land" -> White, "Ocean" -> Black, 
    "Border" -> None}], GeoRange -> {All, {-180, 28}}]

So. 28 degrees east. More or less. I hasten to add that I haven't had my coffee yet and a better version of this answer is welcome.

Answer 3

To build on & correct Carl Lange's image-processing answer, we can use ImageCrop instead of blacking out the pixels. We use a cylindrical equal-area projection for the world map to avoid exaggerating the amount of land area at high latitudes (note that Mathematica's default projection is equirectangular, which is not an equal-area projection.) To increase the precision of the final result, we also use a larger bitmap than the default.

eqari = Binarize@
  GeoGraphics[GeoRange -> "World", 
  GeoBackground -> 
    GeoStyling[{"Coastlines", "Land" -> White, "Ocean" -> Black, "Border" -> Black}], 
  GeoProjection -> "CylindricalEqualArea",
  ImageSize -> 15000]   

We then define

goal = Total@Values@ComponentMeasurements[eqari, "Area"]/2
halfwidth = Round[ImageDimensions[eqari][[1]]/2];
f[x_] := goal - 
  Total@Values@
    ComponentMeasurements[ImageCrop[eqari, {halfwidth, Full}, {-x, 0}], "Area"]

The syntax for this f[x] is a bit different. Here, because of the arguments for ImageCrop, a value of $x$ of -1 returns the land area of the Western Hemisphere while $x=+1$ corresponds to the Eastern Hemisphere. The longitude $\lambda$ of the western boundary of our hemisphere is related to $x$ by $\lambda = 90(x-1)$.

We can plot our new f (this takes some time to render; and I didn't bother to do it for the 15000-pixel-wide image; the image below is for the 10000-pixel-wide image):

Plot[f[x], {x, -1, 1}]

The desired point appears to be around $x \approx -0.65$. We can get it more precisely (though given the crudeness of this method, any precision beyond two decimal places is probably illusory):

root = FindRoot[f[x], {x, -0.7, -0.6}, Method -> "Brent"]
(* {x -> -0.6708} *)

This corresponds to a latitude of about

90 (x - 1) /. root
(* -150.372 *)

so the desired great circle consists of longitude lines at approximately 150°W and 30°E. Here's the map:

ImageCrop[eqari, {halfwidth, Full}, {-x, 0}] /. root

A width of 15000 pixels seems to be the practical limit of the image size on my machine (not sure why), but in principle this estimate could be further refined using larger bitmaps.

Answer 4

Thanks to the ideas proposed by @JimB, I managed to derive a faster and seemingly more accurate approach. Please make sure to upvote their answer first.

To get all the lands on Earth in Mathematica, it suffices to get the list of continents instead of countries. The only drawback is, weirdly enough, some countries in the Pacific Ocean aren't listed in any continent. Although if you type "countries of oceania" in WolframAlpha or as a free-form input, it returns a list of 25 countries including Australia. But in the list of continents, Australia is the fourth continent which is a region including the country of Australia and (part of) Papua New Guinea and doesn't contain the others!

Anyway, working with the list of continents has the advantage of getting a faster answer, and the area of ignored areas are fairly negligible (this can be fixed with a dirty hack which will be discussed later). Here is the code

(* Order of magnitude for number of random points *)
n = 1000;
(* get the list of continents *)
c = List @@ EntityList @ EntityClass["GeographicRegion", "Continents"];
(* Generate uniformly distributed random points on continents.
The number of points for each region is proportional to its total area *)
points = MapThread[RandomGeoPosition, {c, Floor[n Normalize[GeoArea /@ c]]}];

After generating random points, we can see them on the map like this:

GeoGraphics[Point[points], GeoProjection -> "Mollweide"]

(* Get the list of latitudes and longitudes of generated points *)
latlong = Flatten[#[[1]] & /@ points, 1];
(* Sort the points by their longitudes *)
info = SortBy[If[Last[#] > 180, # - {0,180}, #]& /@ latlong, #[[2]]&];

This function returns the difference between the number of points inside the longitude range of [x ... 180+x] and outside this range:

f[x_]:= With[{t= CountsBy[info, #[[2]]> x && #[[2]]<= x+180 &]}, First@t - Last@t]
ListPlot[Table[Abs@f[x], {x, -180, 0}], Joined-> True, DataRange-> {-180, 0}]

The above plot shows that for approximately 100000 points, the difference is minimum around 150 degrees west. We can get a more precise measure by reducing the domain of f around 150 with a smaller step:

ListPlot[Table[Abs@f[x], {x,-151,-149,.05}], Joined-> True, DataRange-> {-151, -149}]

The ignored lands of Oceanea

One can join the list of these countries to the list of continents and then generate random points:

(* List of countries in Oceania, excluding Australia and Papua New Guinea *)
oc = Delete[List@@EntityList@EntityClass["Country", "Oceania"], {{2}, {17}}];
world = Join[c, oc];
points = MapThread[RandomGeoPosition,
  {world, Floor[n Normalize[GeoArea /@ world]]}];

This would barely modify the final result.