[gdal-dev] Resulting geometry from OGRGeometry::difference somehow got simplified

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|

[gdal-dev] Resulting geometry from OGRGeometry::difference somehow got simplified

mikeucfl
I need some help trying to make sense of an issue I'm running into using
OGRGeometry::difference. I have some pretty large geometries (spans across
multiple longitude/latitude values) in which I am utilizing the difference
operation in OGRGeometry to cut holes in and get the resulting geometry.
This has worked for some time, but one of my datasets returns a result that
looks as if the geometry got heavily simplified. I'm not a great artist but
here's a drawn example of what I'm talking about:
<http://osgeo-org.1560.x6.nabble.com/file/t381880/5U0bWIg.png>
(the black is what the resulting geometry should look like, and the red is
what it actually looks like).

I've done testing by trying to reduce the size and or number of interior
rings for the geometry that is being used in the difference operation, which
seems to work but I'm not 100% convinced that those variables are the issue.
Does anyone know of any limitations with geometry sizes, interior ring
sizes, or possibly vert counts with GDAL/GEOS? Or any other ideas as to what
my issue could be? My software is written in c++ using Visual Studio 2015
(x64) and I'm linking to GDAL v2.1.3 with GEOS 3.5.0 (have also tried GDAL
2.0.1 with same results). Thanks.



--
Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Resulting geometry from OGRGeometry::difference somehow got simplified

mikeucfl
Just to add some more information.. upon some further inspection it seems all
the values are snapped to a grid. All my geometry is in a WGS84 (EPSG:4326)
projection, and the values are snapped to a 0.01 lat/long grid. I don't do
anything to set a tolerance or anything for these values, or see where I
could set these from the GDAL layer (as the driver I'm using wraps most of
the functionality for the FGDB database using the ESRI FileGDB driver)



--
Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Resulting geometry from OGRGeometry::difference somehow got simplified

mikeucfl
Got permission to show the actual data:
<http://osgeo-org.1560.x6.nabble.com/file/t381880/RUklfeY.png>  (left side
is what the data should look like, red is what the data looks like after
calling the difference operation).



--
Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Resulting geometry from OGRGeometry::difference somehow got simplified

Even Rouault-2
In reply to this post by mikeucfl

On lundi 15 janvier 2018 16:49:57 CET mikeucfl wrote:

> Just to add some more information.. upon some further inspection it seems

> all the values are snapped to a grid. All my geometry is in a WGS84

> (EPSG:4326) projection, and the values are snapped to a 0.01 lat/long grid.

> I don't do anything to set a tolerance or anything for these values, or see

> where I could set these from the GDAL layer (as the driver I'm using wraps

> most of the functionality for the FGDB database using the ESRI FileGDB

> driver)

 

The FileGDB format using integer storage with a precision for coordinates. Perhaps this is a missconfiguration of this precision when you have created your FGDB database ?

I guess you would have the same issue if you just read or ogr2ogr the source layer without doing any processing on it.

 

Even

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Resulting geometry from OGRGeometry::difference somehow got simplified

mikeucfl
So the geometry before using the difference operation looks good (geodetic
values go up to 8-9 decimal places), and then the resulting geometry is
truncated to two (this is before placing in any file geodatabase) as well. I
have checked the source feature classes and verified that the precision is
set to 0.00000001, so looks good to me there.

One odd thing I noticed too is I output the two problematic geometries to a
new file geodatabase so I could iterate on the problem quicker, but when I
did that it tells me one of the geometries is invalid and the resulting
difference geometry is NULL. I had to do a geometry->buffer(0.0) to repair
it, and got valid geometry coming out.. but when I run it from the original
source all the geometries say they're valid and even forcing a buffer(0.0)
doesn't seem to fix it.

Really odd problem, not sure what's going on.. especially since my small
test case doesn't seem to re-create what the large dataset does. I'm
wondering if writing the geometry out is not the exact same as it was in
memory and possibly modifying the values? That would make some sense as to
why I see different results.. but again, kinda lost on what's going on.



--
Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Resulting geometry from OGRGeometry::difference somehow got simplified

mikeucfl
Just to update with what I've figured out so far.

1. The issue surrounds an OGRMultiGeometry that has intersecting interior
rings.
2. The projection of the geometries are in EPSG 4326 (lat/long).
3. Internally in GEOS when the difference operator runs, there is a call to
check if the geometry is valid and it threw an exception saying it wasn't
and spit out a warning to the console that I didn't notice. Geos tried to
fix the geometry and in doing so the values got truncated or rounded to the
hundredths place. I'm guessing some algorithm is thinking the projection is
in meters and is what causes my data to get screwed up.

I have tried fixing the geometry prior to calling difference with a
OGRGeometry::Buffer(0.0) but that still doesn't fix the issue. What seemed
to resolve it was putting the geometry in a different projection that is in
meters (used web mercator), calling Buffer(0.0) and then calling Difference.

So if anyone sees a similar issue, it's something internally with GEOS but I
didn't have enough time to track it down. Re-projecting to something using
meters does help though.



--
Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Resulting geometry from OGRGeometry::difference somehow got simplified

jratike80
mikeucfl wrote
> So if anyone sees a similar issue, it's something internally with GEOS but
> I
> didn't have enough time to track it down. Re-projecting to something using
> meters does help though.

Hi,

If you could share couple of problematic geometries for example as WKT then
somebody else could have a try as well.

-Jukka Rahkonen-






--
Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev