[Dachs-support] Spherepoly line segment overlap

[Markus Demleitner]

Hi Sebastian,

On Wed, Jan 03, 2018 at 06:03:47PM +0100, Sebastian Walter wrote:
> The polygons come unaltered from the data's labels, and are stored in a PG
> database and PostGIS. I made an experiment with the postGIS ST_simplify
> function (a Douglas-Peucker implementation), which reduces the coordinates
> of our example to the following:
> 160.186005 44.483501 158.231003 44.4636 158.253998 44.465 158.395996
> 41.412899 158.539993 36.9044 160.156006 36.912498 160.289993 44.473301
> 160.186005 44.483501
> 
> Here's what I get out from gavo.utils.pgsphere.SPoly.asPgSphere for this
> simplified geometry:
> spoly '{(2.7957732029,0.7763835553),(2.7616519811,0.7760362173),(2.7620533195,0.7760606519),(2.7645316522,0.7227914403),(2.7670448739,0.6441032885),(2.7952496215,0.6442446252),(2.7975881358,0.7762055317),(2.7957732029,0.7763835553)}';

That looks much better (and I'd say at least for discovery purposes
any precision possibly lost is definitely worth the simplification in
inspection, display, and computation).  However, it seems the "odd"
points are still there (see attached PDF).

To investigate this a bit closer, I've written a quick script to
compute the steps between the points:

  points = [
  (2.7957732029,0.7763835553),
  (2.7616519811,0.7760362173),
  (2.7620533195,0.7760606519),
  (2.7645316522,0.7227914403),
  (2.7670448739,0.6441032885),
  (2.7952496215, 0.6442446252),
  (2.7975881358,0.7762055317),
  (2.7957732029,0.7763835553)]
for i, (p1, p2) in enumerate(zip(points[:-1], points[1:])):
	print i, p1[0]-p2[0], p1[1]-p2[1], \
		math.sqrt((p1[0]-p2[0])**2+(p1[1]-p2[1])**2), \
		math.atan2(p1[1]-p2[1], p1[0]-p2[0])

The output of this is

index  dx        dy          length         direction

0 0.0341212218 0.000347338 0.034122989623 0.0101791783835
1 -0.0004013384 -2.44345999999e-05 0.000402081535253 -3.08078492555
2 -0.0024783327 0.0532692116 0.0533268322466 1.617287476
3 -0.0025132217 0.0786881518 0.0787282764768 1.60272448352
4 -0.0282047476 -0.0001413367 0.0282051017237 -3.13658159946
5 -0.0023385143 -0.1319609065 0.13198162559 -1.58851573787
6 0.0018149329 -0.0001780236 0.00182364301157 -0.0977754789759


As you can see, from 0 to 1 the path essentially reverses its
direction, which is bound to confuse pgsphere.

A full diagnosis of such a situation is probably complex, but as a
simple approximation I tried the following piece of python that
throws out points that are "too close" to each other in some sense:

newpoints = [points[0]]
for p2 in points[1:]:
	p1 = newpoints[-1]
	if math.sqrt((p1[0]-p2[0])**2+(p1[1]-p2[1])**2)>0.01:
		newpoints.append(p2)

The resulting polygon is pgsphere-clean:

select spoly '{(2.7957732029, 0.7763835553), (2.7616519811, 0.7760362173), (2.7645316522, 0.7227914403), (2.7670448739, 0.6441032885), (2.7952496215, 0.6442446252), (2.7975881358, 0.7762055317)}';

If something like this can work for you, you can put the cleanup code
right into your RD.  As a child of the rowmaker, you can put an apply
like this:

  <apply name="simplify_polygon">
    <doc>
      An apply that takes polygon vertex coordinates from @poly,
      throws out vertices too close to each other, and makes
      an spoly in s_region out of the result.

      We needed this because pgsphere would barf on the polygons that
      came out of postgis -- they had odd kinks in them.
    </doc>
    <setup>
      <par name="min_dist" description="Minimal distance two
        points must have to be retained in the polygon">0.5</par>
      <code>
        from gavo.base.coords import getGCDist
      </code>
    </setup>
    <code>
      # assuming data comes in as [(long, lat), ...] in degrees in "poly"
      # change as necessary.
      points = @poly 
      dali_poly = [points[0][0], points[1][1]]
      for p in points[1:]:
        if getGCDist((dali_poly[-1][0], dali_poly[-1][1]), p)>min_dist:
      		dali_poly.extend(p)
      @s_region = pgsphere.SPoly.fromDALI(newpoints)
    </code>
  </apply>

What's going on here?

apply elements are procs, a general thing in DaCHS to stick python
code into RDs (other types of procs include rowfilters and
sourceFields in grammars, iterators in embedded grammars, regTests
and almost everything to do with datalink).  All of these can have
setup elements, which contain various sorts of initialisation stuff.
Here, we're making explicit that we have a parameter we're pulling
from thin air; do yourself and the folks that come after you a favour
and add a <doc> element to the apply and a description attribute to
the par.

I'm also pulling getGCDist ("get great circle distance") from the
bowels of DaCHS in the setup.  I'm sure there's an analog in astropy
or somewhere else, so you could have taken that from there as well.
Whatever names are defined in setup are available as python variables
in the proc code, as are the pars.

Within <code>, there's simple, generic python code.  Indentation
matters, and DaCHS takes the indentation of the first like off all
the following ones.  In case of doubt, just start the code in column
1.

The code itself is fairly generic, except that we're producing the
polygon in DALI style; DALI is a VO standard containing general
conventions, and one is that you can write a polygon as a flat float
array.  Which is what I'm doing here.

This code is untested and written off the top of my head.  Take with
a grain of salt and shout if you get incomprehensible gibberish from
DaCHS.

          -- Markus
-------------- next part --------------
A non-text attachment was scrubbed...
Name: poly2.pdf
Type: application/pdf
Size: 2814 bytes
Desc: not available
URL: <http://lists.g-vo.org/pipermail/dachs-support/attachments/20180104/a36c583e/attachment-0002.pdf>