Skip to content

Instantly share code, notes, and snippets.

@arbakker
Last active January 21, 2023 23:43
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save arbakker/47bbee69ce4194a41e32593731047f3e to your computer and use it in GitHub Desktop.
Save arbakker/47bbee69ce4194a41e32593731047f3e to your computer and use it in GitHub Desktop.
Spatialite ST_IsPolygonCCW Issue

Spatialite ST_IsPolygonCCW Issue

To reproduce the issue use this GeoPackage containing 6 geometries (valid) with a counterclockwise (CCW) winding order.

Description: When determining polygon winding order for the set of geometries in the provided GeoPackage the function ST_IsPolygonCCW returns 0 for all geometries while it should return 1. The function ST_IsPolygonCW does produce the correct result

The winding order of the polygons can be visualized with the following QGIS Layer definition file, which shows the winding order of the polygons is in fact CCW (warning: this QGIS style uses a geometry generator, which is quite computionally heavy when viewing larger geometries):

image

See the below steps to reproduce the issue:

Run in Bash shell:

wget https://www.dropbox.com/s/h2tf7vob53wpib3/ccw-issue-spatialite.gpkg
sqlite3 ccw-issue-spatialite.gpkg

Run following SQL commands in sqlite shell:

select load_extension('mod_spatialite');
select 'spatialite_version', spatialite_version();
select 'geos_version', geos_version();
select 'count ST_IsValid', count(*) from protected_site where  ST_IsValid(GeomFromGPB(geom)) = 1;
select 'count ST_IsPolygonCCW = 0', count(*) from  protected_site where ST_IsPolygonCCW(GeomFromGPB(geom)) == 0;
select 'count ST_IsPolygonCCW = 1', count(*) from  protected_site where ST_IsPolygonCCW(GeomFromGPB(geom)) == 1;
select 'count ST_IsPolygonCW = 0', count(*) from  protected_site where ST_IsPolygonCW(GeomFromGPB(geom)) == 0;
select 'count ST_IsPolygonCC = 1',  count(*) from  protected_site where ST_IsPolygonCW(GeomFromGPB(geom)) == 1;

This will output the following result:

spatialite_version|5.0.1
geos_version|3.10.2-CAPI-1.16.0
count ST_IsValid|6
count ST_IsPolygonCCW = 0|6
count ST_IsPolygonCCW = 1|0
count ST_IsPolygonCW = 0|6
count ST_IsPolygonCCW = 1|0

Addendum 22-01-23

Forcing the geoms a CCW orientation does not work using ST_ForcePolygonCCW:

update protected_site set geom = AsGPB(ST_ForcePolygonCCW(GeomFromGPB(geom)));

Rerunning select count(*) from protected_site where ST_IsPolygonCCW(GeomFromGPB(geom)) == 1; still returns a count of 0.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment