overlaps->intersections (#283)

import shapely

import sqlalchemy

def compute_overlaps(sql='(SELECT * FROM images) AS images'):

    """

    For the candidate graph, compute the overlapping polygons that

    comprise the entire candidate graph / footprint map. Each overlap

    includes an 'overlaps' attribute/column that includes a list of the

    footprint polygons that have contributed to given overlap.

    """

    query = f"""

SELECT ST_AsEWKB(geom) AS geom FROM ST_Dump((

    SELECT ST_Polygonize(the_geom) AS the_geom FROM (

        SELECT ST_Union(the_geom) AS the_geom FROM (

            SELECT ST_ExteriorRing((ST_DUMP(footprint_latlon)).geom) AS the_geom

            FROM {sql}) AS lines

    ) AS noded_lines

)

)"""

    if not Session:

        warnings.warn('This function requires a database connection configured via an autocnet config file.')

        return 

    session = Session()

    oquery = session.query(Overlay)

    iquery = session.query(Images)

    srid = config['spatial']['srid']

    rows = []

    for q in engine.execute(query).fetchall():

        overlaps = []

        b = bytes(q['geom'])

        qgeom = shapely.wkb.loads(b)

        res = iquery.filter(Images.footprint_latlon.ST_Intersects(geoalchemy2.shape.from_shape(qgeom,

                                                                                               srid=srid)))

        for i in res:

            fgeom = geoalchemy2.shape.to_shape(i.footprint_latlon)

            area = qgeom.intersection(fgeom).area

            if area < 1e-6:

                continue

            overlaps.append(i.id)

        o = Overlay(geom=f'srid={srid};{qgeom.wkt}', overlaps=overlaps)

        res = oquery.filter(Overlay.overlaps == o.overlaps).first()

        if res is None:

            rows.append(o)

    session.bulk_save_objects(rows)

    session.commit()

    # If an overlap has only 1 entry, it is a sliver and we want to remove it.

    res = oquery.filter(sqlalchemy.func.array_length(Overlay.overlaps, 1) <= 1)

    res.delete(synchronize_session=False)

    session.commit()

    session.close()

def place_points_in_overlaps(cg, size_threshold=0.0007, reference=None, height=0,

                             iterative_phase_kwargs={'size':71}):

    lla = pyproj.Proj(proj='latlon', a=semi_major, b=semi_minor)   

    # TODO: This should be a passable query where we can subset.

    for o in session.query(Overlay.id, Overlay.geom, Overlay.overlaps).\

    for o in session.query(Overlay.id, Overlay.geom, Overlay.intersections).\

             filter(sqlalchemy.func.ST_Area(Overlay.geom) >= size_threshold):

        valid = compgeom.distribute_points_in_geom(geoalchemy2.shape.to_shape(o.geom))

        if not valid:

            continue

        overlaps = o.overlaps

        overlaps = o.intersections

        if reference is None:

            source = overlaps[0]