Tuesday 19 November 2019

Calculating Centroid of all lands of Earth using SQL Server Spatial?


I want to find the center point of all lands of the earth.


What I have tried up to now is downloading shapefile of all countries and polar parts of the earth. Then I imported that shapefile into a SQL Server 2008 database and then tried to run some queries to find that point but all of them were in vain.


How can I calculate the Centroid of all lands of the earth using SQL Server?



Answer



This perhaps is most easily done with raster calculations, because the centroid of the land surface is found (by definition) by averaging the geocentric coordinates (X, Y, Z), weighting them by land area. Averaging is done using zonal statistics.


To obtain these coordinates and the area scale factor, compute the following grids in geographic coordinates covering all the land areas at any resolution desired (I used 0.1 degree cells, which form a 3600 by 1800 grid covering the entire earth):





  • [Longitude] is the first coordinate.




  • [Latitude] is the second coordinate.




  • [X] is Cos([Longitude]) * Cos([Latitude]).





  • [Y] is Sin([Longitude]) * Cos([Latitude]).




  • [Z] is Sin([Latitude]).




  • [F] is Cos([Latitude]). This is the area scale factor.





(Note that it is unnecessary to multiply any of [X], [Y], or [Z] by the earth's radius. In effect we are computing the centroid of a perfectly scaled spherical model of the earth. The position of that centroid on the model will be the scaled position of the desired centroid.)


Compute zonal means of [X] * [F], [Y] * [F], and [Z] * [F]. When I did this with a shapefile showing all major land areas (including Antarctica), with longitude ranging from 0 to 360 degrees, the zonal means were -0.1204, -0.0661, and 0.1356, respectively.


(Alternatively you can reproject [X], [Y], and [Z] into any equal-area projection and compute their zonal means directly without resorting to the area weighting by [F]. Do whichever is easier in your GIS.)


As we should expect, this centroid lies much closer to the earth's center at (0, 0, 0) than it does to its surface. If we wish to map it we can compute its latitude and longitude and ignore its depth. The longitude of a point (X, Y, Z) is determined by (X, Y) alone and usually is computed as their arctangent; the latitude is determined by Sqrt(X^2+Y^2) and Z and usually is computed as their arctangent. I obtain the location (208.77, 44.63) (which is 28.77 degrees east longitude, 44.63 degrees north latitude). This lies near the eastern coast of Romania.


World centroid map


This map is in geographic coordinates: the horizontal coordinate is longitude and the vertical coordinate is latitude. The background depicts the X coordinate grid. The centroid is located beneath the red marker at the center right.


If that doesn't look quite like the center of land masses to you, inspect it from a different location, such as straight overhead:


Orthographic projection


The invisible portions of the United States, Central, and South America will more than counterbalance the visible landmasses; Antarctica and Australia restore the balance. This result looks plausible to me. Your results may vary a little depending on the scale and accuracy of your land mass data, the resolution of your grid, and whether you use a spherical or ellipsoidal model of the earth. (The area factor [F] requires a more involved computation in an ellipsoidal model.)





The advantage of a grid-based solution becomes evident in generalizations of this question. For instance, by replacing [F] (the "land amount" grid) by a grid representing the total amount of something else per cell, you can compute weighted centroids. If you have a population density grid [P] then [F] * [P] is proportional to total population; using [F] * [P] in place of [F] in these calculations would give a world population centroid.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...