4

I am not normally a GIS developer, I am just helping out with a bug and I am stumped on what is going on here, could I get some help please.

I have two calculations, one in 4326 (WGS84) and one in 3857 (Google Mercator). The coordinates were converted from 4326 -> 3857 and I have confirmed those are correct with ArcGIS and various online site conversion tools. The results I have are...

-114.50,38.0=-12746081.696,4579425.813

-114.00,40.0=-12690421.950,4865942.280

If I do an estimate of the area for 3857 by calculating the difference in the x and y points then multiplying the differences and dividing by 1 mil to get sq km, I get the result just under 16,000 sq km. Running SqlGeometry.STArea on a polygon bbox using these points this gives me a similar number.

If I convert the lat long to UTM I receive the numbers ...

-114.50,38.0=719510.335753775,4208764.46330096

-114.00,40.0=243900.352021924,4432069.05663229

Excel Calculations
Easting         Northing        sq m                    sq km
719,510.3358    4,208,764.4633      
243,900.3520    4,432,069.0566  
--------------------------------------------------------------------
475,609.9837    223,304.5933    106,205,894,001.5610    106,205.8940

Sites used to convert to UTM

http://www.latlong.net/lat-long-utm.html

http://www.rcn.montana.edu/Resources/Converter.aspx

http://www.earthpoint.us/Convert.aspx

Running the same calculation of the difference and the area gives me the result of 106,205 sq km.

The code currently calculates the lat long by doing the following on a polygon of the bbox ...

SqlGeography geog = SqlGeography.STGeomFromWKB(aoi.STAsBinary(), (int)aoi.STSrid);
area = (double)geog.STArea();

This gives me about 9,600 sq km.

What I need to know...

I am trying to find out which calculation is correct. I am leaning toward 3857 being correct as I the calculations give me approximatly the same numbers no matter how I do them. Since I feel that 4326 is incorrect, I would need a better way to calculate this in the current application (uses C# and SQLGeometry). If someone sees an error in my calculations please let me know and how to do this better if so.

Improve this question
14
  • 1
    Neither. You can't scale areas from degrees to hectares (square miles etc..). If you want an area in metres (or feet) project to a suitable projected coordinate system and calculate the area there. WGS84 UTM has zones that cover the whole world; you can determine from the coordinates whether it's north or south of equator and which zone it's in.
    – Michael Stimson
    Commented Sep 17, 2015 at 21:57
  • 1
    Yes, but probably not one so distorted. 3587 is good because it covers the whole world (like WGS84) but flattening out a sphere is going to lead to severe distortion away from the centre point. A projection like UTM still distorts but it does so from a central meridian at intervals of 6 degrees so suffers much less distortion from being flattened.
    – Michael Stimson
    Commented Sep 17, 2015 at 22:09
  • 2
    You may also have a problem that the polygon doesn't have enough vertices to adequately reflect its shape in a projected coordinate reference system (ProjCRS). That is, if the longitude line might be curved in the ProjCRS, then two points aren't enough.
    – mkennedy
    Commented Sep 17, 2015 at 23:02
  • 3
    The 106km2 value is wrong because your UTM coordinates are wrong. Whatever you're using is treating each coordinate separately. The two longitude values are in different zones which is why the 'left' value has a X value of 700k while the right one is 200k. The true separation in X is more like 44km. Also it's not a rectangle in UTM so a quick rectangular calc will have errors too.
    – mkennedy
    Commented Sep 17, 2015 at 23:06
  • 3
    EPSG::3857 does NOT maintain areas so the 16km2 calculation is wrong. If you instead project it into UTM (still not equal area but designed for the data's location) or an equal area projection, the areas are 9621 km2 versus 9615 km2 for an undensified polygon.
    – mkennedy
    Commented Sep 17, 2015 at 23:08

1 Answer 1

Reset to default
2

Apart from the fact that the EPSG:3857 projection does not preserve areas, the units of it are not real meters. Only at the equator they match real meters, while they get stretched the more to the poles you go.

UTM does not preserve areas either, but it minimizes distance errors as long as you are within the same UTM zone (which is not the case in the example). So it uses real meters as units.

The only correct way to calculate areas is using an equal area projection. For Northern America, you can try ESRI:102008

+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs

Using QGIS and the built-in $area function of the field calculator, I get for your 0.5x2.0 degree area:

EPSG:2163  8105.49 km² (based on a sphere)
EPSG:3857 15947.43 km² (google square meters)
EPSG:32107 9617.65 km² (NAD83 Nevada East)
EPSG:32611 9622.15 km² (WGS 84 / UTM zone 11N)
ESRI:12008 9616.38 km² (North_America_Albers_Equal_Area_Conic)

As you see, EPSG:3857 is worth for nothing. It is just a popular visualization, nothing more.

Improve this answer
5
  • With the understanding that the meters are stretched, I still see results that are too widely different from what they should be. I would expect a small % difference with the stretching. I am seeing with the SQLGeometry results about 60% reduction in the area 15.9K sq km vs 9.6K sq km. That seems more like an error in the SQLGeometry results (could be user input). With my manual calculations I see 15.9K vs 106K sq km. That is almost certainly an incorrect calculation somewhere.
    – JabberwockyDecompiler
    Commented Sep 18, 2015 at 14:39
  • The UTM calculation across two zones (11N and 12N) is obviously wrong. If you shift your area half a degree to the East or West, you will be within one UTM zone and get reasonable values. For the others, see my extended answer.
    – AndreJ
    Commented Sep 18, 2015 at 15:59
  • Mercator projection maintains angles and shapes not distances nor areas. The linear distortion is 1.269 at 39N.
    – mkennedy
    Commented Sep 18, 2015 at 16:06
  • The updated numbers are very helpful to see. From what I understand from your updates, the 15K km² is incorrect and the 9.6K km² is correct. In that case I should be converting the 3857 to something else to calc area. Does that sound correct?
    – JabberwockyDecompiler
    Commented Sep 18, 2015 at 16:53
  • You can take the original EPSG:4326 degree values and transform them to something useful.
    – AndreJ
    Commented Sep 18, 2015 at 17:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.