Understanding Differences in Area Calculations between SQL Server and QGIS

Question

I'm trying to get accurate calculations of area from SQL Server. When I perform comparisons between it and QGIS, I'm finding that my numbers do not match.

Here's what I'm doing in SQL Server and QGIS, and the results I find from each:

0. Environment

First, let me state what my environment is:

1. SQL Server: I'm running SQL Server 2016 Developer Edition. SELECT @@VERSION gives this version string:

   Microsoft SQL Server 2016 (SP1-GDR) (KB3207512) - 13.0.4199.0 (X64)   Nov 18 2016 15:56:54   Copyright (c) Microsoft Corporation  Developer Edition (64-bit) on Windows 10 Home 6.3 <X64> (Build 14393: )
   

2. QGIS: Version 2.18.3 'Las Palmas', downloaded from the QGIS Downloads page.

3. ogr2ogr.exe: --version reports GDAL 2.1.3, released 2017/20/01. Downloaded from the GISInternals Downloads page.

1. Computing Area in SQL Server

My starting point is a geometry that I have in a SQL Server database table. Note that I'm specifying 4326 as the SRID (WGS84).

CREATE TABLE dbo.MyTable(
    blockid [varchar](15) NOT NULL,
    geog [geography] NOT NULL
);

INSERT INTO MyTable (blockid, geog)
VALUES ('060730053003', geography::STGeomFromText('POLYGON(( -117.166 32.7199, -117.166 32.7136, -117.165 32.7136, -117.162 32.7126, -117.162 32.7115, -117.161 32.7115, -117.161 32.7157, -117.161 32.7199, -117.166 32.7199 ))', 4326));

I want to compute the area in square kilometers. Based on other examples I've seen online, here's how I can do so in SQL:

SELECT blockid, geog, geog.STArea()/1E6 AS area_squarekm FROM dbo.MyTable WHERE blockid = '060730053003';

This gives me:

╔══════════════╦═════════════════════╦═══════════════════╗
║   blockid    ║        geom         ║   area_squarekm   ║
╠══════════════╬═════════════════════╬═══════════════════╣
║ 060730053003 ║ 0xE6100000010409... ║ 0.364933290578713 ║
╚══════════════╩═════════════════════╩═══════════════════╝

I can also verify the shape looks as expected in SQL Server Management Studio's Spatial Results tab:

Note that for the area, I got 0.364933290578713 km². Let's remember this number. This is the result that I want to compare to QGIS in the steps that follow.

So let's perform a similar computation in QGIS...

2. Computing Area in QGIS for the Same Shape

To start, I'll export the geography from SQL Server into a shapefile that I can pull into QGIS. From the command-line:

ogr2ogr -f "ESRI Shapefile" -nln "060730053003" ".\060730053003" "MSSQL:server=(local);database=mydb;tables=MyTable;trusted_connection=yes;" -sql "SELECT blockid, geog FROM MyTable WHERE blockid='060730053003'"

I check the details of the exported shapefile to make sure things look as expected:

D:\repro\060730053003>ogrinfo -al 060730053003.shp
INFO: Open of `060730053003.shp'
      using driver `ESRI Shapefile' successful.

Layer name: 060730053003
Metadata:
  DBF_DATE_LAST_UPDATE=2017-02-02
Geometry: Polygon
Feature Count: 1
Extent: (-117.166000, 32.711500) - (-117.161000, 32.719900)
Layer SRS WKT:
(unknown)
blockid: String (15.0)
OGRFeature(060730053003):0
  blockid (String) = 060730053003
  POLYGON ((-117.166 32.7199,-117.161 32.7199,-117.161 32.7157,-117.161 32.7115,-117.162 32.7115,-117.162 32.7126,-117.165 32.7136,-117.166 32.7136,-117.166 32.7199))

The points of the polygon appear to be in the reverse order that I specified when placing them in SQL Server, but everything appears intact/as expected.

Now I turn to QGIS and import this shape. I keep WGS84 selected in the initial dialog that comes up:

Then I select the shape and go to Plugins | Python Console. I enter the following series of commands to compute the shape area (I'm using this post as a guide):

wb=iface.activeLayer()
feat = wb.selectedFeatures()
geom = feat[0].geometry().asPolygon()
print geom[0]   # check that the geometry is as expected
area = QgsDistanceArea()
area.setEllipsoid('WGS84') 
area.setEllipsoidalMode(True)
area.computeAreaInit()
area.measurePolygon(geom[0])/1e6

...and the result I get back is 0.36697665124273515 km², as seen in the screenshot below:

3. Summary

Compare the QGIS result, 0.36697665124273515 km², to the one I got from SQL Server before: 0.364933290578713 km². This isn't a massive difference, but it does appear significant, and this is a small shape. I worry whether I'd see differences of larger magnitude for larger shapes that I'm working with. I'm also worried whether it indicates I am doing something wrong in either or both tools to compute area in square meters correctly.

Is there something I'm missing, or does this just reflect inconsistency in calculation and/or floating point error between SQL Server and QGIS?

Answer 1

I have a hard time understanding how a valid area can be returned on spatial data in WGS84 - even the SQL Server documentation confuses me:

Returns the total surface area of a geography instance. Results for STArea() are returned in the square of the unit of measure used by the spatial reference identifier of the geography instance; for example, if the SRID of the instance is 4326, STArea() returns results in square meters.

Any area calculations we make in SQL Server is done on GEOMETRY that is stored in a local projected coordinate system (COLORADO CENTRAL STATE PLANE FEET).

In PostGIS, our data is stored in WGS84, but transformed to State Plane in a query using ST_Transform - which SQL Server does not have the ability to do hence we store our data in projected coordinate system.

I have read that PostGIS Geography will 'pick the best projected coordinate system for the data in the query' ie. it will 'guess' what projected coordinate system to use in order to get valid units.

But - WGS84 is decimal degrees, so back to the main point - is that value you're getting even in M or is it in decimal degrees (even in QGIS?)

Answer 2

Reprojecting your polygon in to EPSG:2771 (based on the Projestions plugin) gives an area of about 0.364941 km2. Any more decimal places is probably spurious depending on the source of the data.