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:
First, let me state what my environment is:
SQL Server: I'm running SQL Server 2016 Developer Edition.
SELECT @@VERSIONgives 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.
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 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.geometry().asPolygon() print geom # check that the geometry is as expected area = QgsDistanceArea() area.setEllipsoid('WGS84') area.setEllipsoidalMode(True) area.computeAreaInit() area.measurePolygon(geom)/1e6
...and the result I get back is 0.36697665124273515 km2, as seen in the screenshot below:
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?