[postgis-users] Image boundary cross-secting
Guido Lemoine
guido.lemoine at jrc.it
Thu Oct 27 07:05:50 PDT 2005
Dear List,
I am using postgres 8.0.3 with postgis 1.0.3 and GEOS 2.1.3. My quest
is to find fractal areas that are covered by one or more images out of
a monthly time series (of 44 images). Each image can have different size
and location in an area that is much larger than an individual image scene.
An fractal area can be covered by 0, 1 or more than one image. An image
scene
is neither rectangular, nor square in comparison to a geographical grid
(each one is actually a SAR image frame from either RADARSAT ScanSAR
narrow B or standard image mode or from ENVISAT ASAR image mode.
All images are from descending orbits. Image GCPs are taken, using
gdalinfo,
from the headers and converted to RHR polygons in the data base. All have
SRID=4326).
I am trying to use the intersection and difference functions to check
for geometric
overlap for pairs of images. If two images A and B partially overlap, I
store
the intersection and replace A with difference(A, intersection(A,B)) and
replace B
with difference(B, intersection(A,B)). All is done inside two loops in a
PLPGSQL
routine. The loops exit when no more intersections are found. I have
tested the
routine on a simple set of rectangular polygons that are varyiously
overlapping
(including full enclosures), and all goes well (and sufficiently fast as
well).
However, if I do the same for my "real world" image polygons, I get various
GEOS Intersection errors, such as:
TopologyException: no outgoing dirEdge found (some_lon, some_lat)
TopologyException: side location conflict (some_lon, some_lat)
(some_lon, some_lat = a coordinate pair that varies with polygon pair
that causes the
error).
These errors typically occur where one member of the pair is a newly
created
difference(A, intersection(A,B)) polygon (cf. the routine above). I have
inspected
these cases with JUMP and it turns out that there are some spurious
single point
vertices that form small lines at the crossing points of the cross
sections.
I googled for clues on this. The geometries usually pass the isValid
constraint but
their use in a new intersection throws the exceptions. Apart from this,
there is
little I can find.
Is there a way to avoid generating the spurious vertices or clean up for
them (they
usually compose a sub-segment of the polygon that has zero area).
Any hints welcome, as I have lost some serious time trying to solve it
(I even reverted
to using Java2D constructive geometry, which runs into some other
trouble, but I
may have to turn back to that).
I am willing to post further details, including routine, scripts and
offending sets of
polygons, if that helps,
Guido Lemoine
More information about the postgis-users
mailing list