[gdal-dev] Clipping shapefile with another produces invalid shapes

[Paul Meems]

Hi list,

I have a very large shapefile, with almost 700k shapes which I want to clip
with a border shapefile (in red):
https://ibb.co/gKth1G

I'm using this command:
ogr2ogr -clipsrc "border.shp" -overwrite -explodecollections -f "ESRI
Shapefile" ogr_clipped.shp Fishnet.shp
with GDAL 2.1.3, released 2017/20/01
Eventually, I'll transform the command to a call to ogr.VectorTranslate().

The result is visually correct:
https://ibb.co/hOvJab

But at the red asterix, I have an invalid shape. The reason is according to
GEOS:
Polygon must be clockwise
I have more of these errors.
They all seem to occur when the shape due to the clipping is triangular
shaped:
https://ibb.co/cKuY8w
https://ibb.co/eXcWvb

Is this an error/warning which can be expected or do I need to add an
additional parameter?
If needed I can try to make a test set with as little data as possible.

Thanks,

Paul

[Even Rouault]

Paul,

Post by Paul Meems
But at the red asterix, I have an invalid shape. The reason is according to
Polygon must be clockwise

I can't see such error emitted by GEOS source code.

***@even-i700:~$ grep -r "must be clockwise" geos/
***@even-i700:~$


Are you sure this comes from it ? As far as I know, GDAL and GEOS are pretty much
indifferent regarding the winding order of polygons. Other tools might be more sensitive.

You could possibly use the ST_ForceLHR() function of spatialite (through the SQL SQLite
dialect) to post process your result:
http://www.gaia-gis.it/gaia-sins/spatialite-sql-4.2.0.html

Even

[Paul Meems]

Thanks Even for your quick reply.

I was using MapWinGIS to perform the validation and I assumed it was just
using GEOS.
I did a closer look at the code and it is first doing some internal checks
like the minimum number of points and the order of the points.
The code is checking if the polygon shape only has 1 part. If so it should
be clockwise.
After that, it is using the GEOS:IsValid method.
As I understand you correctly I should interpret this message as a warning,
not as an error.

I also tried your suggestion:
ogr2ogr -sql "select ST_ForceLHR(geometry) from ogr_clipped" -dialect
sqlite ogr_clipped_lhr.shp ogr_clipped.shp
It is very fast, only a few seconds on my 700k shapes, but the result is
the same. I still get the 'Polygon must be clockwise' error.
Of course, MapWinGIS also has a fix method, but it is very strict and not
so fast with 700k shapes.

I'll send the resulting shapefile to my client and ask if his application
can process the data.

Thanks,



Paul

*Paul Meems *
Release manager, configuration manager
and forum moderator of MapWindow GIS.
www.mapwindow.org

Owner of MapWindow.nl - Support for
Dutch speaking users.
www.mapwindow.nl


*The MapWindow GIS project has moved to GitHub
<https://github.com/MapWindow>!*


Download the latest MapWinGIS mapping engine.
<https://github.com/MapWindow/MapWinGIS/releases>

Download the latest MapWindow 5 open source desktop application.
<https://github.com/MapWindow/MapWindow5/releases>

Post by Even Rouault
Paul,

Post by Paul Meems
But at the red asterix, I have an invalid shape. The reason is according

to

Post by Paul Meems
Polygon must be clockwise

I can't see such error emitted by GEOS source code.
Are you sure this comes from it ? As far as I know, GDAL and GEOS are
pretty much indifferent regarding the winding order of polygons. Other
tools might be more sensitive.
You could possibly use the ST_ForceLHR() function of spatialite (through
http://www.gaia-gis.it/gaia-sins/spatialite-sql-4.2.0.html
Even
--
Spatialys - Geospatial professional services
http://www.spatialys.com

...

[Even Rouault]

Post by Paul Meems
Thanks Even for your quick reply.
I was using MapWinGIS to perform the validation and I assumed it was just
using GEOS.
I did a closer look at the code and it is first doing some internal checks
like the minimum number of points and the order of the points.
The code is checking if the polygon shape only has 1 part. If so it should
be clockwise.
After that, it is using the GEOS:IsValid method.
As I understand you correctly I should interpret this message as a warning,
not as an error.

Different software, standard, formats have differen rules on if there must be a winding order
and which one it is...

A return of false by GEOS:IsValid() should be examined a bit more closely (but apparently you
get this error before calling IsValid())

Post by Paul Meems
ogr2ogr -sql "select ST_ForceLHR(geometry) from ogr_clipped" -dialect
sqlite ogr_clipped_lhr.shp ogr_clipped.shp

Jukka answered on the LHR vs RHR confusion, but I completely missed the fact that you were
outputting to shapefile. In that case, whatever the original orientation, or if you modified it
with a SQL statement, the OGR shapefile writer will force the polygons to follow the
shapefile specification, which requires the outer ring to be in clockwise order. So you
shouldn't get this error... Which means that either there's an issue in some case in the
shapefile driver (numerical precision issue), or in the check in MapWinGIS.

Even

[Paul Meems]

Thanks Even and Jukka for the clarification.

Due to historical reasons, MapWinGIS is not using GDAL/OGR to read/write
shapefiles. It is using its own implementation.
For now, I will clip with OGR first and then fix with MapWinGIS, this will
produce valid shapes.
It does take a bit longer.

Should I create a test dataset to reproduce this issue. I'm not capable of
debugging GDAL myself.




Paul

*Paul Meems *
Release manager, configuration manager
and forum moderator of MapWindow GIS.
www.mapwindow.org

Owner of MapWindow.nl - Support for
Dutch speaking users.
www.mapwindow.nl


*The MapWindow GIS project has moved to GitHub
<https://github.com/MapWindow>!*


Download the latest MapWinGIS mapping engine.
<https://github.com/MapWindow/MapWinGIS/releases>

Download the latest MapWindow 5 open source desktop application.
<https://github.com/MapWindow/MapWindow5/releases>

Post by Paul Meems

Post by Paul Meems
Thanks Even for your quick reply.
I was using MapWinGIS to perform the validation and I assumed it was just
using GEOS.
I did a closer look at the code and it is first doing some internal

checks

Post by Paul Meems
like the minimum number of points and the order of the points.
The code is checking if the polygon shape only has 1 part. If so it

should

Post by Paul Meems
be clockwise.
After that, it is using the GEOS:IsValid method.
As I understand you correctly I should interpret this message as a

warning,

Post by Paul Meems
not as an error.

Different software, standard, formats have differen rules on if there must
be a winding order and which one it is...
A return of false by GEOS:IsValid() should be examined a bit more closely
(but apparently you get this error before calling IsValid())

Post by Paul Meems
ogr2ogr -sql "select ST_ForceLHR(geometry) from ogr_clipped" -dialect
sqlite ogr_clipped_lhr.shp ogr_clipped.shp

Jukka answered on the LHR vs RHR confusion, but I completely missed the
fact that you were outputting to shapefile. In that case, whatever the
original orientation, or if you modified it with a SQL statement, the OGR
shapefile writer will force the polygons to follow the shapefile
specification, which requires the outer ring to be in clockwise order. So
you shouldn't get this error... Which means that either there's an issue in
some case in the shapefile driver (numerical precision issue), or in the
check in MapWinGIS.
Even
--
Spatialys - Geospatial professional services
http://www.spatialys.com

...

[jratike80]

Post by Even Rouault
You could possibly use the ST_ForceLHR() function of spatialite (through the SQL SQLite
http://www.gaia-gis.it/gaia-sins/spatialite-sql-4.2.0.html

Interesting, that ST_ForceLHR in Spatialite gives the same result than
ST_ForceRHR in PostGIS.

Description for Spatialite:
Any Polygon Ring will be oriented accordingly to Left Hand Rule (Exterior
Ring will be clockwise oriented, and Interior Rings will be
counter-clockwise oriented).

Description for PostGIS:
Forces the orientation of the vertices in a polygon to follow a
Right-Hand-Rule, in which the area that is bounded by the polygon is to the
right of the boundary. In particular, the exterior ring is orientated in a
clockwise direction and the interior rings in a counter-clockwise direction.
This function is a synonym for ST_ForcePolygonCW

Perhaps the rule of thumb would be not to talk just about right hand/left
hand rules in GIS because they are ambiguous, but mention also what it
really means like in the descriptions above, or in the GeoJSON
specification:

A linear ring MUST follow the right-hand rule with respect to the
area it bounds, i.e., exterior rings are counterclockwise, and
holes are clockwise.

-Jukka Rahkonen-




--
Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html