Qt World Mapper with Boost.Geometry

Last year I created a simple world mapper using WxWidgets. For various reasons I decided to make the Qt equivalent as well.

Both Qt and WxWidgets are well-known cross-platform windows frameworks. Qt is used in Quantum GIS.

The World Mapper I created with WxWidgets did read world countries from a WKT (well-known text) file, and displayed them. Moving the mouse did highlight the country. For Qt it is similar (the highlighting is not yet been done).

The nice thing is that, by setting Boost.Geometry in this context, its force is shown:

Look below to lines marked with // Adapt. The QPointF and QPolygonF are just registered. From now on they act as normal Boost.Geometry geometries. So they can be used in area calculation, distance, transform, etc. We only use transform in this example.

So also look at lines commented with // This is the essention. Here a QPolygonF is declared, a ring (a normal Boost.Geometry part of a polygon, part of a multi-polygon) is transformed using the transformer to a QPolygonF. And the last one is of course displayed without any problem.

Delivering this application:

Complete C++ Code

#include <fstream> 

#include <QtGui> 
#include <QWidget> 
#include <QObject> 
#include <QPainter> 

#include <boost/foreach.hpp> 

#include <boost/geometry/geometry.hpp> 
#include <boost/geometry/geometries/register/point.hpp>
#include <boost/geometry/geometries/register/ring.hpp>

#include <boost/geometry/multi/multi.hpp> 
#include <boost/geometry/extensions/algorithms/selected.hpp>
#include <boost/geometry/extensions/gis/io/wkt/wkt.hpp> 

// Adapt a QPointF such that it can be handled by Boost.Geometry
BOOST_GEOMETRY_REGISTER_POINT_2D_GET_SET(QPointF, double, cs::cartesian, x, y, setX, setY)

// Adapt a QPolygonF as well. 
// A QPolygonF has no holes (interiors) so it is similar to a Boost.Geometry ring
BOOST_GEOMETRY_REGISTER_RING(QPolygonF) 

typedef boost::geometry::model::d2::point_xy<double> point_2d; 
typedef boost::geometry::model::multi_polygon 
    < 
        boost::geometry::model::polygon<point_2d> 
    > country_type; 


class WorldMapper : public QWidget
{ 
public: 
    WorldMapper(std::vector<country_type> const& countries, boost::geometry::model::box<point_2d> const& box)
        : m_countries(countries) 
        , m_box(box) 
    { 
        setPalette(QPalette(QColor(200, 250, 250)));
        setAutoFillBackground(true); 
    } 

protected: 
    void paintEvent(QPaintEvent*) 
    { 
        map_transformer_type transformer(m_box, this->width(), this->height()); 

        QPainter painter(this); 
        painter.setBrush(Qt::green); 

        BOOST_FOREACH(country_type const& country, m_countries)
        { 
          typedef boost::range_value<country_type>::type polygon_type;
          BOOST_FOREACH(polygon_type const& polygon, country)
          { 
           typedef boost::geometry::ring_type<polygon_type>::type ring_type;
           ring_type const& ring = boost::geometry::exterior_ring(polygon);

           // This is the essention: 
           // Directly transform from a multi_polygon (ring-type) to a QPolygonF 
           QPolygonF qring; 
           boost::geometry::transform(ring, qring, transformer); 

           painter.drawPolygon(qring); 
          } 
       } 
    } 

private: 
    typedef boost::geometry::strategy::transform::map_transformer
        < 
          point_2d, QPointF, 
          true, true 
        > map_transformer_type; 

     std::vector<country_type> const& m_countries;
     boost::geometry::model::box<point_2d> const& m_box;
}; 


class MapperWidget : public QWidget
{ 
public: 
     MapperWidget(std::vector<country_type> const& countries, boost::geometry::model::box<point_2d> const& box, QWidget *parent = 0)
         : QWidget(parent) 
     { 
         WorldMapper* mapper = new WorldMapper(countries, box);

         QPushButton *quit = new QPushButton(tr("Quit")); 
         quit->setFont(QFont("Times", 18, QFont::Bold)); 
         connect(quit, SIGNAL(clicked()), qApp, SLOT(quit())); 

         QVBoxLayout *layout = new QVBoxLayout; 
         layout->addWidget(mapper); 
         layout->addWidget(quit); 
         setLayout(layout); 
     } 

}; 


template <typename Geometry, typename Box>
inline void read_wkt(std::string const& filename, std::vector<Geometry>& geometries, Box& box)
{ 
    std::ifstream cpp_file(filename.c_str()); 
    if (cpp_file.is_open()) 
    { 
        while (! cpp_file.eof() ) 
        { 
          std::string line; 
          std::getline(cpp_file, line); 
          if (! line.empty()) 
          { 
              Geometry geometry; 
              boost::geometry::read_wkt(line, geometry); 
              geometries.push_back(geometry); 
              boost::geometry::combine(box, boost::geometry::make_envelope<Box>(geometry)); 
          } 
        } 
    } 
} 

int main(int argc, char *argv[])
{ 
    std::vector<country_type> countries; 
    boost::geometry::model::box<point_2d> box; 
    boost::geometry::assign_inverse(box); 
    read_wkt("../data/world.wkt", countries, box); 

    QApplication app(argc, argv); 
    MapperWidget widget(countries, box); 
    widget.setWindowTitle("Boost.Geometry for Qt - Hello World!"); 
    widget.setGeometry(50, 50, 800, 500);
    widget.show(); 
    return app.exec(); 
}

Conclusions

With less than 160 lines of code, Qt and Boost.Geometry we showed how it can quite conveniently combine together.
The source is shown above and in SVN it can be found here.

Actually this application (and its WxWidgets counterpart) was for me the reason to start the shapefile research, which final conclusions still have to be drawn.

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; 
        } 

    } 
} 

Precision, the cause of spikes

This is a follow-up on my last blog. The spikey effects we have seen are perfectly explainable and understandable.

What we did was:

1. cut a parcel in two pieces
1. create a cutline and make a (nearly linear) cutting polygon of it (using buffer)
2. subtract that cutline from the geometry, using STDifference
2. keep one piece, using STGeometryN
3. subtract that piece from the parcel, using STDifference

The result should be the discarded piece, but we got back two free spikes extra.

This short blog will explain this a little more and, as a side note, use just another SQL syntax for the same result.

Side note: with

The query from my previous blog can be rewritten much more elegantly using the with syntax:

with myquery as 
(
  select
    geometry::STGeomFromWKB(0x0103000000010000001A00000023DBF97EE316064146B6F3FDD32513415A643BDFD216064175931804E225134196438B6C52150641C976BE9F34261341F6285C8F06150641022B871641261341F853E3A5D3140641E17A14AE50261341B81E85EBFB120641A4703D8A142713414E6210584D120641AC1C5A64602713414E621058431106414E621058C9271341A01A2FDD1B11064121B07268D8271341F4FDD478571006411904560ECF271341448B6CE7DD0F06418195438BC1271341F6285C8F6A0F06413BDF4F0DA72713410E2DB29D360F06416F1283C091271341E5D022DB070F0641F6285C0F7227134160E5D022DA0E06416F1283404F271341448B6CE7D30E0641F853E3A52427134154E3A59BE60E06411904568E07271341643BDF4F0C0F0641D7A3703DFC2613414A0C022B100F064125068115FB26134191ED7C3F310F0641B6F3FDD4F42613414C378941F11006414A0C022BA0261341EC51B81ECC1206413BDF4F0D40261341022B87167514064125068115F1251341B4C876BE8C160641C74B37897F25134121B07268C6160641508D976E7525134123DBF97EE316064146B6F3FDD3251341,28992)
      as parcel,
    geometry::STGeomFromText('LINESTRING(180955 313700,180920 313740)', 28992).STBuffer(0.1)
      as cutline
)
select parcel.STDifference((select parcel.STDifference(cutline).STGeometryN(2) from myquery))
from myquery

The with clause is useful to avoid repeating pieces (as having the WKB twice), and is essential to do recursive queries. You define an inlined viewy thingy (here "myquery") and use it.

Both SQL Server and PostGreSQL support with queries. See also this link which explains with in depth.

Floating point precision

Floating points are really awful. Everyone knows that if you divide one by six (1/6) and multiply the result with six ((1/6) * 6) you should get one back. However, with floating points this is not always the case. You might get 0.9999999 back. If you try it, you might also get the correct 1.0 back. That is (according to this book) because intermediate results are stored in the floating point register of the CPU. So it is by chance that it is right. In reality it will be wrong.

See also this link and scroll to "Rounding errors can play havoc with math-intense programs". Adding ten 0.1 values will result into 0.99999999999999989. I repeated this and indeed this is the case with a double. With a float the result is 1.0000001192092896. With ttmath it is 1.0, but as soon as you do the same trick with eleven 1/11 additions, you might get 0.9999999999999999999999999999999 back, or 1.0, depending on which precision you specify. Using Boost.Rational you get always 1/1, guaranteed, because it calculates in another way. It uses fractions and keeps fractions all the time. That is great in this use case.

What then caused the spikes

Now that we realize again that floating points are awful, we investigate how imprecision causes the spikes. The damage happens in the first phase. The first STDifference is used to cut the parcel into two pieces.

See the picture below. The grid is the so-called floating point grid. The smallest difference in a 32 bits floating point value (1.175e-38) is depicted. Any point (with float) must be located on one of the gridpoints. (In reality the FP-grid is irregular, but the idea will be clear).

So in this case (in most cases), the intersection point happens to fall between the gridpoints... All the vertices of both input polygons are located on the gridpoints, of course. But the real intersection point does not. It is rounded to one of the four gridpoints. It is just by chance to which one it is rounded. So the rounded intersection point is located on a gridpoint.




Suppose it rounds to the point labeled "3". We then get a resulting polygon as partly depicted here:




The resulting polygon (in blue) is too small. There is a small piece of the original polygon which is not covered by this piece, at the north side, upperleft of point 3. Subtract this resulting polygon from the original and of course, you will get the spike. There is also a similar small piece at the east side but that will not create a visual effect in this use case.

There is no obvious solution. It is the same as adding ten times 0.1. You get too much or too little. These spikes are too little in the first phase, and after the second difference, the two resulting spikes are a bit too much. This explains that sometimes you get zero, sometimes you get one and sometimes you get two spikes. It is by chance. If the rounding would go to point "1" in the picture above, there would be no spike there.

The spikes are not the result of flawed logic, or wrong decisions somewhere. It is just the rounding.

Use ttmath if you want to minimize the chance on it. Or work around it, e.g. by buffering the result with a small value e.g. -0.000001. Boost.Geometry now has some algorithms to find spikes and remove them.

By the way, the nice trick on cutting a parcel into two pieces was introduced in our project by our partners from B3Partners. The project where this study is from extracted is the Portaal Natuur en Landschap

Conclusions

In this use case, spikes are the geometric visualization of floating point imprecision.

SQL Server, PostGIS, STDifference and Precision

Today I encountered a problem related to the OGC difference function. I will describe this problem in this blog. I created a sort of minimal sample, without database, and present that here.

The problem

The problem I encountered is precision: both SQL Server (where I encountered it) and PostGIS (where I verified it) do not handle the situation described here sufficiently. The results is a polygon with two spikes which should not be there. The cause must be imprecision of floating point types. The message is: be warned!

The situation is as following: there is a parcel. That parcel is cut into two pieces (I will explain below how). One piece is thrown away, the other piece is kept. That kept piece is subtracted from the original parcel. What you would expect is having the original piece that was earlier thrown away. But what is happening is that you get one or two spikes extra...

The parcel and the cutline are shown here:

Cutting into two pieces

We cut (for this experiment) the polygon into two pieces by buffering the cutline, and subtracting (with STDifference) that from the parcel. The cutline is shown in the figure above as the dotted red line. Then we keep one of the two resulting polygons using the STGeometryN function. By trying we get the left piece which we want to keep by index 2.

So this piece of SQL is, for SQL Server:

select geometry::STGeomFromWKB(0x0103000000010000001A00000023DBF97EE316064146B6F3FDD32513415A643BDFD216064175931804E225134196438B6C52150641C976BE9F34261341F6285C8F06150641022B871641261341F853E3A5D3140641E17A14AE50261341B81E85EBFB120641A4703D8A142713414E6210584D120641AC1C5A64602713414E621058431106414E621058C9271341A01A2FDD1B11064121B07268D8271341F4FDD478571006411904560ECF271341448B6CE7DD0F06418195438BC1271341F6285C8F6A0F06413BDF4F0DA72713410E2DB29D360F06416F1283C091271341E5D022DB070F0641F6285C0F7227134160E5D022DA0E06416F1283404F271341448B6CE7D30E0641F853E3A52427134154E3A59BE60E06411904568E07271341643BDF4F0C0F0641D7A3703DFC2613414A0C022B100F064125068115FB26134191ED7C3F310F0641B6F3FDD4F42613414C378941F11006414A0C022BA0261341EC51B81ECC1206413BDF4F0D40261341022B87167514064125068115F1251341B4C876BE8C160641C74B37897F25134121B07268C6160641508D976E7525134123DBF97EE316064146B6F3FDD3251341,28992)
.STDifference(geometry::STGeomFromText('LINESTRING(180955 313700,180920 313740)', 28992).STBuffer(0.1)).STGeometryN(2)

And for PostGIS it is:

select ST_GeometryN(ST_Difference(
ST_SetSRID('0103000000010000001A00000023DBF97EE316064146B6F3FDD32513415A643BDFD216064175931804E225134196438B6C52150641C976BE9F34261341F6285C8F06150641022B871641261341F853E3A5D3140641E17A14AE50261341B81E85EBFB120641A4703D8A142713414E6210584D120641AC1C5A64602713414E621058431106414E621058C9271341A01A2FDD1B11064121B07268D8271341F4FDD478571006411904560ECF271341448B6CE7DD0F06418195438BC1271341F6285C8F6A0F06413BDF4F0DA72713410E2DB29D360F06416F1283C091271341E5D022DB070F0641F6285C0F7227134160E5D022DA0E06416F1283404F271341448B6CE7D30E0641F853E3A52427134154E3A59BE60E06411904568E07271341643BDF4F0C0F0641D7A3703DFC2613414A0C022B100F064125068115FB26134191ED7C3F310F0641B6F3FDD4F42613414C378941F11006414A0C022BA0261341EC51B81ECC1206413BDF4F0D40261341022B87167514064125068115F1251341B4C876BE8C160641C74B37897F25134121B07268C6160641508D976E7525134123DBF97EE316064146B6F3FDD3251341'
 ::geometry, 28992)
, (ST_Buffer(ST_GeomFromText('LINESTRING(180955 313700,180920 313740)',28992),0.1))), 2)

Note that the 28992 is just an SRID of a coordinate system (the Dutch), it probably works (or better stated: does not work) with any coordinate system. The result is this:

Subtracting one of the pieces from the original

After cutting we subtract the piece we kept from the original parcel. This will result in the other piece. In fact we get the area covered by the small buffer as well, but that is hardly visible, and is not important for this story.

This is how it should look like:



And this is the complete query in SQL Server:

select
geometry::STGeomFromWKB(0x0103000000010000001A00000023DBF97EE316064146B6F3FDD32513415A643BDFD216064175931804E225134196438B6C52150641C976BE9F34261341F6285C8F06150641022B871641261341F853E3A5D3140641E17A14AE50261341B81E85EBFB120641A4703D8A142713414E6210584D120641AC1C5A64602713414E621058431106414E621058C9271341A01A2FDD1B11064121B07268D8271341F4FDD478571006411904560ECF271341448B6CE7DD0F06418195438BC1271341F6285C8F6A0F06413BDF4F0DA72713410E2DB29D360F06416F1283C091271341E5D022DB070F0641F6285C0F7227134160E5D022DA0E06416F1283404F271341448B6CE7D30E0641F853E3A52427134154E3A59BE60E06411904568E07271341643BDF4F0C0F0641D7A3703DFC2613414A0C022B100F064125068115FB26134191ED7C3F310F0641B6F3FDD4F42613414C378941F11006414A0C022BA0261341EC51B81ECC1206413BDF4F0D40261341022B87167514064125068115F1251341B4C876BE8C160641C74B37897F25134121B07268C6160641508D976E7525134123DBF97EE316064146B6F3FDD3251341,28992)
.STDifference
(
  (
 select 
 geometry::STGeomFromWKB(0x0103000000010000001A00000023DBF97EE316064146B6F3FDD32513415A643BDFD216064175931804E225134196438B6C52150641C976BE9F34261341F6285C8F06150641022B871641261341F853E3A5D3140641E17A14AE50261341B81E85EBFB120641A4703D8A142713414E6210584D120641AC1C5A64602713414E621058431106414E621058C9271341A01A2FDD1B11064121B07268D8271341F4FDD478571006411904560ECF271341448B6CE7DD0F06418195438BC1271341F6285C8F6A0F06413BDF4F0DA72713410E2DB29D360F06416F1283C091271341E5D022DB070F0641F6285C0F7227134160E5D022DA0E06416F1283404F271341448B6CE7D30E0641F853E3A52427134154E3A59BE60E06411904568E07271341643BDF4F0C0F0641D7A3703DFC2613414A0C022B100F064125068115FB26134191ED7C3F310F0641B6F3FDD4F42613414C378941F11006414A0C022BA0261341EC51B81ECC1206413BDF4F0D40261341022B87167514064125068115F1251341B4C876BE8C160641C74B37897F25134121B07268C6160641508D976E7525134123DBF97EE316064146B6F3FDD3251341,28992)
 .STDifference(geometry::STGeomFromText('LINESTRING(180955 313700,180920 313740)', 28992).STBuffer(0.1)).STGeometryN(2)
 )
)

And in PostGIS:

select ST_Perimeter(ST_Difference(
ST_SetSRID('0103000000010000001A00000023DBF97EE316064146B6F3FDD32513415A643BDFD216064175931804E225134196438B6C52150641C976BE9F34261341F6285C8F06150641022B871641261341F853E3A5D3140641E17A14AE50261341B81E85EBFB120641A4703D8A142713414E6210584D120641AC1C5A64602713414E621058431106414E621058C9271341A01A2FDD1B11064121B07268D8271341F4FDD478571006411904560ECF271341448B6CE7DD0F06418195438BC1271341F6285C8F6A0F06413BDF4F0DA72713410E2DB29D360F06416F1283C091271341E5D022DB070F0641F6285C0F7227134160E5D022DA0E06416F1283404F271341448B6CE7D30E0641F853E3A52427134154E3A59BE60E06411904568E07271341643BDF4F0C0F0641D7A3703DFC2613414A0C022B100F064125068115FB26134191ED7C3F310F0641B6F3FDD4F42613414C378941F11006414A0C022BA0261341EC51B81ECC1206413BDF4F0D40261341022B87167514064125068115F1251341B4C876BE8C160641C74B37897F25134121B07268C6160641508D976E7525134123DBF97EE316064146B6F3FDD3251341'::geometry,28992)
,
(select ST_GeometryN(ST_Difference(
ST_SetSRID('0103000000010000001A00000023DBF97EE316064146B6F3FDD32513415A643BDFD216064175931804E225134196438B6C52150641C976BE9F34261341F6285C8F06150641022B871641261341F853E3A5D3140641E17A14AE50261341B81E85EBFB120641A4703D8A142713414E6210584D120641AC1C5A64602713414E621058431106414E621058C9271341A01A2FDD1B11064121B07268D8271341F4FDD478571006411904560ECF271341448B6CE7DD0F06418195438BC1271341F6285C8F6A0F06413BDF4F0DA72713410E2DB29D360F06416F1283C091271341E5D022DB070F0641F6285C0F7227134160E5D022DA0E06416F1283404F271341448B6CE7D30E0641F853E3A52427134154E3A59BE60E06411904568E07271341643BDF4F0C0F0641D7A3703DFC2613414A0C022B100F064125068115FB26134191ED7C3F310F0641B6F3FDD4F42613414C378941F11006414A0C022BA0261341EC51B81ECC1206413BDF4F0D40261341022B87167514064125068115F1251341B4C876BE8C160641C74B37897F25134121B07268C6160641508D976E7525134123DBF97EE316064146B6F3FDD3251341'::geometry,28992)
,(ST_Buffer(ST_GeomFromText('LINESTRING(180955 313700,180920 313740)', 28992),0.1))),2)
)))

Results

To my surprise I got two spikes in SQL Server. To verify this, I looked into PostGIS. I got two spikes as well. First I thought it was caused by precision, going to WKT and back. But also if you use WKB, or if you use the geometry from the underlying database (not shown here), you will get the same. So it is the imprecision of the handling. When varying the cutline a bit (changing 180955 to 180960), I get one spike in SQL Server and still two in PostGIS.

The spikes are shown here, bordering the blue result:


The perimeter of the result is interesting. The result, without spikes, should be 92.5 meter. But including the spikes it is 151.792871568541 (in SQL Server, using the function .STLength() and 151.792880859545 (in PostGIS, using the function ST_Perimeter(...).

The area is everywhere about 295.4 square meter, with or without spikes, indicating that the spikes have an area about zero (so those are real spikes).


Why would you do such a thing... Creating a parcel, cutting, throwing away one, subtracting and getting the other one...

Well, I was just testing our application. Our application has WebGIS functionality to create parcels. And to cut parcels into two pieces. And to throw away pieces. And to fill-up space which is left from a parcel and a piece. So just this scenario. And to my annoyance, often you get spikes! This project is done with SQL Server, and we get spikes. But this little research shows that, would we have done it with PostGIS, we would have had resulting spikes as well.

It is as a simple mathematical exercise. We cut 5 into two pieces, 3 and 2. We throw away the 2 and keep the 3. We subtract 3 from our 5. And we get the 2 back, of course. In geometry we subtract (difference) twice but that is just because the cut is not an OGC operation. So we trick it.

Note that I didn't get this just once. I got it several times for different geometries. Note also that it can easily be solved using a STBuffer of -0.00001 on the resulting polygon, which will eliminate the spikes and have no visual impact (it has an impact though...).

With Boost.Geometry

Of course I was curious what Boost.Geometry would do with this.

The buffer of Boost.Geometry is not completely implemented, so I create the buffer in SQL Server, saving the WKB.

select geometry::STGeomFromText('LINESTRING(180955 313700,180920 313740)', 28992).STBuffer(0.1).STAsBinary()

and draft some C++ code:

template <typename T> 
void test_difference_parcel_precision() 
{
    namespace bg = boost::geometry; 
    typedef bg::model::d2::point_xy<T> point_type; 
    typedef bg::model::polygon<point_type> polygon_type;
    typedef std::vector<boost::uint8_t> byte_vector;

    polygon_type parcel, buffer; 

    { 
        byte_vector wkb; 
        bg::hex2wkb("0103000000010000001A00000023DBF97EE316064146B6F3FDD32513415A643BDFD216064175931804E225134196438B6C52150641C976BE9F34261341F6285C8F06150641022B871641261341F853E3A5D3140641E17A14AE50261341B81E85EBFB120641A4703D8A142713414E6210584D120641AC1C5A64602713414E621058431106414E621058C9271341A01A2FDD1B11064121B07268D8271341F4FDD478571006411904560ECF271341448B6CE7DD0F06418195438BC1271341F6285C8F6A0F06413BDF4F0DA72713410E2DB29D360F06416F1283C091271341E5D022DB070F0641F6285C0F7227134160E5D022DA0E06416F1283404F271341448B6CE7D30E0641F853E3A52427134154E3A59BE60E06411904568E07271341643BDF4F0C0F0641D7A3703DFC2613414A0C022B100F064125068115FB26134191ED7C3F310F0641B6F3FDD4F42613414C378941F11006414A0C022BA0261341EC51B81ECC1206413BDF4F0D40261341022B87167514064125068115F1251341B4C876BE8C160641C74B37897F25134121B07268C6160641508D976E7525134123DBF97EE316064146B6F3FDD3251341", std::back_inserter(wkb)); 
        bg::read_wkb(wkb.begin(), wkb.end(), parcel); 
    } 
    { 
        byte_vector wkb; 
        bg::hex2wkb(
"01030000000100000083000000000032FCD716064100009E998F25134100000706D81606410040A4998F2513410000DA0FD816064100C0E6998F2513410000A819D81606410080659A8F25134100806A23D816064100C0209B8F25134100801E2DD81606410080189C8F2513410000BE36D816064100404D9D8F25134100004540D816064100C0BE9E8F2513410000AF49D816064100806DA08F2513410000F752D8160641008059A28F2513410000195CD816064100C082A48F25134100800F65D81606410080E9A68F2513410000D66DD816064100408EA98F25134100006876D816064100C070AC8F2513410000C17ED8160641000091AF8F2513410080DC86D816064100C0EFB28F25134100009E8ED816064100C081B68F2513410080EC95D816064100803ABA8F2513410080C79CD8160641000018BE8F25134100002FA3D8160641008017C28F251341008022A9D8160641000037C68F2513410080A1AED8160641000074CA8F2513410000ACB3D81606410040CCCE8F251341000042B8D816064100403DD38F251341000062BCD81606410000C5D78F25134100000DC0D8160641000061DC8F251341000042C3D816064100000FE18F251341000001C6D81606410080CCE58F251341008049C8D8160641004097EA8F25134100001BCAD816064100006DEF8F251341008075CBD816064100804BF48F251341008058CCD8160641004030F98F2513410000C4CCD8160641000019FE8F2513410080B7CCD81606410080030390251341008032CCD81606410000ED0790251341000035CBD81606410000D40C902513410080BEC9D81606410040B511902513410000CFC7D816064100408F1690251341008065C5D816064100005F1B90251341008082C2D81606410080222090251341000025BFD81606410080D7249025134100004DBBD816064100807B29902513410080FAB6D816064100800C2E9025134100002DB2D816064100C08732902513410080E3ACD81606410000EB369025134100801EA7D81606410000343B902513410000DEA0D81606410080603F902513410080209AD816064100406E43902513410080209AC015064100406E43302613410080FC92C015064100004F473026134100008B8BC01506410040F64A302613410000D083C015064100C0634E302613410000D17BC0150641008097513026134100009273C015064100409154302613410000186BC015064100C050573026134100806762C01506410000D6593026134100808559C01506410000215C3026134100007650C01506410000315E3026134100003E47C015064100800660302613410000E23DC01506410000A1613026134100006734C015064100800063302613410080D12AC015064100C024643026134100002621C015064100800D653026134100006917C015064100C0BA653026134100809F0DC015064100402C66302613410000CE03C015064100006266302613410000F9F9BF15064100C05B6630261341000026F0BF1506410040196630261341000058E6BF15064100809A6530261341008095DCBF1506410040DF64302613410080E1D2BF1506410080E76330261341000042C9BF15064100C0B262302613410000BBBFBF1506410040416130261341000051B6BF1506410080925F30261341000009ADBF1506410080A65D302613410000E7A3BF15064100407D5B302613410080F09ABF150641008016593026134100002A92BF15064100C071563026134100009889BF15064100408F533026134100003F81BF15064100006F503026134100802379BF1506410040104D3026134100006271BF15064100407E49302613410080136ABF1506410080C5453026134100803863BF1506410000E841302613410000D15CBF1506410080E83D302613410080DD56BF1506410000C9393026134100805E51BF15064100008C35302613410000544CBF15064100C03331302613410000BE47BF15064100C0C22C3026134100009E43BF15064100003B28302613410000F33FBF15064100009F23302613410000BE3CBF1506410000F11E302613410000FF39BF1506410080331A302613410080B637BF15064100C06815302613410000E535BF150641000093103026134100808A34BF1506410080B40B302613410080A733BF15064100C0CF063026134100003C33BF1506410000E7013026134100804833BF1506410080FCFC2F2613410080CD33BF150641000013F82F2613410000CB34BF15064100002CF32F26134100804136BF15064100C04AEE2F26134100003138BF15064100C070E92F26134100809A3ABF1506410000A1E42F26134100807D3DBF1506410080DDDF2F2613410000DB40BF150641008028DB2F2613410000B344BF150641008084D62F26134100800549BF1506410080F3D12F2613410000D34DBF150641004078CD2F26134100801C53BF150641000015C92F2613410080E158BF1506410000CCC42F2613410000225FBF15064100809FC02F2613410080DF65BF15064100C091BC2F2613410080DF65D716064100C091BC8F2513410080036DD71606410000B1B88F25134100007574D716064100C009B58F2513410000307CD716064100409CB18F25134100002F84D7160641008068AE8F25134100006E8CD716064100C06EAB8F2513410000E894D71606410040AFA88F2513410080989DD716064100002AA68F25134100807AA6D71606410000DFA38F25134100008AAFD71606410000CFA18F2513410000C2B8D71606410080F99F8F25134100001EC2D716064100005F9E8F251341000099CBD71606410080FF9C8F25134100802ED5D71606410040DB9B8F2513410000DADED71606410080F29A8F251341000097E8D71606410040459A8F251341008060F2D716064100C0D3998F251341000032FCD716064100009E998F251341", std::back_inserter(wkb)); 
        bg::read_wkb(wkb.begin(), wkb.end(), buffer); 
    } 
    bg::correct(parcel); 
    bg::correct(buffer); 

    std::vector<polygon_type> pieces; 
    bg::difference(parcel, buffer, pieces); 

    std::vector<polygon_type> filled_out; 
    bg::difference(parcel, pieces.back(), filled_out); 

    std::cout << std::setprecision(16) << bg::wkt(filled_out.front()) << std::endl; 
    std::cout << bg::area(filled_out.front()) << std::endl; 
    std::cout << bg::perimeter(filled_out.front()) << std::endl; 
}

We are using WKB here to not loose precision. I run this with three different numeric types (float, double and ttmath high precision). And that is the essence of Boost.Geometry! Not one type! Many types! All types!

int main(int, char* [])
{ 
    test_difference_parcel_precision<float>(); 
    test_difference_parcel_precision<double>(); 
    test_difference_parcel_precision<ttmath_big>(); 
    return 0; 
}

Using these three types we get:

• float: one spike, perimeter 136.28
  
• double: two spikes, perimeter 151.793
  
• ttmath: zero spikes, correct result, perimeter 92.50548050416428971188559805054187252

Conclusions

• SQL Server is not perfect, it can cause spikey features
  
• PostGIS is not perfect, it can cause spikey features
  
• The cause is not WKT, nor it is WKB
  
• Boost.Geometry, using double, is not perfect because it causes the same spikey features
  
• Boost.Geometry, using ttmath, gives the correct results
  
• Boost.Geometry, using float, gives one spikey feature
  
• Changing the cutline a bit causes only one spikey feature within SQL Server, still two in PostGIS, and zero in Boost.Geometry (with all numeric types).
  
• Be aware of floating point precision

Free Shapefiles of Countries of the World (4)

There was the suggestion to look at gadm (Global Administrative Areas), thanks for this suggestion. This blog handles this new source in detail. Note that some sources are apparently hard to find, i.e. not found immediately by googling to "world country shapefiles" or similar. I'm glad that this blog appears now in the first 10 results, so hope to provide by this an index to world country shapefiles.

Look also at  part 1,  part 2, and part 3 of this blog.

Large

This newly found shapefile is very interesting. The first observation is that it is rather large, the zipfile is 158MB and takes several minutes of downloading. Unzipped the triplet (shp,shx,dbf) is almost half a gigabyte. When we zoom in, we see how this is caused, it is a shapefile taken from a grid:




The image above is a piece of the island Elba, the magenta line is the kml-ified shapefile displayed in Google Earth (yes, I like Google Earth).

I don't know if I like all those cornery borders... At least, it is not necessary to keep them, why are they not dropped or generalized away?

Anyway, if we zoom out a little again, the borders seem to be reasonably accurate at the first sight, here is the whole island Elba:

Vatican City and Liechtenstein

In last blog I described the small countries in and around Italy. So I shortly redo this for this new datasource. The gadm contains really the most accurate approximatino of Vatican City, compared with Google Earth as still our Single Point of Truth...



It is much better than all other sources found, and it therefore raises an interesting question: what is that little block where it differs, just to the southwest of Saint Peter's Square? Who is right there? And who is wrong? Interesting questions in the neighbourhood of the heart of the Roman Church... But we mean it geographically. Is Google failing here, as in Nicaragua and Costa Rica, nearly causing a war? Or is it the gadm.org website, a particular initiative providing freely available shapefiles for academic and non-commercial use (e.g. blogging)? According to Wikipedia, there indeed is a border dispute near Vatican City, but that is a little more to the east...

And look, we see a similar issue with Liechtenstein, probably a disputed valley...



But for the rest also here this shapefile is relatively accurate. For its size, accuracty is expected, of course...

Monaco

In last blog we saw that the borders of Monaco were correct in only one of the compared shapefiles... Alas, wrong again here:



It seems hard to create a correct shapefile... What is the source of all those errors? Who is really right here? Of course Google! For all certainty I checked this at Wikipedia (because I've never been in Monaco). Yes, the yellow line of Google is following Monaco's portrayal within the Wiki page of Monaco. So alas, gadm is totally wrong here.

The Netherlands

I'm living there, in that little country which is so difficult to be stored in free shapefiles... There is really a large lake, only bordered from the sea by a dike. That lake should be displayed, and it is not. If very small and uninhabited islands are displayed, the lake should really be there. And it is not the only error, in Zeeland the (smaller) water bodies behind the dikes are omitted as well.

The Caspian Sea

This was another requirement (CS1). It is not as bad as it was in some other shapefiles, where the Caspian Sea was divided between surrounding countries. But the Caspian Sea is included as a separate country in gadm... Interesting. Do they issue passports?

Validity

Finally there is validity. All countries should be modelled as multi-polygons and be valid in both SQL Server and PostGIS. We've already seen that this requirement is met in only one caes, until now (11). And also for this case, 13, alas, there are many multi-polygons invalid in SQL Server. In PostGIS it is OK, though the query runs for a long time due to the excessive size of  the geometries. The values are post-editted in the previous blog.

Actually it is nice to show a variant of the query, using the with syntax:

with valid as (select ST_IsValid(the_geom) as v from world13)
select v,count(1) from valid group by v;

Conclusion

As always I'm glad with new shapefiles. This one is quite accurate, and in many aspects the best until now. But also, it is way too large, too cornery, and contains errors, as all others do...

So I think a good and freely available shapefile of the countries of the world is still very welcome.

If there are more shapefiles, I'm still interested, but especially I'm interested in a the Google country borders in any format? Are they available?