SDE MSSQL SqlSpatial Geometry Self Intersects problem

[ChrisSmith]

Hi

I have created a SqlGeometry using the Microsoft.SqlServer.Types SqlGeometryBuilder, the geometry IsValid and IsSimple but when trying to view it in ArcMap it wont draw, when trying to export to a non SqlSpatial SDE table it complains that the polygon is self intersecting. The geometry uses srid 27700 (BNG) and is viewable in SQL management Studio, I have tried various methods that are suggested to correct the geometry including MakeValid(), STUnion(STStartPoint()) on the SqlSpatial side but since it is valid and doesnt self intersect then it doesnt fix the issue.

Any help is appreciated, Geometry details below



Regards

chris

Geometry as Text...(after trying MakeValid()/StUnion(STStartPoint))

POLYGON ((531405.73999977112 104148.89000034332, 531414.71999979019 104179.65999984741, 531407.15999984741 104182.27999973297, 531406.94000005722 104182.35000038147, 531404.44000005722 104177.02000045776, 531401.59000015259 104168.47000026703, 531401.59999990463 104168.5, 531401.48999977112 104168.19999980927, 531395.61999988556 104152.19999980927, 531405.73999977112 104148.89000034332))

Geometry as Gml...

<Polygon xmlns="http://www.opengis.net/gml">
  <exterior>
    <LinearRing>
      <posList>531405.73999977112 104148.89000034332 531414.71999979019 104179.65999984741 531407.15999984741 104182.27999973297 531406.94000005722 104182.35000038147 531404.44000005722 104177.02000045776 531401.59000015259 104168.47000026703 531401.59999990463 104168.5 531401.48999977112 104168.19999980927 531395.61999988556 104152.19999980927 531405.73999977112 104148.89000034332</posList>
    </LinearRing>
  </exterior>
</Polygon>



Original gml entry from file I am importing...

<gml:Polygon srsName='osgb:BNG'>
<gml:outerBoundaryIs>
<gml:LinearRing>
<gml:coordinates>531406.94,104182.35 531404.44,104177.02 531401.59,104168.47 531401.60,104168.50 531401.49,104168.20 531395.62,104152.20 531405.74,104148.89 531414.72,104179.66 531407.16,104182.28 531406.94,104182.35</gml:coordinates>
</gml:LinearRing>
</gml:outerBoundaryIs>
</gml:Polygon>

[ChrisSmith]

Hi

I have loaded the geometry with the problem interior ring into a shapefile, it is clearly self intersecting when viewed in arcmap.

points from geometry in shapefile....
Part 0
529785.170 180867.590
529767.100 180902.360
529827.800 180940.950
529850.040 180907.640
529786.560 180868.450
529786.560 180868.450

Part 1
529786.560 180893.410
529811.730 180908.530
529807.560 180915.420
529786.330 180902.810
529790.100 180895.480

I then export the feature into sde and the process adds extra points I assume to fix the self intersecting issue.

The exported geometry is viewable in arcmap but if you look at point 1 & 4 they have the same values; is this not considered self intersecting?

points for geometry in sde...
Part 0
529785.170 180867.590
529767.100 180902.360
529827.800 180940.950
529850.040 180907.640
529786.560 180868.450

Part 1
529786.330 180902.810
529790.078 180895.523
529786.560 180893.410
529790.100 180895.480
529790.078 180895.523
529811.730 180908.530
529807.560 180915.420

cheers

chris

[VinceAngelo]

The SE_shape_generate_polygon is reasonably adaptive. In this case, if you add the
missing closing vertex to the array, the API will nominate the ring that begins and
ends at 529790.078,180895.523 as a subpart (hole), then complete the primary ring.
Having determined that the subpart is not within its original parent, it is then promoted
to an independent part. Parts and subparts are permitted to touch at a single point
(but not share boundaries or area). The API will also demote rings which *are* within
other parts, making them subparts.

I was the course author for the (now retired) "Accessing SDE with 'C'" 3-day training
class. We spent the second half of the first day on Clemetini topology and relationships.
There was exercise code which permitted generation and comparison of shapes from
command-line arguments. Those exercises are gone, but the core parsing functions
live on, exposed in the se_toolkit DAT module, which is how I was able to assemble
and test your geometry.

- V

[ChrisSmith]

Hi Vince

Once again thanks for the help!

I look forward to resurrecting my C-Api skills, been about 5 years but should be fun 😉

Do you still have the course notes that accompanied the exercises in particular the Clementini Topolgy section - or have a pointer to some reading material available on the web?

Happy New Year to you

chris

[VinceAngelo]

I can't give away even old training materials (mostly bullet lists on which the instructor
needed to expound, anyway). All the concepts are available in the documentation, or
in User Conference proceeedings (1997 and/or 1998 -- I *did* say this was mature
technology, right?)

- V

BTW: The original "Clementini" paper, which I'll often go back to, is:

E. Clementini, P. Di Felice, and P. van Oosterom, "A Small Set of Formal Topological Relationships Suitable for End-User Interaction," in Advances in Spatial Databases - Third International Symposium, SSD '93. vol. 692, D. Abel and B. C. Ooi, Eds. Berlin: Springer-Verlag, 1993, pp. 277-295.

[EdKatibah]

Hi Vince -

In the SQL Server Spatial world, this geometry does not self-intersect.  Michael Kallay of the SQL Server Spatial development team has posted an explanation on this issue to the SQL Server Spatial Forum (http://social.msdn.microsoft.com/Forums/en-US/sqlspatial/thread/d844071c-68b6-4ce8-bbda-05f19ac4b9e0...).

Here is what the geometry looks like around Point 7, the suspected self-intersecting point, using SQL Server and the FME Universal Viewer:



I will work with the the SDE engineering team to look at how SQL Server Spatial and ArcGIS can more gracefully handle these kinds of situations.

Thanks,

- Ed Katibah

SQL Server Spatial

[VinceAngelo]

The principal problem with using arbitrary precision with real-world coordinate data
is the implicit assertion of accuracy not present at input. The orginal data was
collected at centimeter scale. Varying the locations by a few hundred nanometers
can change topology, but if it does, then it becomes a data quality issue -- *should*
we trust spatial data for which the Heisenberg uncertainty principle applies?

There are a number of ways to process this data to remove the problematic vertices.
The simplest to code would be to perform angle calculation over three coordinate
sets -- vertices that belong in a shape at that scale would not bend the line by more
than 165 degrees (+/-) -- Whether this line bends by 180 degrees or 179.999999
degrees wouldn't matter in this sort of algorithm, just that you remove the vertex
that produces the shortest segment.

A more robust processing step would be to apply a full Douglas-Peucker generalization
to the polygon ring as a line, with a tolerance distance of, say, 1/20th the diagonal
of the bounding box (or just a fixed value, like 0.5 meters), then attempt to cast
the rings back into a polygon.

I have a tool that may be able to do this sort of processing (in a timely fashion with
millions of rows), but I'd have to play with an ASCII dump of the well-known polygon
text to figure out if this was a 20-minute hack or something more involved.

- V

[MichaelKallay]

The last post suggests that the MS SQL Server exact computation is akin to splitting hair, and the ESRI computation with a margin of error is better in practice.  Although the suggested margins (165 degrees or .5m or 1/20 of the bounding rectangle's diagonal) are much too broad for my liking, it may be a valid choice, but it is not necessarily better.  Computations with margins cannot be done in a fully consistent way, and therefore cannot be made fully robust. The problem is articulated at http://people.csail.mit.edu/seth/pubs/taskforce/section3_10.html.  The only known way to achieve full robustness is by using exact arithmetic.  Since we at MS consider robustness paramount, we chose the (more sophisticated and more difficult) exact arithmetic approach, and we are very proud of its clean and fast implementation.  Moreover, in our next relase (already available in a preview version) we have increased the resolution (of the integer space in which we perform our exact computations) from 27 to 48 bits of accuracy.

Michael Kallay, SQL Server spatial team

[VinceAngelo]

ArcSDE has been using exact arithmetic for geometry computation since the technology
was aquired in 1995.  What is now known as BASIC precision uses 31 bits; HIGH precision
(released with 9.0, and made default since 9.2) uses 64-bit integers (but is 53-bit, due to
double precision mantissa size).  The system is very flexible for the different data units found
in international databases; the default values chosen by ArcGIS are optimized by coordinate
system (projection) for preservation of real-world spatial data to millimeter precision (which
is roughly one hundred times the accuracy of most of the spatial datasets I've ever used.) .

The distance parameter in Douglas-Peucker doesn't correspond to precision, but to deflection;
value selection is data dependent (and, to some degree, more art than science).

- V

[ChrisSmith]

Hi Vince

Have been away for a week and it appears I have missed quite a bit!

Anyway what I need to be able to do is have a reasonably fast process to scan through 20+ million features to find out which are considered invalid by ESRI software.

Are you able to point me in the right direction? A function I can include in my code, an external library to reference or an executable I can run to select the problem geometries out would be desirable.

cheers

chris

[VinceAngelo]

It looks like the 'sdelayer -o feature_info' capability isn't functional for SQL-Server.

What does 'sdelayer -o describe_long' report for the coordinate reference on your layer?
What command did you use to register the layer?

I prototyped a workaround, but it's not going to be pretty. You'll need a copy of se_toolkit
for the 'sdequery' and 'ascinfo' utilities, and you'll need to create a view of the layer that
dumps the well-known text of the geometries as either string (if it will fit in varchar(8000))
or text. Then you can query the registered rowid column and WKT with 'sdequery', and
pipe it through 'ascinfo' to reject the rows that fail SE_shape_generate_from_text, producing
an ASCII file with ID and problematic well-known text. If you want to go this route, it would
probably be easier to take this offline.

Unfortunately, I didn't have time to whip up some data massage classes to automate removal
of spurious jagged lines in geometry.

- V