SqlGeometry and Boxes

I wanted to write a blog about partitioning (calling an algorithm using an on-the-fly quadtree or octree) and how it is implemented in Boost.Geometry, and how it could also be used / implemented for the SqlGeometry type. But somehow that blog will be too large and I blog here just about a box/rectangle in SqlGeometry. It is not too difficult, but I just googled and did not find other blogs saying the same things, so I hope it adds something.

STEnvelope

The OGC/ISO STEnvelope function, supported by spatial databases and many geometry libraries, returns the envelope of a geometry. An envelope is: a bounding box or bounding rectangle, also known as an axis aligned bounding box (aabb), a bbox, a minimum bounding rectangle (mbr) or an extent.

More in detail, STEnvelope returns a Geometry object describing a rectangle. OGC does not know the notion of a box. That is sometimes a bit inconvenient, but it is how it is. Boost.Geometry created a Box Concept, because Boost.Geometry is strongly typed, and envelopes returning polygons with always four or five points would confuse the non-OGC users of our library.

I blogged earlier about the spatial extent of a query in SQL Server, here.

A box is a simple geometry, with (in 2D) a minimum x and y and a maximum x and y. But how to get that from a geometry in a spatial database... I mean, the envelope returns a polygon containing five points (assuming it is closed), is the first point lower left? And is that guaranteed? The OGC specifications say: The minimum bounding box for this Geometry, returned as a Geometry. The polygon is defined by the corner points of the bounding box [(MINX, MINY), (MAXX, MINY), (MAXX, MAXY), (MINX, MAXY), (MINX, MINY)]. So yes, it should be guaranteed.

Let's look at it in more detail.

SQL Server

In SQL Server it indeed seems the case that the first point of the polygon returned by STEnvelope is the lower left point. Then, the third point must be the upper right point. Because, if you consider the polygon clockwise, it is the third point, and if it is counterclockwise, it is also the third point.

We use this query:

select geometry::STGeomFromText('POLYGON((0 0,2 5,4 1,0 0))',0).STEnvelope().STPointN(3).STAsText()

... to get the third point (which should be upper right), and it indeed gives me: POINT(4 5)

PostGIS

In PostGIS, the first point of the polygon returned by ST_Envelope is the lower left point as well. Great.

As you (maybe) know, SQLGeometry uses methods, PostGIS uses functions. So the corresponding nested function calls would be:

select ST_AsText(ST_PointN((ST_Envelope(ST_GeomFromText('POLYGON((0 0,2 5,4 1,0 0))',0))), 3))

But... no result, we get a null value back. Looking at the OGC specs, this is correct: NumPoints is defined for a LineString and not for a Polygon. SqlServer is somewhat more relaxed here (as is Boost.Geometry), returning the total number of points of a polygon (including interior rings, if any).

OK, behaviour is correct, so we try again, now inserting the ST_Boundary function (making a linestring of the boundary of a polygon). Our new function call is:

select ST_AsText(ST_PointN(ST_Boundary((ST_Envelope(ST_GeomFromText('POLYGON((0 0,2 5,4 1,0 0))',0)))), 3))

... and we have our upper right point. Using a 1 for the last value (the index of ST_PointN is one-based) would return the lower left point.

Boost.Geometry

Boost.Geometry, as said, returns a box (conforming a Box Concept) for the envelope (or assigns one). So this is more or less outside the discussion... Having a box you can get the minimum-corner coordinates (might be 2D, might be 3D, or more) and the maximum-corner coordinates.

SqlGeometry

I blogged about the useful SqlGeometry type here. The SqlGeometry type is polymorphic, it can be any OGC geometry. But, therefore, it cannot be a Rectangle (because, remember, a rectangle is not an OGC Geometry). So if it contains a Rectangle, it is a Polygon.

Using C#, it would be useful to get the coordinates from a Polygon returned by STEnvelope and that goes like this. Note that (in the light of defensive programming) we don't even assume that the first point is the lower left point. We add this code as a function (we could add it to a class GeometryExtensions, but for now we don't). It is called BoxToCoordinates, see the code below.

And, for symmetry, and because it is convenient, or often necessary (last reason most important ;-) ) we add a counterpart function as well, which creates a rectangle for you given four coordinate values. We call that one CoordinatesToBox, see code below.

The code

OK, this was some preparation for next blog, not too difficult, and now we can convert an STEnvelope result to coordinate values, and back. The full code is displayed below. Between BoxToCoordinates and CoordinatesToBox, there will be some more interesting calculations in the next blog. Stay tuned.

using System;
using Microsoft.SqlServer.Types;

namespace BlogSqlGeometryBox
{
    class Program
    {
        public static SqlGeometry FromWkt(string text, int srid)
        {
           return SqlGeometry.STGeomFromText(new System.Data.SqlTypes.SqlChars(text), srid);
        }

        public static SqlGeometry CoordinatesToBox(double x1, double y1, double x2, double y2, int srid)
        {
           SqlGeometryBuilder builder = new SqlGeometryBuilder();
           builder.SetSrid(srid);
           builder.BeginGeometry(OpenGisGeometryType.Polygon);
           builder.BeginFigure(x1, y1);
           builder.AddLine(x1, y2);
           builder.AddLine(x2, y2);
           builder.AddLine(x2, y1);
           builder.AddLine(x1, y1); // Yes, we close it neatly
           builder.EndFigure();
           builder.EndGeometry();
           return builder.ConstructedGeometry;
        }

        private static void MakeFirstSmaller<T>(ref T a, ref T b) where T : IComparable
        {
           if (a.CompareTo(b) > 0)
           {
               // Exchange the two values
               T t = a;
               a = b;
               b = t;
           }
        }

        public static void BoxToCoordinates(SqlGeometry box, out double x1, out double y1, out double x2, out double y2)
        {
           SqlGeometry lower_left = box.STPointN(1);
           SqlGeometry upper_right = box.STPointN(3);
           x1 = lower_left.STX.Value;
           y1 = lower_left.STY.Value;
           x2 = upper_right.STX.Value;
           y2 = upper_right.STY.Value;

           // Defensive programming:
           MakeFirstSmaller(ref x1, ref x2);
           MakeFirstSmaller(ref y1, ref y2);
        }

        static void Main(string[] args)
        {
           SqlGeometry triangle = FromWkt("POLYGON((0 0,2 5,4 1,0 0))", 0);
           SqlGeometry box = triangle.STEnvelope();

           double x1, y1, x2, y2;
           BoxToCoordinates(box, out x1, out y1, out x2, out y2);

           // Do something with these coordinates. Let's create a larger box
           x1 -= 1;
           y1 -= 1;
           x2 += 1;
           y2 += 1;

           box = CoordinatesToBox(x1, y1, x2, y2, 0);

           System.Console.WriteLine(box.STAsText().ToSqlString().ToString());
        }
    }
}

Writing:

POLYGON ((-1 -1, -1 6, 5 6, 5 -1, -1 -1))

SqlGeometry types and LINQ

Last year I started to use the Microsoft .NET SqlGeometry types and I like them very much. Thanks to Bert Temme who attended us (at Geodan) on their existance. See also this page with an introduction.

In short, the SqlGeometry is the type used for the SQL Server Spatial geometry objects. But it can also be used outside of SQL Server environments: it is contained in a pure .NET assembly. The download page says: "This component can be installed separately from the server to allow client applications to use these types outside of the server."

So the SqlGeometry type is comparable to other Geometry libraries, such as Boost.Geometry, GEOS and NTS (the JTS for .NET).

SqlGeometry functionality follows the OGC/ISO specifications quite closely. This means there is not much more than OGC functionality. Sometimes there is the wish to do other things than that. For example: select only the polygons from a geometry collection. Or: select only the polygons larger than X square meters from a MultiPolygon. There comes LINQ into place.

Installing the assemblies to use SqlGeometry

To use SQLGeometry, you need to have a reference to Microsoft.SqlServer.Types, which can be found in C:\Program Files\Microsoft SQL Server\100\SDK\Assemblies (on my (64 bit) machine in c:\Program Files (x86)\..., indicating 32 bits). If not there, download it from here (there is a 64 bits version there).
Note that you have to add the reference manually, by browsing on the file system, it does not appear in the reference tab (at least not on my system, using Microsoft Visual C# 2010 Express).

After that, just add:

using Microsoft.SqlServer.Types;

The geometry collection

A geometry collection is one geometry consisting of several sub-geometries, which possibly can have different types (e.g. linestrings and polygons).

A geometry collection is sometimes not convenient. For example: MapServer cannot display them. And there are probably more GIS packages which do not display geometry collections, because most mapping software is based on layers, and most layers are defined either on a layer with points (or multi-points), or a layer with lines (or multi-linestrings), or on a layer with polygons (or multi-polygons), but in most cases not on a layer with polygons and lines, and in even more cases not on a layer with geometry collections.

Geometry collections can come into existence during an intersection. If you overlay, for example, two valid polygons which share a part of their boundaries in opposite directions, that shared boundary will form a linestring, while the rest of the intersection (if any) will form a polygon. The polygon and the linestring form together a geometry collection. The next SQL statement generates a (probably unwanted) geometry collection:

select geometry::STGeomFromText('POLYGON((0 0,0 8,8 8,8 4,4 4,4 0,0 0))', 0) 
    .STIntersection(geometry::STGeomFromText('POLYGON((4 0,4 8,8 8,8 0,4 0))', 0))

The result is: GEOMETRYCOLLECTION (POLYGON ((4 4, 8 4, 8 8, 4 8, 4 4)), LINESTRING (4 4, 4 0))

Because geometry collections are often not displayed, functionality to remove all linestrings from a geometry collection would be very useful, or crucial. SqlGeometry does not offer this functionality, but it is possible to build it. It can be built by creating WKT's of all sub-geometries, and parsing and recombining them, but LINQ offers a more elegant way.

Selecting polygons from a geometry collection

We design our algorithm as:

• create a list of geometries from the geometry collection
• remove the linestrings and points from this list
• build a Polygon or a MultiPolygon from the polygons kept in the list

Creating the list

Creating a list of geometries can be done using the STGeometryN and STNumGeometries functions, which are available. We can iterate through the single geometries and add them to the list, using:

var list = new List<SqlGeometry>(); 
for (int i = 1; i <= geometry.STNumGeometries().Value; i++)
{ 
    list.Add(geometry.STGeometryN(i)); 
} 

Note that STGeometryN uses a one-based index, instead of the usual zero-based index.

Because we think this functionality is general, and will often be reused, we make an extension method of it (I really like them!). So we make a static class called GeometryExtensions and we add a static method containing the piece above, and we can just call:

IList<SqlGeometry> list = geom.ToList();

on any variable geom of type SqlGeometry. See source code below for the extension method.

Creating a geometry from a list

Of course, we want to have the reverse as well, making an SqlGeometry from a list of polygons. The key is that an STUnion is required to make a valid result. If two polygons of the list mutually overlap, this overlap will be resolved by STUnion. Even if there is no overlap, the STUnion has to be called because there is no STAdd function available. Alternatively we could make use of the SqlGeometryBuilder type which is available in the SqlGeometry library. However, using STUnion guarantees a valid result. The next loop does this for us:

SqlGeometry result = null; 
foreach (SqlGeometry g in list)
{ 
   result = result == null ? g : result.STUnion(g);
} 

Using LINQ

So having programmed step 1 and step 3 as extension methods, and having a convenient list, we can now implement step 2 using LINQ:

var result = (from g in geom.ToList() where g.STGeometryType().ToString().ToLower() == "polygon" select g).ToSqlGeometry(); 

The magic. This one line statement with LINQ selects the polygons from an on-the-fly list, and put them via another list back into an SqlGeometry. In the listing below it is formatted as three lines.

The complete listing

using System; 
using System.Collections.Generic; 
using System.Linq; 

using Microsoft.SqlServer.Types; 

namespace SqlGeometryBlog1 
{ 
    public static class GeometryExtensions
    { 
        // Converts a multi-geometry to a list of single-geometries 
        public static List<SqlGeometry> ToList(this SqlGeometry geometry)
        { 
          var list = new List<SqlGeometry>(); 
          for (int i = 1; i <= geometry.STNumGeometries().Value; i++)
          { 
              list.Add(geometry.STGeometryN(i)); 
          } 
          return list; 
        } 

        // Converts a list of geometries to a multi-geometry 
        public static SqlGeometry ToSqlGeometry(this IEnumerable<SqlGeometry> list)
        { 
          SqlGeometry result = null; 
          foreach (SqlGeometry g in list)
          { 
              result = result == null ? g : result.STUnion(g);
          } 
          return result; 
        } 

        // For convenience, we add two extra methods 
        public static SqlGeometry FromWkt(string text, int srid)
        { 
          return SqlGeometry.STGeomFromText(new System.Data.SqlTypes.SqlChars(text), srid);
        } 
        public static String ToWkt(this SqlGeometry geometry)
        { 
          return geometry.STAsText().ToSqlString().Value; 
        } 
    } 

    class Program 
    { 
        static void Main(string[] args) 
        { 
          var geom1 = GeometryExtensions.FromWkt("POLYGON((0 0,0 8,8 8,8 4,4 4,4 0,0 0))", 0); 
          var geom2 = GeometryExtensions.FromWkt("POLYGON((4 0,4 8,8 8,8 0,4 0))", 0); 

          var geom = geom1.STIntersection(geom2); 
          System.Console.WriteLine(geom.STIsValid().Value); 
          System.Console.WriteLine(geom.ToWkt()); 

          var list = geom.ToList(); 
          var selection = from g in list where g.STGeometryType().ToString().ToLower() == "polygon" select g; 
          geom = selection.ToSqlGeometry(); 
          System.Console.WriteLine(geom.ToWkt()); 
        } 
    } 
} 

The output

True
GEOMETRYCOLLECTION (POLYGON ((4 4, 8 4, 8 8, 4 8, 4 4)), LINESTRING (4 4, 4 0))
POLYGON ((4 4, 8 4, 8 8, 4 8, 4 4))