Previously, we&#8217;ve demonstrated using the <a href="https://www.sisense.com/blog/calculating-distance-between-data-centers-on-a-globe/">Haversine Formula to calculate distances</a> on a globe in SQL. Today we&#8217;ll extend that to look at areas for irregular polygons defined by latitude-longitude coordinates.



Computing the area for polygons on a globe is challenging for a number of reasons:



<ul><li>A polygon drawn on the Earth&#8217;s surface cannot be computed using simple Pythagorean math. To deal with this, we&#8217;ll project our latitude-longitude coordinates to a plane.</li><li>For polygons with more than three points, the drawing order is important! The correct area is calculated using the coordinates of the points in relation to their neighbors. An improperly ordered table of polygon vertices will return wildly inaccurate results.</li></ul>



<figure class="wp-block-image fancybox"><img decoding="async" src="https://cdn.sisense.com/wp-content/uploads/Red-globe.png" alt="Sinusoidal globe projection" class="wp-image-75315"/></figure>



The Data



We will use a table &#8220;points&#8221; containing several lat-long points in the Bay Area that will be the vertices of our polygon:



<figure class="wp-block-image fancybox"><img decoding="async" src="https://cdn.sisense.com/wp-content/uploads/CA-cities-table.png" alt="California cities table" class="wp-image-75320"/></figure>



Which looks like this when plotted on Google Maps:



<figure class="wp-block-image fancybox"><img decoding="async" src="https://cdn.sisense.com/wp-content/uploads/Google-Map.png" alt="Google maps image" class="wp-image-75325"/></figure>



<h2 class="wp-block-heading">Converting Latitude and Longitude to (x, y) Coordinates</h2>



First we&#8217;ll convert our latitude and longitude to Cartesian coordinates. There are a lot of options when projecting a sphere to a flat surface — each of which has tradeoffs depending on whether you want to preserve area, shape or <a href="https://www.wired.com/2013/07/projection-mercator/" target="_blank" rel="noreferrer noopener" aria-label=" (opens in a new tab)">be able to plot a ship&#8217;s course on a constant bearing</a> with a straight line.



For our purpose, we want an equal-area projection. We&#8217;ll use a <a href="https://en.wikipedia.org/wiki/Sinusoidal_projection" target="_blank" rel="noreferrer noopener" aria-label=" (opens in a new tab)">sinusoidal projectio</a>n, which can be projected from lat-long to Cartesian coordinates, and is defined by:



<pre class="wp-block-code"><code>x = (longitude - prime_meridian) * cosine(latitude)
y = latitude</code></pre>



Here is a picture of a sinusoidal projection. Please don&#8217;t use this for navigation:



<figure class="wp-block-image fancybox"><img decoding="async" src="https://cdn.sisense.com/wp-content/uploads/Sinusoidal-globe-min.png" alt="Sinusoidal globe" class="wp-image-75352"/></figure>



We&#8217;ll convert our map into a new unit more appropriate for measuring area in the process. I&#8217;ll use kilometers, which is factored into our conversion through the number 6371, the Earth&#8217;s radius in km.



When we project our points to x,y Cartesian coordinates, we&#8217;ll grab a few more useful numbers:



<ul><li>The geographic mean avgx, avgy of the vertices using a window function. We&#8217;ll need this later to define pathing for our polygon; i.e., the order in which the vertices are connected.</li><li>Our points normalized to the mean, so we can do some math centered around 0,0.</li></ul>



<pre class="wp-block-code"><code>with sinusoidal as (
 select
 -- y-axis variables
 lat * 3.14159 * 6371 / 180 as y
 , avg(lat) over() * 3.14159 * 6371 / 180 as avgy
 , avg(lat) over() * 3.14159 * 6371 / 180 - (lat * 3.14159 * 6371 / 180) as normy
 -- x-axis variables
 , lng * 3.14159 * 6371 / 180 * cos(radians(lat)) as x
 , avg(lng) over() * 3.14159 * 6371 / 180 * cos(radians(avg(lat) over())) as avgx
 , avg(lng) over() * 3.14159 * 6371 / 180 * cos(radians(avg(lat) over())) - (lng * 3.14159 * 6371 / 180 * cos(radians(lat))) as normx

 from
 points)</code></pre>



<h2 class="wp-block-heading">Drawing a Polygon</h2>



We need to work out the pathing between the vertices. We are going to compute the area using the angles between points, so knowing which points come before and after the current point is critical. To think of this another way, four vertices could define both of the following shapes:



<figure class="wp-block-image fancybox"><img decoding="async" src="https://cdn.sisense.com/wp-content/uploads/purple-shapes.png" alt="Purple shapes" class="wp-image-75335"/></figure>



We want our path to follow the outside perimeter of our polygon, touching each vertex and never intersecting itself. One way to do this is to compute the angle between the geographic mean of our polygon and each point, and then order our table by the angle:



<figure class="wp-block-image fancybox"><img decoding="async" src="https://cdn.sisense.com/wp-content/uploads/purple-square.png" alt="Purple square" class="wp-image-75341"/></figure>



We can determine our angles using atan, like so:



<pre class="wp-block-code"><code>, angles as (
 select
 *
 , (atan(normy / normx) * 180) / 3.14159 
 + ((normx / abs(normx)) * -1) * 90 + 180 
 as angle
 from
 sinusoidal)</code></pre>



The table now looks like:



<figure class="wp-block-image fancybox"><img decoding="async" src="https://cdn.sisense.com/wp-content/uploads/flat-table.png" alt="Flat table" class="wp-image-75346"/></figure>



Now we have a column to order our vertices and thus define our polygon!



Each vertex will need to have values calculated between itself and its neighbors in a complete loop. We can union all the table to itself so the first row can look at the last row using a lag function. An &#8220;ordr&#8221; column can be used to grab the second half of the table later, which will have fully populated data.



<pre class="wp-block-code"><code>, unioned as (
select
 1 as ordr
 , *
from 
 angles
union all
select
 2
 , *
from 
 angles )</code></pre>



Ok, our last CTE!&nbsp;



We use the lag window function to append the neighboring vertex&#8217;s x,y coordinates. I&#8217;ve slightly modified this query and my final equation from the standard form to avoid using both a lead and a lag function, instead opting to use lag(x,1)lag(y,2)order bylag and . This is the step where table order is critical, and we need to&nbsp; both our ordr and angle columns to to the correct data.



<pre class="wp-block-code"><code>, final as (
select
 *
 , lag(y, 2) over(order by ordr, angle) as lagy2
 , lag(x, 1) over(order by ordr, angle) as lagx
from
 unioned
order by
 ordr
 , angle)</code></pre>



With our table fully populated, we can now calculate our polygon&#8217;s area. More information on this algorithm, as well as an implementation in C++, can be found <a href="http://geomalgorithms.com/a01-_area.html" target="_blank" rel="noreferrer noopener" aria-label=" (opens in a new tab)">here</a>.



Any projection and measurement will have errors, related to both the projection&#8217;s limitations and the Earth being an oblate spheroid rather than a true sphere, so it&#8217;s important to validate any geospatial algorithm you use in code.&nbsp;



In this case. our query returns a result of 766.93 sq. km, and Google Earth returns 767 sq. km (without the option to see decimals) — making this method very robust at the city-regional scale!&nbsp;

Geographic Analysis in SQL: Measuring Polygon Area from Latitude and Longitude

SQL Superstar

LinkedIn

Twitter

GitHub

curve-image-unique-image-unique

curve

3-dark-2-image-unique-image-unique

3 DARK 2

Get the latest in analytics right in your inbox.

Article