Showing posts with label Extent. Show all posts
Showing posts with label Extent. Show all posts

Friday, May 6, 2011

SqlGeometry and Boxes


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

Friday, April 8, 2011

Extent of a SQL Server Spatial table


Extent of a SQL Server Spatial table


PostGIS does have a convenient spatial aggregation function: ST_Extent. SQL Server Spatial does not have that functionality. This little blog shows how to achieve its effect.

PostGIS


Given a table (e.g. called world1) with a geometry column (e.g. called geom), a call to the PostGIS function ST_Extent will give a neat box containing all geometries. For example:

select ST_Extent(geom) from world1;

returns a bounding box:

BOX(-180 -89.9,180 83.674733)

SQL Server Spatial


SQL Server Spatial follows OGC closely and therefore does not implement this functionality. Neither it has implemented it as a so-called extended geometry method. I've seen various creative solutions on the web, e.g. here and here and here. Most of these efforts are based on iterating through the rows, which is possible but often inconvenient (because of the required loop) and (also therefore) probably less performant. I will not say that my solution is performing better, but at least it is one SQL statement.

Be prepared on use of with, my current favorite SQL clause.  I described it earlier here.

So what happens: we define a sort of on-the-fly view of all the envelopes (select geom.MakeValid().STEnvelope() as envelope from world1). (The call to MakeValid() is only necessary if there might be invalid geometries in the table). Then we ask for the corner points of that viewy thingy. An envelope is a rectangle, always containing five points, where the fifth point is closing point and therefore the same as the first point. So either point 1 and 3, or point 2 and 4 are the opposite corners of a rectangle. So we take point 1 and point 3 (as done in first creative solution I referenced above). We define a second on-the-fly view of this, using the first view as input. To get both diagonally opposite points, we union them (union all is required). Note that I'm referring to an SQL union here, not to a spatial union.

Those on-the-fly views can be (very conveniently) implemented using with. I name them then viewy but that's up to you. So envelope_viewy is the first one and corner_viewy is the second one.

The complete SQL statement is then, for a column called geom in a table called world:

with
  envelope_viewy as
  (
    select geom.MakeValid().STEnvelope() as envelope from world1
  ),
  corner_viewy as
  (
    select envelope.STPointN(1) as point from envelope_viewy
    union all select envelope.STPointN(3) from envelope_viewy
  )
select min(point.STX) as min_x,min(point.STY) as min_y,max(point.STX) as max_x,max(point.STY) as max_y
from corner_viewy

And voilĂ , you will get:
min_x    min_y    max_x    max_y
-180    -89.9    180    83.674733

The same! (Of course, I used the same table, world1, used from my shapefile research-still-to-be-finished).

You can also create a stored procedure of this (e.g. called ST_Extent).

PostGIS remark


Writing this blog, I noticed a strange thing of the PostGIS ST_Extent. It only works for 2D, there is explicitly written in the documentation that it does not work for 3D (use a ST_Extent3D for that). But it delivers a 3D box. And indeed it is explicitly typed as a 3D box. I wonder why.