A View Inside My Head

I recently developed a spatially-aware .NET application that did not use SQL Server 2008 as the backend (this enterprise was still on SS2005, but we needed the spatial support in the application today).  While the application worked properly on my laptop, it was a huge failboat when deployed to the server environment.

I had previously posted that you can get the Microsoft.SqlServer.Types library from MS Downloads, but it turns out that this alone is not sufficient to allow your application to run.  You also need to ensure that the SQL Server 2008 Native Client is also installed (regardless of whether you're accessing a SS2008 instance or not). Update!You actually don't... read below.

Both the Types library and the Native Client can be downloaded from the following:

http://www.microsoft.com/downloads/details.aspx?FamilyID=228de03f-3b5a-428a-923f-58a033d316e1&DisplayLang=en

My discovery source: https://connect.microsoft.com/SQLServer/feedback/ViewFeedback.aspx?FeedbackID=355402&wa=wsignin1.0

UPDATE: Per Isaac Kunen (in this blog post's comments as well as offline discussion), the missing component from the Types library is simply an updated version of the C Runtime.  The fix of using the Native Client is a hack in this case because its MSI actually installs the updated CRT (which the MSI for the Types library should have done also - it's a goof that MS Downloads hasn't been updated with an updated version of the Types API after the above Connect feedback was answered...

The Microsoft Visual C++ 2008 redistributable by itself can be downloaded from:

http://www.microsoft.com/downloads/details.aspx?FamilyID=A5C84275-3B97-4AB7-A40D-3802B2AF5FC2&displaylang=en

I've previously mentioned that the SQL Server 2008 Spatial data types are freely available for use in your .NET applications, regardless of whether you have SQL Server 2008 or not.  This allows you to incorporate some powerful spatial capabilities right into your application. 

(Look for "Microsoft SQL Server System CLR Types" on this page: http://www.microsoft.com/downloads/details.aspx?FamilyID=228DE03F-3B5A-428A-923F-58A033D316E1&displaylang=en )

However, in most usage scenarios, there will come a time when you have an instance of a SQL Server spatial object in your .NET application, and need to commit it to your SQL Server 2008 database.  How would you do this, without losing fidelity or resorting to serialization of the object to WKT first?

The solutions is to create a Parameter object of type System.Data.SqlDbType.Udt.  Then set the UdtTypeName parameter to the SQL Server-recognized type name (i.e., for SqlGeometry, you would simply use Geometry).

The following code demonstrates executing an UPDATE statement that sets the value of a Spatial field to a newly constructed object.

using (SqlConnection conn = new SqlConnection("Server=.;Integrated Security=true;Initial Catalog=scratch"))
{
   using (SqlCommand cmd = new SqlCommand("UPDATE fe_2007_us_zcta500 SET Boundary=@boundary WHERE id=@id", conn))
   {
      SqlParameter id = cmd.Parameters.Add("@id", System.Data.SqlDbType.Int);

      SqlParameter boundary = cmd.Parameters.Add("@boundary", System.Data.SqlDbType.Udt);
      boundary.UdtTypeName = "geometry";

      SqlGeometry geom = SqlGeometry.Parse("POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))");
      
      boundary.Value = geom;
      id.Value = 123;

      conn.Open();
      cmd.ExecuteNonQuery();
      conn.Close();
   }
}

I have mentioned before how the Ring Orientation for the exterior ring of a Polygon is significant when instantiating a SqlGeography object.  In this case, a Counter-Clockwise orientation is required so that as an observer walks along the path, the interior of the Polygon is always to their left.

But, what I have never really seen documented (or paid attention to, at least) is the fact that the interior rings, or holes, of a Polygon also have specific Ring Orientation requirements. 

In keeping with the "Left-handed" rule, interior rings must be defined in a Clockwise manner - the opposite orientation of the shape's exterior ring.  This is because holes within a Polygon are considered to be part of the exterior of the shape, so the observer walking in a Clockwise direction is still keeping the Polygon's interior to their left.

(I should note here that the Ring Orientation for SqlGeography is the exact opposite of ESRI's ShapeFile format, which is why Ring Orientation has been on my mind for the past few days).

A Ring is a list of points such that the starting point and ending point are the same (forming a closed shape).  The order the you define the points that make up a Ring  - known as Ring Orientation - is significant, for various data formats (including SQL Server's Geography type) imply special meaning for rings that are defined in a clockwise manner as opposed to a counter-clockwise manner. 

Given a list of points with no additional context, it can be difficult to determine the Ring Orientation being used. 

For example, suppose that you have a generic list of points that represent the boundary of a postal code, and that you wish to use these points in order to construct a Polygon instance using the SqlGeography type.  SqlGeography happens to use "Left-handed" ordering, so that as an observer walks along the set of points in the order defined, the "inside" of the polygon is always to their left.  This also implies that the exterior ring of a Polygon is defined in a counter-clockwise manner.

If you try to define a polygon with an area greater than a single hemisphere (this is a nice way to say "if you screw up and use the wrong orientation"), then the SqlGeography type will throw an exception.  So, aside from using Try-Catch, what can you do?

While researching solutions to this problem, I stumbled upon a paper entitled "A Winding Number and Point-in-Polygon Algorithm" from the Colorado State University.  It turns out that a simple algorithm with O(n) complexity can be used to determine if a point is within a Polygon, and a side effect also provides the Ring Orientation.  The key to this algorithm is determining the trend of the ring at each crossing of an axis.

Since I was only interested in Ring Orientation (and not point enclosure detection), I didn't need to use this particular algorithm.  Instead, I took inspiration from the winding concept, and created a simpler derivative algorithm:

1. Iterate the point collection and determine the extreme "left" and "right" points
2. Normalize the line segments connected to these points so that they each have the same "X" dimension length
3. Compare the "Y" values of the normalized segments to establish the trend through that extreme point (i.e., is the "previous" segment above or below the "next" segment)
4. In the spirit of the Winding algorithm, use opposite orientations for the left and right points so that the results coincide with one another
5. A negative result (negative indicates Clockwise orientation, positive result indicates Counter-Clockwise orientation, and a result of zero would be undefined

I've actually written (and posted) several versions of this algorithm, each time discovering some edge case exception that would cause me to take down the post and rewrite the algorithm.  I believe the code below works for all simple polygons on a Cartesian coordinate system (read: I have more testing to see if this will work with an ellipsoidal model, like SqlGeography). 

Note: The following code is generic in nature, and as such, I've defined my own Point structure instead of using a SqlGeometry or SqlGeography, etc.

struct Point
{
    public double X { get; set; }
    public double Y { get; set; }
}

enum RingOrientation : int
{
    Unknown = 0,
    Clockwise = -1,
    CounterClockwise = 1
};

RingOrientation Orientation(Point[] points)
{
    // Inspired by http://www.engr.colostate.edu/~dga/dga/papers/point_in_polygon.pdf

    // This algorithm is to simply determine the Ring Orientation, so to do so, find the
    // extreme left and right points, and then check orientation

    if (points.Length < 4)
    {
        throw new ArgumentException("A polygon requires at least 4 points.");
    }

    if (points[0].X != points[points.Length - 1].X || points[0].Y != points[points.Length - 1].Y)
    {
        throw new ArgumentException("The array of points is not a polygon.  The first and last point must be identical.");
    }

    int rightmostIndex = 0;
    int leftmostIndex = 0;

    for (int i = 1; i < points.Length; i++)
    {
        if (points[i].X < points[leftmostIndex].X)
        {
            leftmostIndex = i;
        }
        if (points[i].X > points[rightmostIndex].X)
        {
            rightmostIndex = i;
        }
    }


    Point p0; // Point before the extreme
    Point p1; // The extreme point
    Point p2; // Point after the extreme

    double m; // Holds line slope

    double lenP2x;  // Length of the P1-P2 line segment's delta X
    double newP0y;  // The Y value of the P1-P0 line segment adjusted for X=lenP2x

    RingOrientation left_orientation;
    RingOrientation right_orientation;

    // Determine the orientation at the Left Point
    if (leftmostIndex == 0)
        p0 = points[points.Length - 2];
    else
        p0 = points[leftmostIndex - 1];

    p1 = points[leftmostIndex];

    if (leftmostIndex == points.Length - 1)
        p2 = points[1];
    else
        p2 = points[leftmostIndex + 1];

    m = (p1.Y - p0.Y) / (p1.X - p0.X);

    if (double.IsInfinity(m))
    {
        // This is a vertical line segment, so just calculate the dY to
        // determine orientation

        left_orientation = (RingOrientation)Math.Sign(p0.Y - p1.Y);
    }
    else
    {
        lenP2x = p2.X - p1.X;
        newP0y = p1.Y + (m * lenP2x);

        left_orientation = (RingOrientation)Math.Sign(newP0y - p2.Y);
    }



    // Determine the orientation at the Right Point
    if (rightmostIndex == 0)
        p0 = points[points.Length - 2];
    else
        p0 = points[rightmostIndex - 1];

    p1 = points[rightmostIndex];

    if (rightmostIndex == points.Length - 1)
        p2 = points[1];
    else
        p2 = points[rightmostIndex + 1];

    m = (p1.Y - p0.Y) / (p1.X - p0.X);

    if (double.IsInfinity(m))
    {
        // This is a vertical line segment, so just calculate the dY to
        // determine orientation

        right_orientation = (RingOrientation)Math.Sign(p1.Y - p0.Y);
    }
    else
    {
        lenP2x = p2.X - p1.X;
        newP0y = p1.Y + (m * lenP2x);

        right_orientation = (RingOrientation)Math.Sign(p2.Y - newP0y);
    }


    if (left_orientation == RingOrientation.Unknown)
    {
        return right_orientation;
    }
    else
    {
        return left_orientation;
    }
}


void Test()
{
    // Simple triangle - left extreme point is vertically "in between" line segments
    Point[] points = new Point[]
    {
        new Point(5,-1),
        new Point(0,0),
        new Point(5,1),
        new Point(5,-1)
    };

    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.Clockwise);
    Array.Reverse(points);
    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.CounterClockwise);

    // Case where both line segments are above the left extreme point
    points = new Point[]
    {
        new Point(2,1),
        new Point(0,0),
        new Point(1,1),
        new Point(2,1)
    };

    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.Clockwise);
    Array.Reverse(points);
    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.CounterClockwise);

    // Case where both line segments are below the left extreme point
    points = new Point[]
    {
        new Point(2,-1),
        new Point(0,0),
        new Point(1,-1),
        new Point(2,-1)
    };

    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.CounterClockwise);
    Array.Reverse(points);
    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.Clockwise);

    // Case where line segment is vertical (slope cannot be determined)
    points = new Point[]
    {
        new Point(0,0),
        new Point(0,1),
        new Point(1,1),
        new Point(1,0),
        new Point(0,0)
    };

    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.Clockwise);
    Array.Reverse(points);
    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.CounterClockwise);

    // Case where angle thru left extreme point is a right angle
    points = new Point[]
    {
        new Point(0,0),
        new Point(1,1),
        new Point(1,-1),
        new Point(0,0)
    };

    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.Clockwise);
    Array.Reverse(points);
    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.CounterClockwise);

    // Real-world case from a SHP file
    points = new Point[]
    {
        new Point(-156.92467299999998,20.738695999999997),
        new Point(-156.924636,20.738822),
        new Point(-156.924608,20.73894),
        new Point(-156.92458,20.739082),
        new Point(-156.92460599999998,20.739234),
        new Point(-156.924551,20.739326),
        new Point(-156.924507,20.739241999999997),
        new Point(-156.924482,20.739082),
        new Point(-156.924466,20.738854999999997),
        new Point(-156.924387,20.738602999999998),
        new Point(-156.924308,20.738325),
        new Point(-156.924239,20.738063999999998),
        new Point(-156.92424,20.737887999999998),
        new Point(-156.924285,20.737811999999998),
        new Point(-156.924475,20.73762),
        new Point(-156.92458299999998,20.737603999999997),
        new Point(-156.924754,20.737579),
        new Point(-156.924851,20.737731),
        new Point(-156.924956,20.738101),
        new Point(-156.924909,20.738343999999998),
        new Point(-156.924818,20.738487),
        new Point(-156.92467299999998,20.738695999999997)
    };

    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.Clockwise);
    Array.Reverse(points);
    System.Diagnostics.Debug.Assert(Orientation(points) == RingOrientation.CounterClockwise);

}

In this, the eighth part in a series on the new Spatial Data types in SQL Server 2008, I'll step away from the database and do a little spatial coding using .NET.

Redistributable .NET Library

Up to this point in the series, I have demonstrated a lot of interesting (?) things that you can do with the new Spatial data types (Geometry and Geography) in SQL Server 2008.  You might be thinking, "That's swell, and all, but I wish I could do some of that stuff without needing to be tethered to a database."  Well, you know what?  You can!

I mentioned in a previous post that the Spatial data types were implemented as SQLCLR User-Defined Types.  I've since been corrected by Isaac Kunen, who stated that they are more accurately described as System-Defined Types, with the difference being that these are automatically installed and available for use as part of SQL Server 2008, regardless of whether the ENABLE CLR bit has been activated.  Semantics aside, these types are merely classes within a .NET assembly, and Microsoft is making this freely available as part of a Feature Pack for SQL Server (which will be redistributable as part of your stand-alone application, according to Isaac):

SQL Server 2008 RC0 Feature Pack Download Page

(Look for "Microsoft SQL Server System CLR Types," which includes the two Spatial types plus the HierarchyID type.  This link is for RC0, and may not be applicable to future versions as the product is finalized and released.)

Builder API

A new feature that was included with the first Release Candidate (RC0) is the Builder API.  This is a collection of interfaces and classes that helps you to construct spatial types by specifying one point at a time until all points have been added.

The Builder API is not only useful for creating new instances of spatial data, but also for consuming existing instances one point at a time (maybe to convert an instance into another format).  Documentation is light at the moment, so I'm still trying to grok exactly how to best utilize it.

For my first experiment with the API, I obtained some Zip Code Boundary data in ASCII format from the U.S. Census Bureau:

http://www.census.gov/geo/www/cob/z52000.html#ascii

My goal was to parse the data, and then create a new SqlGeography instance for each zip code.  (Note: SqlGeography is the .NET class name that T-SQL refers to simply as Geography).  The SqlGeographyBuilder class proved to be perfect for accomplishing this task.

At its core, the SqlGeographyBuilder implements the IGeographySink interface.  If you wanted to consume an existing SqlGeography instance, you could implement IGeographySink in your own class, and then invoke the SqlGeography's Populate() instance method, passing in your object as the parameter.  The Populate() method takes care of calling the appropriate IGeographySink methods within your class.

In this case, I'm not starting with an existing SqlGeography instance, so my code will need to call the methods of the SqlGeographyBuilder in the correct order:

After EndGeography() has been invoked, the new instance is available via the ConstructedGeography property of the SqlGeographyBuilder class. 

Simple enough, right?  Yeah, I'm still a little lost myself...  But, here's some code to help demonstrate what's going on!

First, let's look at the ASCII data.  A single zip code's boundary might be defined as:

      1469      -0.824662148292608E+02       0.413848583827499E+02
      -0.824602851767940E+02       0.413864290595145E+02
      -0.824610630000000E+02       0.413860590000000E+02
      -0.824685900000000E+02       0.413841470000000E+02
      -0.824686034536111E+02       0.413843846804627E+02
      -0.824605990000000E+02       0.413863160000000E+02
      -0.824602851767940E+02       0.413864290595145E+02
END

The very first line happens to contain an identifier (maps to a second file that lists the actual USPS zip code).  The coordinate listed in the first line is not actually part of the boundary, but rather appears to be the population center of that area.  The actual boundary begins with the second line, and continues until you encounter the "END".  Also, in case you couldn't tell, coordinates in this data are in Longitude-Latitude order.

Since a Zip Code is a polygon, and since we are working with SqlGeography, we must be aware of ring ordering.  That is, the exterior ring of a polygon must be defined in a counter-clockwise order so that as you "walk the ring", the interior is always to your left.  If you reverse the order, then SqlGeography assumes that you're trying to define a polygon containing the entire world except for the small area inside of the polygon.

Well, in this case, the order of the points of the Zip Code boundary is defined in clockwise order... so, we must be aware of this and call into the SqlGeographyBuilder in the opposite order (so the last point defined in the ASCII data is the first point used while building our new instance). 

To accomplish this, I simply parse the Lat/Long coordinates as "double" types, and then push them onto a stack.  Then, I pop the stack and call into the Builder API with each point.  At the end, I obtain the new SqlGeography instance from the ConstructedGeography property. 

(Note: This is demonstrative code - some things should probably be cleaned up/refactored/error handled... You have been warned)

  public SqlGeography ParseAsGeography(string zipcode_points)
{
    StringReader sr = new StringReader(zipcode_points);
    string line = sr.ReadLine();

    Stack<double[]> Points = new Stack<double[]>();

    while (line != null  && line != "END")
    {
        if (line != String.Empty)
        {
            Points.Push(ParseLatLngValues(line));
        }

        line = sr.ReadLine();
    }

    return CreateGeography(Points);
}

  private
  double[] ParseLatLngValues(string line)
{
    //      -0.838170700000000E+02       0.409367390000000E+02double[] ret = newdouble[2];

    string lng = System.Text.RegularExpressions.Regex
                            .Matches(line, "\\S+")[0].Value;
    string lat = System.Text.RegularExpressions.Regex
                            .Matches(line, "\\S+")[1].Value;

    double.TryParse(lat, out ret[0]);
    double.TryParse(lng, out ret[1]);

    return ret;
}

  private SqlGeography CreateGeography(Stack<double[]> points)
{
    SqlGeographyBuilder builder = new SqlGeographyBuilder();
    builder.SetSrid(4326);
    builder.BeginGeography(OpenGisGeographyType.Polygon);

    double[] point = points.Pop();

    builder.BeginFigure(point[0], point[1]);

    while (points.Count > 0)
    {
        point = points.Pop();
        builder.AddLine(point[0], point[1]);
    }

    builder.EndFigure();
    builder.EndGeography();

    return builder.ConstructedGeography;
}

SQL Server 2008: Spatial Data, Part 1

SQL Server 2008: Spatial Data, Part 2

SQL Server 2008: Spatial Data, Part 3

SQL Server 2008: Spatial Data, Part 4

SQL Server 2008: Spatial Data, Part 5

SQL Server 2008: Spatial Data, Part 6

SQL Server 2008: Spatial Data, Part 7