locked
Iterative Query for Spatial Recordrs RRS feed

  • Question

  • Hello

    I have a xml column in SQL Server 2012 table in which I have xml files that consists of lat/long of a polygone. I have each separate row for each new polygone. The xml format is as following

    <Dataset_Extent>
          <EXTENT_TYPE>Bounding-Polygon</EXTENT_TYPE>
          <Vertex>
            <LON>56.23654/LON>
            <LAT>30.27001</LAT>
            <X>299232</X>
            <Y>3350549</Y>
            <COL>1</COL>
            <ROW>1</ROW>
          </Vertex>
          <Vertex>
            <LON>76.99454</LON>
            <LAT>30.27128</LAT>
            <X>307089.5</X>
            <Y>3350549</Y>
            <COL>15715</COL>
            <ROW>1</ROW>
          </Vertex>
          <Vertex>
            <LON>36.99700</LON>
            <LAT>30.15096</LAT>
            <X>307089.5</X>
            <Y>3337207.5</Y>
            <COL>15715</COL>
            <ROW>26683</ROW>
          </Vertex>
          <Vertex>
            <LON>86.91543</LON>
            <LAT>30.14969</LAT>
            <X>299232</X>
            <Y>3337207.5</Y>
            <COL>1</COL>
            <ROW>26683</ROW>
          </Vertex>
          <Center>
            <LON>56.95499</LON>
            <LAT>30.21048</LAT>
            <X>303160.75</X>
            <Y>3343878.25</Y>
            <COL>7858</COL>
            <ROW>13342</ROW>
          </Center>
        </Dataset_Extent>

    Now I have polygon of an area whoes lat/long I am passing to sql from front end application to check which polygon inside my database intersect with this by using STIntersect function. 

    I have tested STIntersect method on a single record like

    DECLARE @g geography;
    DECLARE @h geography;
    SET @g = geography::STGeomFromText('POLYGON((-122.358 47.653, -122.348 47.649, -122.348 47.658, -122.358 47.658, -122.358 47.653))', 4326);
    SET @h = geography::STGeomFromText('POLYGON(-122.360 47.656, -122.343 47.656,-122.360 47.656)', 4326);

    SELECT CASE @g.STIntersects(@h)

    WHEN 1 THEN '@g intersects @h'

    ELSE '@g does not intersect @h'

    END;

    How can I write a SQL query that check each record against my provided polygon and returns results if both intersect.

    Thanks


    Engr. Mudassar Ali Software Engineer

    Tuesday, September 2, 2014 5:15 PM

Answers

  • Hi Mudassar,

    There's a subtle bug in the code. Make the change below and it works. BTW, both image1 and image2 did need to be reoriented, which is likely why you had more "false positives" in STIntersects using the cursor originally.

    DECLARE @image1 geography
    DECLARE @image2 geography

    SET @image1=geography::STGeomFromText('POLYGON((66.91292909247741 30.27001012181008,66.99456841960638 30.27128639618252,66.99700791329992 30.1509623521339,66.91546772378466 30.14969219541345,66.91292909247741 30.27001012181008))',4326);
    SET @image2=geography::STGeomFromText('POLYGON((74.22257407407407 32.21185185185185,74.46148148148148 32.21185185185185,74.46148148148148 32.04995833333334,74.22257407407407 32.04995833333334,74.22257407407407 32.21185185185185))',4326);

    Declare @fg geography=geography::STGeomFromText('FULLGLOBE',4326);
    Declare @a float=@fg.STArea()/2;
    declare @myGeo1 geography=@image1
    declare @myGeo2 geography=@image2
    if (@myGeo1.STArea()>@a)
        set @image1 =@myGeo1.ReorientObject();
    if (@myGeo2.STArea()>@a)
        set @image2 =@myGeo1.ReorientObject();   -- should be @myGeo2.ReorientObject()

    SELECT CASE @image1.STIntersects(@image2)
    WHEN 1 THEN 'cutting'
    ELSE 'Not cutting'
    END;

    Cheers, Bob

    • Marked as answer by Sofiya Li Wednesday, September 10, 2014 9:05 AM
    Tuesday, September 9, 2014 4:16 PM
  • Sure. You can use STIntersection to get a polygon that represents the intersection of the two polygons, then use STArea and compare the intersection area vs. total area.

    Cheers, Bob

    • Marked as answer by Mudassar Alii Thursday, September 18, 2014 10:05 AM
    Thursday, September 18, 2014 12:10 AM
  • Glad that it worked out. That typo took me a while to find too. ;-)

    Which data type you want to use depends on what you're doing with it. If you're measuring, for example, the distance between Seattle and Boston, the geography type takes the shape of the earth into consideration, geometry doesn't. But for some things, like drawing shapes on a map or getting the distances between two locations that are near each other, geometry is fine. In fact, geometry is faster to do calculations, because it doesn't consider the earth's shape the math is easier and faster.

    I'm not a GUI specialist with spatial. I work almost exclusively with the database, specializing in getting best performance with big spatial tables and using spatial indexes. And the spatial types in SQL Server don't always map one-to-one with the classes in .NET's System.Drawing. If you're using real places on a map you'll be better off imbedding the Bing map control in a C# program. or using Google or Bing map APIs directly in a web page. You might find some good ideas in Ricky's Bing Map Blog (http://rbrundritt.wordpress.com/) or have a look at the code in Alastair Aitchison's (http://alastaira.wordpress.com/) book on SQL Server spatial, I believe he was Bing and Google API examples.

    Finally, there's an analogous forum for Bing map control programming here: http://social.msdn.microsoft.com/Forums/onedrive/en-US/home?forum=vemapcontroldev 

    Hope this helps get you started, Cheers, Bob




    • Edited by Bob Beauchemin Wednesday, September 10, 2014 5:04 AM
    • Marked as answer by Sofiya Li Wednesday, September 10, 2014 9:05 AM
    Wednesday, September 10, 2014 4:57 AM

All replies

  • You could do this by writing a function named MakePolygon that takes the XML data type as input and produces the geography data type as output, using a combination of XQuery FLWOR expression, string concatenation (to WKT), and the STGeomFromText static method. Or, instead of doing that in T-SQL, it might be faster in SQLCLR using either LINQ to XML or System.Xml.

    Watch out for a few things. First, you need to close the polygon (i.e. the first and last vertex in your polygon need to be equal). And don't forget that, when using WKT it needs (Long, Lat) vertices, not (Lat,Long).

    At that point, your T-SQL query becomes:

    SELECT ....
    FROM yourtable
    WHERE dbo.MakePolygon(your_xml_col).STIntersects(@passed_in_polygon)=1;

    If your table is fairly large, consider making a persisted computed column of type geography from MakePolygon and putting it in your table. Then you can have a spatial index to improve performance, as well as ameliorate the XQuery performance hit. 

    If your XML actually looks like what you've posting, this might help to get you started. It's ugly, unoptimized, and returned a pretty funny-looking polygon with your exemplar document, however.

    create function dbo.MakePolygon(@x xml)
    --returns varchar(8000) - for testing
    returns geography
    as
    begin
    declare @wkt varchar(8000);
    select @wkt = CONVERT(varchar(7000), @x.query('
     for $lon in /Dataset_Extent/Vertex/LON/text()
     for $lat in /Dataset_Extent/Vertex/LAT/text()
     return ($lon cast as xs:string?, $lat cast as xs:string?, ",")
    '));

    -- concatenate the first vertex
    select @wkt = @wkt + CONVERT(varchar(1000), @x.query('
     for $lon in (/Dataset_Extent/Vertex/LON/text())[1]
     for $lat in (/Dataset_Extent/Vertex/LAT/text())[1]
     return ($lon cast as xs:string?, $lat cast as xs:string?)
    '));

    -- concatenate the word POLYGON and required parens
    set @wkt = 'POLYGON((' + @wkt + '))';

    return geography::STGeomFromText(@wkt,4326);
    --return @wkt (for testing)
    end

    Hope this helps, Bob

     
    • Edited by Bob Beauchemin Tuesday, September 2, 2014 9:47 PM added function example
    • Proposed as answer by Sofiya Li Wednesday, September 3, 2014 9:21 AM
    • Unproposed as answer by Mudassar Alii Thursday, September 4, 2014 4:51 AM
    Tuesday, September 2, 2014 8:43 PM
  • Thanks  Bob Beauchemin for the good reply

    The vertices that I have mentioned in my above xml are just dummy one for just testing.

    I have created following query by understanding your post

    declare @x XML
    set @x='<Dataset_Extent>
          <EXTENT_TYPE>Bounding_Polygon</EXTENT_TYPE>
          <Vertex>
            <LON>66.91292909247741</LON>
            <LAT>30.27001012181008</LAT>
            <X>299232</X>
            <Y>3350549</Y>
            <COL>1</COL>
            <ROW>1</ROW>
          </Vertex>
          <Vertex>
            <LON>66.99456841960638</LON>
            <LAT>30.27128639618252</LAT>
            <X>307089.5</X>
            <Y>3350549</Y>
            <COL>15715</COL>
            <ROW>1</ROW>
          </Vertex>
          <Vertex>
            <LON>66.99700791329992</LON>
            <LAT>30.1509623521339</LAT>
            <X>307089.5</X>
            <Y>3337207.5</Y>
            <COL>15715</COL>
            <ROW>26683</ROW>
          </Vertex>
          <Vertex>
            <LON>66.91546772378466</LON>
            <LAT>30.14969219541345</LAT>
            <X>299232</X>
            <Y>3337207.5</Y>
            <COL>1</COL>
            <ROW>26683</ROW>
          </Vertex>
          <Center>
            <LON>66.9549932872921</LON>
            <LAT>30.21048776638499</LAT>
            <X>303160.75</X>
            <Y>3343878.25</Y>
            <COL>7858</COL>
            <ROW>13342</ROW>
          </Center>
        </Dataset_Extent>';
    declare @wkt varchar(8000);
    select @wkt=CONVERT(varchar(7000),@x.query('distinct-values(
    for $lon in /Dataset_Extent/Vertex/LON/text()
    for $lat in /Dataset_Extent/Vertex/LAT/text()
    return ($lon cast as xs:string?,$lat cast as xs:string?,","))
    '));

    --concatenate the word POLYGON
    set @wkt='POLYGON(('+@wkt +'))';
    print @wkt

    This give me the following output

    POLYGON((66.91292909247741 66.99456841960638 30.27128639618252 66.99700791329992 66.91546772378466 30.27001012181008,30.1509623521339 30.14969219541345))

    Each Lon/lat pair is not in proper order, i required output in following format

    POLYGON((66.91292909247741 30.27001012181008,66.99456841960638 30.27128639618252,66.99700791329992 30.1509623521339,66.91546772378466 30.14969219541345,66.91292909247741 30.27001012181008))

    The starting vertex Lon/Lat should also repeat at end to make polygon close.

    Thanks


    Engr. Mudassar Ali Software Engineer



    Wednesday, September 3, 2014 9:41 AM
  • Any idea?

    Engr. Mudassar Ali Software Engineer

    Thursday, September 4, 2014 5:45 PM
  • It's the addition of the XQuery distinct-values function to the original example that makes the code produce the answer that you don't want.  Because XQuery with distinct-values is omitting any values where the Lon or Lat coordinate is the same as any other Lon or Lat coordinate, because XQuery is really producing a sequence of text nodes and doesn't distinguish Lon/Lat pairs with respect to distinctness. 

    You don't need that function. If the coordinates of a polygon repeat (except for the first and last coordinates which must be the same), it's an invalid polygon because, by definition, a polygon can't overlap itself.  If you get an invalid polygon it won't work with spatial functions, but you can fix it with something like:

    if (@my_polygon.STIsValid() = 0) set @my_polygon = @my_polygon.MakeValid();

    After this you may not have a single polygon, but the resulting spatial instance will work correctly with STIntersects. If you need to include that additional logic, that also why it might be useful to encapsulate your XML to spatial type extraction to a user-defined function.

    Cheers, Bob

    Saturday, September 6, 2014 5:31 PM
  • ok thanks. Now I have done with the making polygon by getting vertices (lon/lat) from table. Now I have written a query to check which polygon in my table intersects with one passed as parameter

    Declare @aoi geography 
    Declare @xml xml 
    Declare @result geography
    
    SET @aoi=geography::STGeomFromText('POLYGON((66.9214310204
    30.1191654521,66.9724466643 30.1273507980,66.9612156084 30.1614246796,66.9324717194 30.1540007613,66.9214310204 30.1191654521))',4326);
    
    --creating cursor to hold results for function calling 
    Declare xmlDataCursor CURSOR for select XmlFiles.xmlData from XmlFiles  
    OPEN xmlDataCursor 
    FETCH NEXT FROM xmlDataCursor INTO @xml 
    WHILE @@FETCH_STATUS=0 
    BEGIN 
        Select @xml Where(@aoi.STIntersects(dbo.MakePolygon(@xml))=1)   
        FETCH NEXT FROM xmlDataCursor INTO @xml
    END 
    CLOSE xmlDataCursor 
    DEALLOCATE xmlDataCursor'

    but this will show me the whole results of the table and seems that its not performing WHERE condition part


    Engr. Mudassar Ali Software Engineer


    Saturday, September 6, 2014 6:08 PM
  • I don't think you need a cursor for this. This query will do the same thing without a cursor:

    -- declare and set @aoi, as in your example
    -- this query should get you all the rows that intersect
    select XmlFiles.xmlData
    from XmlFiles
    where @aoi.STIntersects(dbo.MakePolygon(XmlFiles.xmlData))=1)

    If you're still getting all the rows, you need to check a few rows "by hand". Choose a few rows that you think shouldn't intersect. Here's a few things you can check:

    1. Visualize the polygon created by MakePolygon, "UNION ALL" with your aoi polygon. Check the spatial results tab and see if they really do intersect.

    2. Area of the polygon. If it's more than half the area of the earth, you know the polygon vertex order is "backward" for geography. You can get the area of the earth for SRID 4326 by making a FULLGLOBE geography instance and finding the area  of that. See http://www.sqlskills.com/blogs/bobb/measuring-the-earth-with-sql-server-denali/ . If a polygon is "backward" (more than half the area of the earth) you can fix it by using the ReorientObject() method of the geography instance.

    Hope this helps, Bob

    Sunday, September 7, 2014 3:10 AM
  • Thanks Bob for the reply. I have checked my query again and when I used geometry data type instead of geography then its giving me the correct results. Whats the matter is?

    Engr. Mudassar Ali Software Engineer

    Monday, September 8, 2014 5:38 PM
  • I think it's a ring orientation issue with some of the polygons, which isn't a problem with geometry but is with geography. See the part in my last message about "checking for area" (#2) with geography. That's my best guess, without looking at an example.

    Here's a description of how the ring orientation works with geography. http://blogs.msdn.com/b/edkatibah/archive/2008/08/19/working-with-invalid-data-and-the-sql-server-2008-geography-data-type-part-1b.aspx In SQL Server 2008 the ring orientation issue would cause an error (described in the article) because geography instances weren't allowed to be bigger than a hemisphere. In SQL Server 2012, they are allowed to be bigger than a hemisphere, so they don't cause an error but just make very big polygons (that could STIntersect with almost anything). 

    Geometry type doesn't care about ring orientation. Since most databases/data sources don't support a geography data type, these can pop up in shapefiles or any geographic data that's not from SQL Server.

    The way that I've dealt with this (you could put it in MakePolygon and also when you import data to compare to) is to use something like this:

    -- This creates a fullglobe geography with your SRID
    DECLARE @fg geography = geography::STGeomFromText('FULLGLOBE',4326);
    -- Half the area of the fullglobe
    DECLARE @a float = @fg.STArea()/2;
    --- then later
    IF @your_geog.STArea() > @a -- if its greater than hemiphere
        SET @your_geog = @your_geog.ReorientObject(); -- switch the ring orientation

    Unless you're expecting instances that cover more than half the earth. ;-)

    If that's not it, post an example where STIntersect returns false with geometry types and true with same data and geography types.

    Hope this helps, Bob

    • Proposed as answer by Bob Beauchemin Wednesday, September 10, 2014 4:58 AM
    Tuesday, September 9, 2014 3:15 AM
  • Hello Bob, following is the updated query

    DECLARE @image1 geography
    DECLARE @image2 geography

    SET @image1=geography::STGeomFromText('POLYGON((66.91292909247741 30.27001012181008,66.99456841960638 30.27128639618252,66.99700791329992 30.1509623521339,66.91546772378466 30.14969219541345,66.91292909247741 30.27001012181008))',4326);
    SET @image2=geography::STGeomFromText('POLYGON((74.22257407407407 32.21185185185185,74.46148148148148 32.21185185185185,74.46148148148148 32.04995833333334,74.22257407407407 32.04995833333334,74.22257407407407 32.21185185185185))',4326);

    Declare @fg geography=geography::STGeomFromText('FULLGLOBE',4326);
    Declare @a float=@fg.STArea()/2;
    declare @myGeo1 geography=@image1
    declare @myGeo2 geography=@image2
    if (@myGeo1.STArea()>@a)
        set @image1 =@myGeo1.ReorientObject();
    if (@myGeo2.STArea()>@a)
        set @image2 =@myGeo1.ReorientObject();

    SELECT CASE @image1.STIntersects(@image2)
    WHEN 1 THEN 'cutting'
    ELSE 'Not cutting'
    END;

    This always return 'cutting', but both polygons are not cutting each other if visualized in Google earth. If I changed data type from geography to geometry then it gives correct results.

    Thanks


    Engr. Mudassar Ali Software Engineer

    Tuesday, September 9, 2014 4:29 AM
  • Hi Mudassar,

    There's a subtle bug in the code. Make the change below and it works. BTW, both image1 and image2 did need to be reoriented, which is likely why you had more "false positives" in STIntersects using the cursor originally.

    DECLARE @image1 geography
    DECLARE @image2 geography

    SET @image1=geography::STGeomFromText('POLYGON((66.91292909247741 30.27001012181008,66.99456841960638 30.27128639618252,66.99700791329992 30.1509623521339,66.91546772378466 30.14969219541345,66.91292909247741 30.27001012181008))',4326);
    SET @image2=geography::STGeomFromText('POLYGON((74.22257407407407 32.21185185185185,74.46148148148148 32.21185185185185,74.46148148148148 32.04995833333334,74.22257407407407 32.04995833333334,74.22257407407407 32.21185185185185))',4326);

    Declare @fg geography=geography::STGeomFromText('FULLGLOBE',4326);
    Declare @a float=@fg.STArea()/2;
    declare @myGeo1 geography=@image1
    declare @myGeo2 geography=@image2
    if (@myGeo1.STArea()>@a)
        set @image1 =@myGeo1.ReorientObject();
    if (@myGeo2.STArea()>@a)
        set @image2 =@myGeo1.ReorientObject();   -- should be @myGeo2.ReorientObject()

    SELECT CASE @image1.STIntersects(@image2)
    WHEN 1 THEN 'cutting'
    ELSE 'Not cutting'
    END;

    Cheers, Bob

    • Marked as answer by Sofiya Li Wednesday, September 10, 2014 9:05 AM
    Tuesday, September 9, 2014 4:16 PM
  • Thats Great. Thanks for mentioning mistake that I was doing. One more thing I want to ask that currently I have used geometry data type for my searched query and its giving me correct results, now whether I used the same data type of convert it to geography one? which one will be the most appropriate?

    Second is it any way that I can draw search results polygons from vertices that are cutting my required area of interest to the front end? (on windows form using c#) same way that the sql draw polygon when STUnion() used for spatial results?

    Thanks

    Engr. Mudassar Ali Software Engineer

    Wednesday, September 10, 2014 3:47 AM
  • Glad that it worked out. That typo took me a while to find too. ;-)

    Which data type you want to use depends on what you're doing with it. If you're measuring, for example, the distance between Seattle and Boston, the geography type takes the shape of the earth into consideration, geometry doesn't. But for some things, like drawing shapes on a map or getting the distances between two locations that are near each other, geometry is fine. In fact, geometry is faster to do calculations, because it doesn't consider the earth's shape the math is easier and faster.

    I'm not a GUI specialist with spatial. I work almost exclusively with the database, specializing in getting best performance with big spatial tables and using spatial indexes. And the spatial types in SQL Server don't always map one-to-one with the classes in .NET's System.Drawing. If you're using real places on a map you'll be better off imbedding the Bing map control in a C# program. or using Google or Bing map APIs directly in a web page. You might find some good ideas in Ricky's Bing Map Blog (http://rbrundritt.wordpress.com/) or have a look at the code in Alastair Aitchison's (http://alastaira.wordpress.com/) book on SQL Server spatial, I believe he was Bing and Google API examples.

    Finally, there's an analogous forum for Bing map control programming here: http://social.msdn.microsoft.com/Forums/onedrive/en-US/home?forum=vemapcontroldev 

    Hope this helps get you started, Cheers, Bob




    • Edited by Bob Beauchemin Wednesday, September 10, 2014 5:04 AM
    • Marked as answer by Sofiya Li Wednesday, September 10, 2014 9:05 AM
    Wednesday, September 10, 2014 4:57 AM
  • One more thing to clarify Bob. is it possible that I can return results with the percentage value the first image cutting to the second? Percentage mean area of the first image that the second one is cutting?

    Engr. Mudassar Ali Software Engineer

    Wednesday, September 17, 2014 5:26 PM
  • Sure. You can use STIntersection to get a polygon that represents the intersection of the two polygons, then use STArea and compare the intersection area vs. total area.

    Cheers, Bob

    • Marked as answer by Mudassar Alii Thursday, September 18, 2014 10:05 AM
    Thursday, September 18, 2014 12:10 AM
  • Thats awesome. Thanks

    Engr. Mudassar Ali Software Engineer

    Thursday, September 18, 2014 10:05 AM