Thursday, 6 July 2017

coordinate system - Preventing ogr2ogr from generating invalid geometries in transformation


I have the following GeoPackage: https://drive.google.com/open?id=1COcGcbZQIsW6NrzuD-K-Kijpvg4BYQ54


1


I have checked it with spatialite and it does not have invalid geometries.


It has a dissolved, promoted to single parts and simplified derived layer.


I want to transform it to a custom reference system via ogr2ogr. For that, I have created a batch file named transform.bat with the following lines:


ogr2ogr -skipfailures^
-ct "+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=0 +y=0 +z=0 +rx=0 +ry=109110.996 +rz=-215775 +s=0 +convention=coordinate_frame +step +inv +proj=cart +R=6371007 +step +proj=pop +v_3 +step +proj=longlat +R=6371007 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1"^
-overwrite^

output_dts.gpkg input_dts.gpkg^
-nlt MULTIPOLYGON^
-clipdst -179.99 -89.99 179.99 89.99^
-wrapdateline

The transformed reference system is based on an spheric datum, with a Helmert transformation to rotate it around Y and Z geocentric axes.


When I run it, seems to detect invalid geometries:


C:\GA\GIS\Pruebas>transform.bat

C:\GA\GIS\Pruebas>ogr2ogr -skipfailures -ct "+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=0 +y=0 +z=0 +rx=0 +ry=109110.996 +rz=-215775 +s=0 +convention=coordinate_frame +step +inv +proj=cart +R=6371007 +step +proj=pop +v_3 +step +proj=longlat +R=6371007 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1" -overwrite output_dts.gpkg input_dts.gpkg -nlt MULTIPOLYGON -clipdst -179.99 -89.99 179.99 89.99 -wrapdateline

ERROR 1: TopologyException: Input geom 0 is invalid: Self-intersection at or near point 178.85626517687388 -44.274345298926086 at 178.85626517687388 -44.274345298926086
ERROR 1: TopologyException: Input geom 0 is invalid: Self-intersection at or near point 133.29913840320728 72.555284217735704 at 133.29913840320728 72.555284217735704



The transformation is well done, but some geometries were lost:


2




The problem is with that geometries that crosses the dateline when transformed. The clipping boundary and the promote to multipolygon were my attemps to return the geometries clipped and valid.


If I remove the -nlt, -clipdst and -wrapdateline parameters, the Americas and Antarctica geometries are valid but cross the dateline:


3





I have found that geometries must have a kind of special shape, in a special place, to turn invalid. This behavior can be reproduced with the following geometry:


POLYGON (( -130 60, -120 50, -125 40, -130 45, -135 40, -145 50, -140 60, -135 55, -130 60 ))




I was able to fix by hand the transformed geometries, by creating new nodes outside the extents and clipping the layer.


4


But I can't figure how to solve this kind of problems via ogr2ogr.




No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...