12
$\begingroup$

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?

$\endgroup$
11
$\begingroup$

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.

Proportion of land between x and x+180 degrees longitude

$\endgroup$
5
  • $\begingroup$ Very nice, I really like this. $\endgroup$ – Carl Lange Jun 21 at 21:46
  • $\begingroup$ GeoIdentify["Country", GeoPosition[Entity["GeographicRegion", "Antarctica"]]] returns empty, which means Antarctica and possibly some other land on Earth are ignored in this method. On the other hand, GeoIdentify does count some lakes as land. $\endgroup$ – polfosol Jun 22 at 10:03
  • $\begingroup$ @polfosol Thanks! I'll change "Country" to "Ocean" (and modify the If statement) and add in a warning that it depends on what is considered "land". $\endgroup$ – JimB Jun 22 at 16:08
  • $\begingroup$ You are welcome! By the way, your method of generating random points and then checking them is way too slow on my system. I was just wondering have you tried RandomGeoPosition? It seems to be more straightforward and does the same job faster. $\endgroup$ – polfosol Jun 22 at 16:26
  • $\begingroup$ @polfosol Good suggestion. I'll try that. I know about sampling but next to nothing about GIS and the associated Mathematica functions. But the slowdown for me is the GeoIdentify function. $\endgroup$ – JimB Jun 22 at 16:46
7
$\begingroup$

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

enter image description here

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]

enter image description here

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

enter image description here

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

enter image description here

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]

enter image description here

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

enter image description here

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.

$\endgroup$
9
  • $\begingroup$ This was a very nice try for another, completely different question! +1 by the way $\endgroup$ – polfosol Jun 21 at 12:11
  • $\begingroup$ If you use the correct projection, a similar technique would work reasonably accurately $\endgroup$ – mikado Jun 21 at 12:14
  • $\begingroup$ @polfosol I wouldn't say it's completely different, just (initiall) in the incorrect dimension - I've added longitude now. Unless I've completely misunderstood the question, which is quite plausible! :) $\endgroup$ – Carl Lange Jun 21 at 12:28
  • $\begingroup$ @mikado What projection would you suggest as more correct in this case? $\endgroup$ – Carl Lange Jun 21 at 12:28
  • 1
    $\begingroup$ The portion of the Earth spanning 180°W to 28°E is greater than a hemisphere. To get an actual hemisphere, you should be looking at the span from longitude $\alpha - 180$ to longitude $\alpha$ (with positive values meaning "east" and negative values meaning "west".) $\endgroup$ – Michael Seifert Jun 21 at 14:12
6
$\begingroup$

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

enter image description here

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

enter image description here

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.

$\endgroup$
4
  • $\begingroup$ I think you need to use a cylindrical projection so that you account for variation of area with latitude. I'm not sure you are doing so. $\endgroup$ – mikado Jun 21 at 14:47
  • $\begingroup$ @mikado: Good point; the default projection Mathematica uses is equirectangular, not equal-area. Let me see if I can fix this. $\endgroup$ – Michael Seifert Jun 21 at 14:55
  • $\begingroup$ Nice, I love it! Much better thought out than my answer :) $\endgroup$ – Carl Lange Jun 21 at 15:25
  • $\begingroup$ I think you can speed up the plot of f by finding the column totals, duplicating the column totals, and then using Accumulate: counts = Total[#] & /@ Transpose[ImageData[eqari]]; counts = counts/Total[counts]; counts = Join[counts, counts]; acc = Accumulate[counts]; g = Interpolation[acc, InterpolationOrder -> 0]; f = {2 #/7500 - 3, g[# + 7500] - g[#] - 0.5} & /@ Range[7500, 15000]; ListPlot[f, PlotRange -> {{-1, 1}, Automatic}]. $\endgroup$ – JimB Jun 21 at 23:31
1
$\begingroup$

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

random points on world

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

plot of f(x)

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

zoom in

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.

$\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.