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:
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: )
QGIS: Version 2.18.3 'Las Palmas', downloaded from the QGIS Downloads page.
ogr2ogr.exe
:--version
reportsGDAL 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 km2. 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 km2, as seen in the screenshot below:
3. Summary
Compare the QGIS result, 0.36697665124273515 km2, to the one I got from SQL Server before: 0.364933290578713 km2. 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?