test_geos_sfgal_limits.py

import struct
from fractions import Fraction
import psycopg2

def double_to_hex(f):
    h = hex(struct.unpack('<Q', struct.pack('<d', f))[0])[2:]
    # reverse each byte for endianness
    return h[6:]+h[4:6]+h[2:4]+h[0:2]

def poly_wkb(x,y):
    #polygon((0 0,3 1,1 1,1 0.33333333,0 0))
    poly_wkb_a = "01030000000100000005000000000000000000000000000000000000000000000000000840000000000000F03F000000000000F03F000000000000F03F"
    #000000000000F03F"
    #DA12C1515555D53F
    poly_wkb_b = "00000000000000000000000000000000"
    return poly_wkb_a + double_to_hex(x) + double_to_hex(y) + poly_wkb_b

if __name__ == "__main__":
    conn = psycopg2.connect("dbname=test_postgis")
    cur = conn.cursor()

    cur.execute("select * from no_plan()")
    old_wkb = None
    for k in range(1000):
        x = Fraction(1,1) + k*Fraction(1,2**53)
        y = x * Fraction(1,3)
        fx = float(x)
        fy = float(y)
        wkb=poly_wkb(fx,fy)
        is_valid = fy < y
        if old_wkb is not None and old_wkb == wkb:
            continue

        cur.execute("select st_isvalid('{}')".format(wkb))
        geos_is_valid = cur.fetchone()[0]

        cur.execute("prepare t as select st_extrude('{}',0,0,1)".format(wkb))
        cur.execute("select substring(lives_ok('t'),1,2)='ok'")
        sfcgal_is_valid = cur.fetchone()[0]
        cur.execute("deallocate t")

        print is_valid, geos_is_valid, sfcgal_is_valid
        if sfcgal_is_valid != geos_is_valid:
            print "!!!!!!!!!!!!!!!!!!!!!!!!"
        old_wkb = wkb