[GIS] Understanding Differences in Area Calculations between SQL Server and QGIS

areaqgissql server

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 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?

Best Answer

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?)

Related Solutions

[GIS] How to force sql server to use the spatial index

The issue is with your use of STDistance rather than the index. This function calculates the distance between two points so you need to replace HomeGeog with a feature rather than the name of the table.

E.g. @g.STDistance(@g) should work (but always return 0)

What are you trying to return from your query?

See http://msdn.microsoft.com/en-us/library/bb933952.aspx for more details on how to use the STDistance function.

[GIS] SQL Server to Google Maps

Here is how I did it using MVC and VB.NET

This SchoolDistrictVertices function is in a partial class of a SchoolDistrict object which was generation using EntityFramework from a SQL database. It contains a geography field named DistrictLocation.

Public Structure Coordinate
  Public Property Latitude As Double
  Public Property Longitude As Double
End Structure

Public Function SchoolDistrictVertices() As List(Of List(Of Coordinate))
  If SchoolDistrict Is Nothing OrElse SchoolDistrict.DistrictLocation Is Nothing Then
    Return Nothing
  Else
    'get the points that define the shape by parsing the shape's WKT
    Dim result = New List(Of List(Of Coordinate))
    Dim wkt = SchoolDistrict.DistrictLocation.Location.WellKnownValue().WellKnownText
    Dim startPos = wkt.IndexOf("((")

    Do Until startPos = -1
      Dim list = New List(Of Coordinate)
      Dim endPos = wkt.IndexOf(")", startPos)
      Dim piece = wkt.Substring(startPos + 1, endPos - (startPos + 1))
      Dim coords = piece.Split(","c)

      For Each item In coords
        item = item.Trim().Replace("(", "")
        Dim parts = item.Split(" "c)
        Dim lng = parts(0)
        Dim lat = parts(1)
        Dim coord = New Coordinate()
        coord.Latitude = Convert.ToDouble(lat)
        coord.Longitude = Convert.ToDouble(lng)

        list.Add(coord)
      Next
      result.Add(list)

      startPos = wkt.IndexOf("(", endPos)
    Loop

    Return result
  End If

End Function

Here is the View file which takes a SchoolDistrict object as the model:

@<div style="margin-top: 5px; height: 500px" id="map-canvas">
</div>

@Section Scripts
<script src="https://maps.googleapis.com/maps/api/js?v=3.exp&signed_in=false"></script>

<script>
  $(function () {
    google.maps.event.addDomListener(window, 'load', InitializeMap());
  });

  function InitializeMap()
  {
    var mapOptions = {
      zoom: 12
    };
    var map = new google.maps.Map(document.getElementById('map-canvas'),  mapOptions);


    @code
      Dim count = 0
      Dim shapeList = ""
    End Code
    @If Model.SchoolDistrict IsNot Nothing Then
  @For Each shape In Model.SchoolDistrictVertices()
  @:var coordList@(count) = [
    @For Each coordinate In shape
    @:new google.maps.LatLng(@coordinate.Latitude, @coordinate.Longitude),
    Next

                     If shapeList <> "" Then
                       shapeList += ", "
                     End If
                     shapeList += "coordList" + count.ToString()
                     count += 1
    @:];
         Next

    @:districtShape = new google.maps.Polygon({paths: [@shapeList], strokeColor: '#A366E0', strokeOpacity: 0.8, strokeWeight: 1, fillColor:  '#A366E0', fillOpacity: 0.35});
                    @:districtShape.setMap(map);
                         End If
  }

</script>
End Section