Sunday, February 13, 2011

Spikes: research with ESRI


Spikes: research with ESRI


A short follow-up on my previous blog about spikes and their causes. My collegue Bert recoded the scenario in C# using ESRI. Thanks Bert! The full code is shown below.

The answer is this: with ESRI the spikes are not there. The program produces:

area: 295.44447138901, perimeter: 92.5058167188044

select geometry::STGeomFromText ('POLYGON((180956.436999999
313716.998,180953.978700001 313701.0152,180926.420299999
313732.510699999,180954.359000001
313720.504000001,180956.436999999 313716.998,180956.436999999
313716.998))', 28992)
And this result is right, area and perimeter are without spikes and if you visualize this WKT (e.g. with SQL Server) you get the right picture...

Interesting is that all coordinates seems to be rounded on one millimeter or a tenth of it. That is probably the result of the fuzzy tolerance policy that ESRI uses. Personally I'm normally a bit sceptic about that approach: though it can avoid spikes and slivers, it also brings some fuzzyness into the result. It is sometimes quite hard to determine a tolerance which is sufficient for cleaning and sufficient for not destroying features... Anyway, in this scenario it works well.

It seems interesting to check, if you use rounding in e.g. Boost.Geometry, would the spikes be gone... The answer is (of course), sometimes yes, sometimes no. The rounding is applied on a larger grid, so the effect remains the same. Rounding is not enough. To get rid of spikes, in this scenario, the spike removal functionality is advised. Which uses a gap distance, which can be seen as a sort of fuzzy tolerance...

C# Code for "difference after cut" with ESRI


using System; 
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Text;
using System.Windows.Forms;
using ESRI.ArcGIS.Geometry;
using ESRI.ArcGIS.esriSystem;
using Microsoft.SqlServer.Types;
using System.Data.SqlTypes;

namespace SpikesResearchWithEsri
{
    public partial class Form1 : Form
    {
        public Form1()
        {
          InitializeComponent();
        }

        private void InitializeLicense()
        {
          AoInitialize aoi = new AoInitializeClass();

          //Additional license choices can be included here.
          esriLicenseProductCode productCode =
          esriLicenseProductCode.esriLicenseProductCodeArcEditor;
          if (aoi.IsProductCodeAvailable(productCode) ==
          esriLicenseStatus.esriLicenseAvailable)
          {
          aoi.Initialize(productCode);
          }
        }


        private void button1_Click(object sender, EventArgs e)
        {


          ESRI.ArcGIS.RuntimeManager.Bind(ESRI.ArcGIS.ProductCode.Desktop);
          InitializeLicense();
          IPolygon parcel = this.GetParcel();
          IPolyline line = this.GetCutLine();

          IPolygon bufferedLine=(IPolygon)((ITopologicalOperator)line).Buffer(0.1);
         
          IGeometry result=((ITopologicalOperator)parcel).Difference(bufferedLine);
          IPolygon4 p = (IPolygon4)result;

          // Big part:
          IPolygon firstPolygon = this.GetPolygon(p, true);
          // Small part:
          IPolygon secondPolygon = this.GetPolygon(p, false);


          // use if we do a cutpolygon
          //IGeometry leftPolygon;
          //IGeometry rightPolygon;
          //((ITopologicalOperator)parcel).Cut(line, out leftPolygon, out rightPolygon);
          //IPolyline cutline=(IPolyline)((ITopologicalOperator)(ICurve)parcel).Intersect(line, esriGeometryDimension.esriGeometryNoDimension);
          //IPolygon firstPolygon = (IPolygon)leftPolygon;
          //IPolygon secondPolygon = (IPolygon)rightPolygon;
          //this.textBox1.Text=("from: " + cutline.FromPoint.X.ToString() + ", " +
          // cutline.FromPoint.Y.ToString() + ", " +
          // "to: " + cutline.ToPoint.X.ToString() + ", " +
          // cutline.ToPoint.Y.ToString()) + Environment.NewLine;

          this.textBox1.Text += "Small polygon stats: " + System.Environment.NewLine;
          this.textBox1.Text += GetStats(secondPolygon);
          this.textBox1.Text += GetWkt(secondPolygon);

          // now calculate the difference between original and the big part (result should be the small part)
          IPolygon4 pol2= (IPolygon4)((ITopologicalOperator)parcel).Difference(firstPolygon);
          this.textBox1.Text += "Result Small polygon stats: " + System.Environment.NewLine;
          this.textBox1.Text += GetStats(pol2);
          this.textBox1.Text += GetWkt(pol2);
        }

        private string GetStats(IPolygon pol)
        {
          double area = ((IArea)pol).Area;
          double perimeter = ((ICurve)pol).Length;
          return "area: " + area.ToString() + ", perimeter: " + perimeter.ToString() + System.Environment.NewLine;
        }


        private string GetWkt(IPolygon pol)
        {
          string res = GeomHelper.aoGeomToWkt(pol, false, 28992);
          res = res.Substring(0, res.Length - 1);
          string print = "select geometry::STGeomFromText ('" + res + "', 28992)" + Environment.NewLine;
          return print;
        }

        public IPolygon GetPolygon(IPolygon4 input, bool first)
        {
          IGeometryBag exteriorRings = input.ExteriorRingBag;
          IEnumGeometry exteriorRingsEnum = exteriorRings as IEnumGeometry;
          exteriorRingsEnum.Reset();

          IRing firstRing = exteriorRingsEnum.Next() as IRing;
          IPointCollection pcFirst = (IPointCollection)firstRing;

          IRing secondRing = exteriorRingsEnum.Next() as IRing;
          IPointCollection pcSecond = (IPointCollection)secondRing;

          IPolygon p = new PolygonClass();

          if(first)((IPointCollection)p).AddPointCollection(pcFirst);
          if(!first)((IPointCollection)p).AddPointCollection(pcSecond);
          return p;
        }


        private IPolygon GetParcel()
        {
          string wkt = "POLYGON ((180956.437 313716.998, 180954.359 313720.504, 180906.303 313741.156, 180896.82 313744.272, 180890.456 313748.17, 180831.49 313797.135, 180809.668 313816.098, 180776.418 313842.336, 180771.483 313846.102, 180746.934 313843.764, 180731.738 313840.386, 180717.32 313833.763, 180710.827 313828.438, 180704.982 313820.515, 180699.267 313811.813, 180698.488 313801.162, 180700.826 313793.889, 180705.539 313791.06, 180706.021 313790.771, 180710.156 313789.208, 180766.157 313768.042, 180825.515 313744.013, 180878.636 313724.271, 180945.593 313695.884, 180952.801 313693.358, 180956.437 313716.998))";
          IGeometryFactory factory = new GeometryEnvironmentClass();
          SqlGeometry g = SqlGeometry.Parse(wkt);
          IGeometry geom2 = new PolygonClass();
          int countout;
          factory.CreateGeometryFromWkbVariant(g.STAsBinary().Value, out geom2, out countout);
          geom2.SpatialReference = this.GetProjectedSpatialReference(28992);
          return (IPolygon) geom2;
        }

        private IPolyline GetCutLine()
        {
          string line = "LINESTRING(180955 313700,180920 313740)";
          IGeometryFactory factory = new GeometryEnvironmentClass();
          SqlGeometry g = SqlGeometry.Parse(line);
          IGeometry geom2 = new PolylineClass();
          geom2.SpatialReference = this.GetProjectedSpatialReference(28992);
          int countout;
         
          factory.CreateGeometryFromWkbVariant(g.STAsBinary().Value, out geom2, out countout);
          return (IPolyline)geom2;

        }


        private ISpatialReference GetProjectedSpatialReference(int pcsType)
        {
          ISpatialReferenceFactory pSRF = new SpatialReferenceEnvironmentClass();
          IProjectedCoordinateSystem m_ProjectedCoordinateSystem = pSRF.CreateProjectedCoordinateSystem(pcsType);
          ISpatialReference spatialReference = (ISpatialReference)m_ProjectedCoordinateSystem;
          return spatialReference;
        }

    }
}




3 comments:

  1. the difference from 151 to 92 is perfectly explicable:
    after subtracting the buffered line from the original polygon , one obtains 2 smaller polygons;
    in the case of sqlserver and postgis, you subtract again one smaller polygon from THE ORIGINAL polygon;
    in the case of ESRI , you calculate the length of a smaller polygon

    conclusion: in the case of sqlserver and postgis, the length of the buffered line (intersected with the original polygon),
    is added to the length of the smaller polygon. this explains the difference from 92 to 151: it is the length of the buffered line.

    I have verified this again and again and again , and it computes every time.

    ReplyDelete
  2. This comment has been removed by the author.

    ReplyDelete
  3. In Alastair’s book “ProSpatial with SQL Server 2012”, at page 173 the author rightly and brilliantly says that
    "SQL Server takes a similar approach; every spatial calculation takes place on an integer grid of
    fixed size, equivalent to the fixed size sheet of graph paper...
    The amount of error introduced in this process is determined by two main factors:
    • The overall size of the grid on which calculations are performed.
    • The extent of the geometries involved in a given calculation.
    As for the second factor, the greater the range of coordinate values that must be covered by the
    integer grid, then the more coarse the grid resolution must become in order to accommodate the full
    set of data. As each cell becomes larger and the distance between cells increases, there is a greater
    distance that coordinates might potentially be snapped. This explains the effect demonstrated in the
    following code listing, in which the area of intersection between the square and rectangular Polygon
    became less and less precise as the overall size of the square became larger and larger, because the
    integer grid had to become ever more coarse to accommodate it, leading to less accurate results."
    DECLARE @rectangle geometry = 'POLYGON((-10 5, 10 5, 10 15, -10 15, -10 5))';
    DECLARE @square geometry = 'POLYGON((0 0, 1000 0, 1000 1000, 0 1000, 0 0))';
    SELECT @rectangle.STIntersection(@square).STArea();
    -- 99.9999999999713
    DECLARE @square2 geometry = 'POLYGON((0 0, 100000 0, 100000 100000, 0 100000, 0 0))';
    SELECT @rectangle.STIntersection(@square2).STArea();
    -- 99.9999999962893
    DECLARE @square3 geometry = 'POLYGON((0 0, 1e9 0, 1e9 1e9, 0 1e9, 0 0))';
    SELECT @rectangle.STIntersection(@square3).STArea();
    -- 99.9999690055833
    DECLARE @square4 geometry = 'POLYGON((0 0, 1e12 0, 1e12 1e12, 0 1e12, 0 0))';
    SELECT @rectangle.STIntersection(@square4).STArea();
    -- 99.9691756255925
    DECLARE @square5 geometry = 'POLYGON((0 0, 1e15 0, 1e15 1e15, 0 1e15, 0 0))';
    SELECT @rectangle.STIntersection(@square5).STArea();
    -- 67.03125

    Now, what if i tell you i have managed to obtain , in c#.code, the csharp-EQUIVALENT-of:
    SELECT @rectangle.STIntersection(@square5).STArea();
    -- 100.000000
    ? do you believe me? do you believe that i'm passing this ULTIMATE test ?
    i must say that i am not using a comercial product , just the csharp equivalent of your boost.geometry (angus johnson's clipper)')

    ReplyDelete