[GIS] Converting Longitude and Latitude Coordinates to Square Miles

arealatitude longitudepolygonspherical-geometry

I have a set of Long and Lat co-ordinates making up a five point polygon which I have converted in excel using the following formula:

=((X1*Y2)-(X2*Y1)+(X2*Y3)-(X3*Y2)+(X3*Y4)-(X4*Y3)+(X4*Y5)-(X5*Y4)+(X5*Y1)-(X1*Y5))*0.5

*Replacing X,Y with cell references of course.

This gives me a negative number but how do I convert this number to square miles or kilometres?

Best Answer

Try this formula (assuming your source is WGS1984, if not then you'll need to adjust the ellipsoid used by the second line):

area = rad(x2 - x1) * (2 + sin(rad(y1)) + sin(rad(y2))) + rad(x3 - x2) * (2 + sin(rad(y2)) + sin(rad(y3))) + rad(x4 - x3) * (2 + sin(rad(y3)) + sin(rad(y4))) + rad(x5 - x4) * (2 + sin(rad(y4)) + sin(rad(y5)))

area = abs(area * 6378137.0 * 6378137.0 / 2.0)

rad() is a function that converts Degrees to Radians (i.e. Degrees * PI / 180)

Source: OpenLayers LinearRing

This will result in an area returned in square meters.

Related Solutions

PostGIS – Searching planet_osm_point by Longitude and Latitude

To search for points around that location in London, you'd want a query like this

SELECT * FROM planet_osm_polygon WHERE ST_Intersects(way, ST_Buffer(ST_Transform(ST_SetSRID(ST_Point(-0.12,51.5),4326), 3857), 500); This will take the point, transform it into the web mercator projection (EPSG 3857) used by default in osm2pgsql, buffer it by 500 mercator meters, and find points that lie within that polygon.

A more accurate way to do this would be

SELECT * FROM planet_osm_polygon WHERE ST_Intersects(way, ST_Transform(ST_Buffer( ST_SetSRID(ST_Point(-0.12,51.5),4326)::geography,500)::geometry, 3857) );

This converts to geography so you can expand the point by 500 true meters.

A third, and probably quickest way, is

SELECT * FROM planet_osm_polygon WHERE ST_DWithin(way, ST_Transform(ST_SetSRID(ST_Point(-0.12,51.5),4326), 3857), 500);

Of course in practice you'll have some other WHERE conditions to apply, and be using a more complicated real-world geometry to search against.

Some common mistakes are

Trying to compare geometries of different projections
Reprojecting the geometry in the planet_osm_point table, preventing index usage
Mixing up latitude and longitude.

It's worth noting that PostGIS geometries are completely distinct from PostgreSQL geometric types. The appropriate page to see available functions is Chapter 8 of the PostGIS reference.

SQL Server – Converting NZTM to Latitude/Longitude

The following SQL query is an implementation of the formula provided by LINZ to convert from NZTM to Geographic (Lat/Long). Source

This example uses the centre of the amphitheater at Te Papa as a reference point.

Input: Northing 5427502.0, Easting 1749165.0 Output: Longitude -41.290153, Latitude 174.781428

Begin 

Declare
-- Constants
@a      float = 6378137,
@f      float = 1/298.257222101,
@phizero float =    0,
@lambdazero float = 173,
@Nzero  float = 10000000,
@Ezero  float = 1600000,
@kzero  float = 0.9996,

-- Common variables
@b  float = 0,
@esq float  = 0,
@mzero float = 0,
@A0 float = 0,
@A2 float = 0,
@A4 float = 0,
@A6 float = 0,

--input Northing (Y), Easting (X) variables 
@N      float = 5427502.0, 
@E      float = 1749165.0, 

-- Variables
@Nprime float = 0,
@mprime float = 0,
@smn    float = 0,
@G      float = 0,
@sigma  float = 0,
@phiprime   float = 0,
@rhoprime   float = 0,
@upsilonprime float = 0,
@psiprime   float = 0,
@tprime float = 0,
@Eprime float = 0,
@chi    float = 0,
@term_1 float = 0,
@term_2 float = 0,
@term_3 float = 0,
@term_4 float = 0,
@term1  float = 0,
@term2  float = 0,
@term3  float = 0,
@term4  float = 0,

--Output
@latitude   float = 0, 
@longitude  float = 0; 

--Calculation: To NZTM to NZGD2000/WGS84

--Common

SET @b          = @a*(1-@f)
SET @esq        = 2*@f-POWER(@f,2)
SET @A0         =1-@esq/4-3*POWER(@esq,2)/64-5*POWER(@esq,3)/256 
SET @A2         =0.375*(@esq+POWER(@esq,2)/4+15*power(@esq,3)/128) 
SET @A4         =15*(POWER(@esq,2)+3*POWER(@esq,3)/4)/256  
SET @A6         =35*POWER(@esq,3)/3072 

SET @Nprime     = @N-@Nzero
SET @mprime     = @Nprime/@kzero
SET @smn        = (@a-@b)/(@a+@b)
SET @G          = @a*(1-@smn)*(1-POWER(@smn,2))*(1+9*POWER(@smn,2)/4+225*POWER(@smn,4)/64)*PI()/180.0
SET @sigma      = @mprime*PI()/(180*@G)
SET @phiprime   = @sigma+(3*@smn/2-27*POWER(@smn,3)/32)*SIN(2*@sigma)+(21*POWER(@smn,2)/16-55*POWER(@smn,4)/32)*SIN(4*@sigma)+(151*POWER(@smn,3)/96)*SIN(6*@sigma)+(1097*POWER(@smn,4)/512)*SIN(8*@sigma)
SET @rhoprime   = @a*(1-@esq)/POWER((1-@esq*POWER((SIN(@phiprime)),2)),1.5)
SET @upsilonprime   =@a/SQRT(1-@esq*POWER((SIN(@phiprime)),2)) 

SET @psiprime   = @upsilonprime/@rhoprime
SET @tprime     = TAN(@phiprime)
SET @Eprime     = @E-@Ezero
SET @chi        = @Eprime/(@kzero*@upsilonprime)
SET @term_1     = @tprime*@Eprime*@chi/(@kzero*@rhoprime*2)
SET @term_2     = @term_1*POWER(@chi,2)/12*(-4*POWER(@psiprime,2)+9*@psiprime*(1-POWER(@tprime,2))+12*POWER(@tprime,2))
SET @term_3     = @tprime*@Eprime*POWER(@chi,5)/(@kzero*@rhoprime*720)*(8*POWER(@psiprime,4)*(11-24*POWER(@tprime,2))-12*POWER(@psiprime,3)*(21-71*POWER(@tprime,2))+15*POWER(@psiprime,2)*(15-98*POWER(@tprime,2)+15*POWER(@tprime,4))+180*@psiprime*(5*POWER(@tprime,2)-3*POWER(@tprime,4))+360*POWER(@tprime,4))
SET @term_4     = @tprime*@Eprime*POWER(@chi,7)/(@kzero*@rhoprime*40320)*(1385+3633*POWER(@tprime,2)+4095*POWER(@tprime,4)+1575*POWER(@tprime,6))
SET @term1      = @chi*(1/COS(@phiprime))
SET @term2      = POWER(@chi,3)*(1/COS(@phiprime))/6*(@psiprime+2*POWER(@tprime,2))
SET @term3      = POWER(@chi,5)*(1/COS(@phiprime))/120*(-4*POWER(@psiprime,3)*(1-6*POWER(@tprime,2))+POWER(@psiprime,2)*(9-68*POWER(@tprime,2))+72*@psiprime*POWER(@tprime,2)+24*POWER(@tprime,4))
SET @term4      = POWER(@chi,7)*(1/COS(@phiprime))/5040*(61+662*POWER(@tprime,22)+1320*POWER(@tprime,4)+720*POWER(@tprime,6))

SET @latitude   = (@phiprime-@term_1+@term_2-@term_3+@term_4)*180/PI()
SET @longitude  = @lambdazero+180/PI()*(@term1-@term2+@term3-@term4)

print N' @a:        ' + STR( @a         ,10,8) 
print N' @b:        ' + STR( @b         ,10,8) 
print N' @esq:      ' + STR( @esq       ,10,8) 
print N' @A0:       ' + STR( @A0        ,10,8) 
print N' @A2:       ' + STR( @A2        ,10,8) 
print N' @A4:       ' + STR( @A4        ,10,8) 
print N' @A6:       ' + STR( @A6        ,15,13) 
print N' @Nprime:   ' + STR( @Nprime    ,10,8) 
print N' @mprime:   ' + STR( @mprime    ,10,8) 
print N' @smn:      ' + STR( @smn       ,10,8) 
print N' @G:        ' + STR( @G         ,10,8) 
print N' @sigma:    ' + STR( @sigma     ,10,8) 
print N' @phiprime: ' + STR( @phiprime  ,10,8) 
print N' @rhoprime: ' + STR( @rhoprime  ,10,8) 
print N' @upsilonprime: ' + STR( @upsilonprime ,10,8) 
print N' @psiprime: ' + STR( @psiprime  ,10,8) 
print N' @tprime:   ' + STR( @tprime    ,10,8) 
print N' @Eprime:   ' + STR( @Eprime    ,10,8) 
print N' @chi:      ' + STR( @chi       ,10,8) 
print N' @term_1:   ' + STR( @term_1    ,10,8) 
print N' @term_2:   ' + STR( @term_2    ,10,8) 
print N' @term_3:   ' + STR( @term_3    ,10,8) 
print N' @term_4:   ' + STR( @term_4    ,10,8) 
print N' @term1:    ' + STR( @term1     ,10,8) 
print N' @term2:    ' + STR( @term2     ,10,8) 
print N' @term3:    ' + STR( @term3     ,10,8) 
print N' @term4:    ' + STR( @term4     ,10,8) 

-- annnddd here we are :)   
print N' @latitude: ' + STR( @latitude  ,10,8) 
print N' @longitude:' + STR( @longitude ,10,8) 

End

Another example uses the centre of the amphitheatre at Te Papa as a reference point.

Input: Longitude -41.290153, Latitude 174.781428 Output: Northing 5427502.0, Easting 1749165.0

Begin 
Declare
-- constants
@a      float = 6378137,
@f      float = 1/298.257222101,
@phizero float =    0,
@lambdazero float = 173,
@Nzero  float = 10000000,
@Ezero  float = 1600000,
@kzero  float = 0.9996,

-- Common 
@b  float = 0,
@esq float  = 0,
@mzero float = 0,
@A0 float = 0,
@A2 float = 0,
@A4 float = 0,
@A6 float = 0,

--input Latitude (Y), Longitude (X) and varibles 

@phi    float = -41.290153, 
@lambda float = 174.781428,

-- Variables
@rho    float = 0,
@upsilon float = 0,
@psi    float = 0,
@t      float = 0,
@omega  float = 0,
@m      float = 0,
@t_1    float = 0,
@t_2    float = 0,
@t_3    float = 0,
@t_4    float = 0,
@t1     float = 0,
@t2     float = 0,
@t3     float = 0, 

--Output
@Easting    float = 0, 
@Northing   float = 0, 
@latitude   float = 0, 
@longitude  float = 0
; 

--Calculation: NZGD2000/WGS84 to NZTM

--Common

SET @b          = @a*(1-@f)
SET @esq        = 2*@f-POWER(@f,2)
SET @A0         = 1-@esq/4-3*POWER(@esq,2)/64-5*POWER(@esq,3)/256 
SET @A2         = 0.375*(@esq+POWER(@esq,2)/4+15*power(@esq,3)/128)  
SET @A4         = 15*(POWER(@esq,2)+3*POWER(@esq,3)/4)/256  
SET @A6         = 35*POWER(@esq,3)/3072     

-- Calculation  To NZTM
SET @rho        = @a*(1-@esq)/POWER(1-@esq*POWER((SIN(PI()*@phi/180)),2),1.5) 
SET @upsilon    = @a/SQRT(1-@esq*POWER((SIN(PI()*@phi/180)),2))
SET @psi        = @upsilon/@rho
SET @t          = TAN(@phi*PI()/180)
SET @omega      = (@lambda-@lambdazero)*PI()/180
SET @m          = @a*(@A0*PI()*@phi/180-@A2*SIN(2*PI()*@phi/180)+@A4*SIN(4*PI()*@phi/180)-@A6*SIN(6*PI()*@phi/180))

SET @t_1        = POWER(@omega,2)*@upsilon/2*SIN(@phi*PI()/180)*COS(@phi*PI()/180)
SET @t_2        = POWER(@omega,4)*@upsilon/24*SIN(@phi*PI()/180)*POWER((COS(@phi*PI()/180)),3)*(4*POWER(@psi,2)+@psi-POWER(@t,2))
SET @t_3        = POWER(@omega,6)*@upsilon/720*SIN(PI()*@phi/180)*POWER((COS(PI()*@phi/180)),5)*(8*POWER(@psi,4)*(11-24*POWER(@t,2))-28*POWER(@psi,3)*(1-6*POWER(@t,2))+POWER(@psi,2)*(1-32*POWER(@t,2))-@psi*2*POWER(@t,2)+POWER(@t,4))
SET @t_4        = POWER(@omega,8)*@upsilon/40320*SIN(PI()*@phi/180)*POWER((COS(PI()*@phi/180)),7)*(1385-3111*POWER(@t,2)+543*POWER(@t,4)-POWER(@t,6))
SET @t1         = POWER(@omega,2)/6*POWER((COS(PI()*@phi/180)),2)*(@psi-POWER(@t,2))
SET @t2         = POWER(@omega,4)/120*POWER((COS(PI()*@phi/180)),4)*(4*POWER(@psi,3)*(1-6*POWER(@t,2))+POWER(@psi,2)*(1+8*POWER(@t,2))-2*@psi*POWER(@t,2)+POWER(@t,4))
SET @t3         = POWER(@omega,6)/5040*POWER((COS(PI()*@phi/180)),6)*(61-479*POWER(@t,2)+179*POWER(@t,4)-POWER(@t,6))

SET @Easting    = @Ezero+@kzero*@upsilon*@omega*COS(PI()*@phi/180)*(1+@t1+@t2+@t3)
SET @Northing   = @Nzero+@kzero*(@m+@t_1+@t_2+@t_3+@t_4)

print N' @b:        ' + STR( @b         ,12,10) 
print N' @esq:      ' + STR( @esq       ,12,10) 
print N' @A0:       ' + STR( @A0        ,12,10) 
print N' @A2:       ' + STR( @A2        ,12,10) 
print N' @A4:       ' + STR( @A4        ,12,10) 
print N' @A6:       ' + STR( @A6        ,12,10) 
print N' @rho:      ' + STR( @rho       ,12,10) 
print N' @upsilon:  ' + STR( @upsilon   ,12,10) 
print N' @psi:      ' + STR( @psi       ,12,10) 
print N' @t:        ' + STR( @t         ,12,10) 
print N' @omega:    ' + STR( @omega     ,12,10) 
print N' @m:        ' + STR( @m         ,12,10) 
print N' @t_1:      ' + STR( @t_1       ,12,10) 
print N' @t_2:      ' + STR( @t_2       ,12,10) 
print N' @t_3:      ' + STR( @t_3       ,12,10) 
print N' @t_4:      ' + STR( @t_4       ,12,10) 
print N' @t1:       ' + STR( @t1        ,12,10) 
print N' @t2:       ' + STR( @t2        ,12,10) 
print N' @t3:       ' + STR( @t3        ,12,10) 

print N' @Easting:  ' + STR( @Easting   ,12,10) 
print N' @Northing: ' + STR( @Northing  ,12,10) 

End

Best Answer

Related Solutions

PostGIS – Searching planet_osm_point by Longitude and Latitude

SQL Server – Converting NZTM to Latitude/Longitude

Related Question