fixes to Cross Instrument Matcher (#407)

import shapely.geometry

from shapely.geometry import Polygon, Point

from shapely.affinity import scale

from shapely import wkt

from autocnet.utils import utils

from autocnet import Session

def two_point_extrapolate(x, xs, ys):

    """

    Parameters

    ----------

    x : float

             point where you want corresponding y value

    xs : ndarray

             (1, 2) array of point x coordinates

    ys : ndarray

             (1, 2) array of point x coordinates

    Returns

    -------

    y : float

            extrapolated value associated with x

    """

    m = (ys[1]-ys[0])/(xs[1]-xs[0])

    y = ys[0] + m*(x-xs[0])

    return x, y

def convex_hull_ratio(points, ideal_area):

    """

    """

    return np.argmin(np.sum((search - pt)**2, axis=1))

def find_side(side):

    """

    Parameters

    ----------

    side: str

            describes which extrema you cube you want; can equal 'east' or 'west'

    Returns

    -------

    lon : float

            longitude

    lat : float

            latitude

    """

    side = side.lower()

    func = {'east': 'st_xmax', 'west': 'st_xmin'}

    func = func[side]

    order = {'east': 'desc', 'west': 'asc'}

    order = order[side]

    query = f"""

    select ST_AsText(geom) from images

    order by {func}(geom) {order}

    limit 1 """

    session = Session()

    geom = session.execute(query).first()

    geom = wkt.loads(geom[0])

    session.close()

    #find eastern/wertern most side of east_geom/west_geom

    fp = geom.minimum_rotated_rectangle

    coords = np.column_stack(fp.exterior.xy)

    fp_lon, fp_lat = zip(*coords)

    # always a counter clockwise motion so find minimum/maximum lon index

    # and use i and i+1 lat lons points as return value

    if side == 'east':

        i = np.argmax(fp_lon)

        lon = fp_lon[i:i+2]

        lat = fp_lat[i:i+2]

    elif side == 'west':

        i = np.argmin(fp_lon)

        lon = fp_lon[i:i+2]

        lat = fp_lat[i:i+2]

    return np.array(lon), np.array(lat)

def create_points_along_line(p1, p2, npts):

    """

    Compute a set of nodes equally spaced between

    """

    return geom.contains(Point(x, y))

def distribute_points(geom, nspts, ewpts):

def distribute_points_classic(geom, nspts, ewpts):

    """

    This is a decision tree that attempts to perform a

    very simplistic approximation of the shape

    valid = [p for p in points if xy_in_polygon(p[0], p[1], geom)]

    return valid

def distribute_points_in_geom(geom,

def distribute_points_new(geom, nspts, ewpts):

    """

    This is a decision tree that attempts to perform a

    very simplistic approximation of the shape

    of the geometry and then place some number of

    north/south and east/west points into the geometry.

    Parameters

    ----------

    geom : shapely.geom

           A shapely geometry object

    nspts : int

            The number of points to attempt to place

            in the N/S (up/down) direction

    ewpts : int

            The number of points to attempt to place

            in the E/W (right/left) direction

    Returns

    -------

    valid : list

            of point coordinates in the form [(x1,y1), (x2,y2), ..., (xn, yn)]

    """

    geom_coords = np.column_stack(geom.exterior.xy)

    coords = np.array(list(zip(*geom.envelope.exterior.xy))[:-1])

    ll = coords[0]

    lr = coords[1]

    ur = coords[2]

    ul = coords[3]

    # Find the points nearest the ur and ll aligned // with eastern side of ground_poly

    elon, elat = find_side('east')

    ur_actual = np.array(two_point_extrapolate(ur[1], elat, elon))[::-1]

    lr_actual = np.array(two_point_extrapolate(lr[1],elat, elon))[::-1]

    wlon, wlat = find_side('west')

    ul_actual = np.array(two_point_extrapolate(ul[1], wlat, wlon))[::-1]

    ll_actual = np.array(two_point_extrapolate(ll[1], wlat, wlon))[::-1]

    dt = (ur_actual-ul_actual)*0.025 #some offset to make sure endpoints are within geom

    db = (lr_actual-ll_actual)*0.025

    newtop = create_points_along_line(ul_actual+dt, ur_actual-dt, ewpts)

    newbot = create_points_along_line(ll_actual+db, lr_actual-db, ewpts)

    points = []

    for i in range(len(newtop)):

        top = newtop[i]

        bot = newbot[i]

        line_of_points = create_points_along_line(top, bot, nspts)

        points.append(line_of_points)

    if len(points) < 1:

        return []

    points = np.vstack(points)

    # Perform a spatial intersection check to eject points that are not valid

    valid = [p for p in points if xy_in_polygon(p[0], p[1], geom)]

    return valid

def distribute_points_in_geom(geom, method="classic",

                              nspts_func=lambda x: ceil(round(x,1)*10),

                              ewpts_func=lambda x: ceil(round(x,1)*5)):

    """

            of valid points in the form (x,y) or (lon,lat)

    """

    point_funcs = {

        "classic" :  distribute_points_classic,

        "new" : distribute_points_new

    }

    point_distribution_func = point_funcs[method]

    coords = list(zip(*geom.envelope.exterior.xy))

    short = np.inf

    long = -np.inf

        if nspts == 1 and ewpts == 1:

            valid = single_centroid(geom)

        else:

            valid = distribute_points(geom, nspts, ewpts)

            valid = point_distribution_func(geom, nspts, ewpts)

    elif ew == True:

        # Since this is an LS, we should place these diagonally from the 'lower left' to the 'upper right'

        nspts = ewpts_func(short)

        if nspts == 1 and ewpts == 1:

            valid = single_centroid(geom)

        else:

            valid = distribute_points(geom, nspts, ewpts)

            valid = point_distribution_func(geom, nspts, ewpts)

    else:

        print('WTF Willy')

    return valid

from shapely.geometry.multipolygon import MultiPolygon

from shapely.geometry import Point

from redis import StrictRedis

from plurmy import Slurm

from autocnet import config, engine, Session

from autocnet.io.db.model import Images, Points, Measures

from autocnet.io.db.model import Images, Points, Measures, JsonEncoder

from autocnet.graph.network import NetworkCandidateGraph

from autocnet.matcher.subpixel import iterative_phase

from autocnet.cg.cg import distribute_points_in_geom

ctypes.CDLL(find_library('usgscsm'))

def generate_ground_points(ground_database, nspts_func=lambda x: int(round(x,1)*1), ewpts_func=lambda x: int(round(x,1)*4)):

    Ground_Session, ground_engine = new_connection(ground_database)

def generate_ground_points(ground_db_config, nspts_func=lambda x: int(round(x,1)*1), ewpts_func=lambda x: int(round(x,1)*4)):

    """

    Provided a config file which points to a database containing ground image path and geom information,

    generates ground points on these images within the range of a source database's images. For example,

    if creating a CTX mosaic which is grounded using themis data, the config files would look like:

        CTX_config -> located in config/[config_name].yml

        database:

            type: 'postgresql'

            username: 'jay'

            password: 'abcde'

            host: '130.118.160.193'

            port: 8085

            pgbouncer_port: 8083

            name: 'somename'

            timeout: 500

        Themis_config -> passed in to this function as ground_db_config

        ground_db_config = {'username':'jay',

                            'password':'abcde',

                            'host':'autocnet.wr.usgs.gov',

                            'pgbouncer_port':5432,

                            'name':'mars'}

    THIS FUNCTION IS CURRENTLY HARD CODED FOR themisdayir TABLE QUERY

    Parameters

    ----------

    ground_db_config : dict

                       In the form: {'username':'somename',

                                     'password':'somepassword',

                                     'host':'somehost',

                                     'pgbouncer_port':6543,

                                     'name':'somename'}

    nspts_func       : func

                       describes distribution of points along the north-south

                       edge of an overlap.

    ewpts_func       : func

                       describes distribution of points along the east-west

                       edge of an overlap.

    """

    warnings.warn('This function is not well tested. No tests currently exists \

    in the test suite for this version of the function.')

    Ground_Session, ground_engine = new_connection(ground_db_config)

    ground_session = Ground_Session()

    session = Session()

    ground_poly = wkt.loads(session.query(functions.ST_AsText(functions.ST_Union(Images.geom))).one()[0])

    fp_poly = wkt.loads(session.query(functions.ST_AsText(functions.ST_Union(Images.geom))).one()[0])

    session.close()

    image_fp_bounds = list(ground_poly.bounds)

    fp_poly_bounds = list(fp_poly.bounds)

    # just hard code queries to the mars database as it exists for now

    ground_image_query = f'select * from themisdayir where geom && ST_MakeEnvelope({image_fp_bounds[0]}, {image_fp_bounds[1]}, {image_fp_bounds[2]}, {image_fp_bounds[3]}, {config["spatial"]["latitudinal_srid"]})'

    ground_image_query = f'select * from themisdayir where ST_INTERSECTS(geom, ST_MakeEnvelope({fp_poly_bounds[0]}, {fp_poly_bounds[1]}, {fp_poly_bounds[2]}, {fp_poly_bounds[3]}, {config["spatial"]["latitudinal_srid"]}))'

    themis_images = gpd.GeoDataFrame.from_postgis(ground_image_query,

                                                  ground_engine, geom_col="geom")

    coords = distribute_points_in_geom(ground_poly, nspts_func=nspts_func, ewpts_func=ewpts_func)

    coords = distribute_points_in_geom(fp_poly, nspts_func=nspts_func, ewpts_func=ewpts_func)

    coords = np.asarray(coords)

    sql = """

    SELECT * FROM themisdayir as i WHERE ST_Contains(i.geom, ST_setsrid(ST_Point({}, {}), 949900))

    """

    records = []

    coord_list = []

    coord_id = []

    # throw out points not intersecting the ground reference images

    for i, coord in enumerate(coords):

        formated_sql = sql.format(coord[0], coord[1])

        res = ground_session.execute(formated_sql)

        for record in res:

        # res = ground_session.execute(formated_sql)

        p = Point(*coord)

        res = themis_images[themis_images.intersects(p)]

        for k, record in res.iterrows():

            record["pointid"] = i

            records.append(record)

            coord_list.append(Point(*coord))

            coord_list.append(p)

    ground_session.close()

    # start building the cnet

    ground_cnet = pd.DataFrame(data = records, columns = ['pointid', 'name', 'path', 'footprint', 'serial'])

    ground_cnet = pd.DataFrame.from_records(records)

    ground_cnet["point"] = coord_list

    ground_cnet['line'] = None

    ground_cnet['sample'] = None

    ground_cnet['resolution'] = None

    # generate lines and samples from ground points

    groups = ground_cnet.groupby('path')

    # group by images so campt can do multiple at a time

        lines = []

        samples = []

        resolutions = []

        indices = []

        for i, res in enumerate(point_list):

            if res[1].get('Error') is not None:

                print('Bad intersection')

                print('Bad intersection: ', res[1].get("Error"))

                lines.append(None)

                samples.append(None)

                resolutions.append(None)

        ground_cnet.loc[index, 'resolution'] = resolutions

    ground_cnet = gpd.GeoDataFrame(ground_cnet, geometry='point')

    return ground_cnet

    return ground_cnet, fp_poly, coords

def propagate_control_network(base_cnet):

def propagate_point(lon, lat, pointid, paths, lines, samples, resolutions):

    """

    """

    dest_images = gpd.GeoDataFrame.from_postgis("select * from images", engine, geom_col="geom")

    spatial_index = dest_images.sindex

    groups = base_cnet.groupby('pointid').groups

    # append to list if images, mostly used for working with the network in python

    # after this step, is this uncecceary outside of debugging? Maybe actually should return

    # more info of where everything was sourced in the original DataFrames?

    images = []

    images = gpd.GeoDataFrame.from_postgis(f"select * from images where ST_Intersects(geom, ST_SetSRID(ST_Point({lon}, {lat}), {config['spatial']['latitudinal_srid']}))", engine, geom_col="geom")

    # append CNET info into structured Python list

    constrained_net = []

    dbpoints = []

    dbmeasures = []

    image_measures = pd.DataFrame(zip(paths, lines, samples, resolutions), columns=["path", "line", "sample", "resolution"])

    measure = image_measures.iloc[0]

    # easily parrallelized on the cpoint level, dummy serial for now

    for cpoint, indices in groups.items():

        measures = base_cnet.loc[indices]

        measure = measures.iloc[0]

        p = measure.point

        # get image in he destination that overlap

        matches = dest_images[dest_images.intersects(p)]

    p = Point(lon, lat)

    new_measures = []

    # lazily iterate for now

        for i,row in matches.iterrows():

            res = isis.point_info(row["path"], p.x, p.y, point_type="ground", allow_outside=False)

    for i,image in images.iterrows():

        res = isis.point_info(image["path"], p.x, p.y, point_type="ground", allow_outside=False)

        dest_line, dest_sample = res["GroundPoint"]["Line"], res["GroundPoint"]["Sample"]

        try:

            dest_resolution = res["GroundPoint"]["LineResolution"].value

        except:

                warnings.warn(f'Failed to generate ground point info on image {row["path"]} at lat={p.y} lon={p.x}')

            warnings.warn(f'Failed to generate ground point info on image {image["path"]} at lat={p.y} lon={p.x}')

            continue

            dest_data = GeoDataset(row["path"])

        dest_data = GeoDataset(image["path"])

        dest_arr = dest_data.read_array()

        # list of matching results in the format:

        # [measure_index, x_offset, y_offset, offset_magnitude]

        match_results = []

        for k,m in image_measures.iterrows():

            # dynamically set scale based on point resolution

            dest_to_base_scale = dest_resolution/measure["resolution"]

            dest_to_base_scale = dest_resolution/m["resolution"]

            scaled_dest_line = (dest_arr.shape[0]-dest_line)*dest_to_base_scale

            scaled_dest_sample = dest_sample*dest_to_base_scale

            dest_arr = imresize(dest_arr, dest_to_base_scale)[::-1]

            scaled_dest_arr = imresize(dest_arr, dest_to_base_scale)[::-1]

            # list of matching results in the format:

            # [measure_index, x_offset, y_offset, offset_magnitude]

            match_results = []

            for k,m in measures.iterrows():

            base_arr = GeoDataset(m["path"]).read_array()

            sx, sy = m["sample"], m["line"]

            dx, dy = scaled_dest_sample, scaled_dest_line

            try:

                # not sure what the best parameters are here

                    ret = iterative_phase(sx, sy, dx, dy, base_arr, dest_arr, size=10, reduction=1, max_dist=1, convergence_threshold=1)

                ret = iterative_phase(sx, sy, dx, dy, base_arr, scaled_dest_arr, size=10, reduction=1, max_dist=2, convergence_threshold=1)

            except Exception as ex:

                match_results.append(ex)

                continue

        match_results = np.asarray([res for res in match_results if isinstance(res, list)])

        if match_results.shape[0] == 0:

            # no matches

            print("No Mathces")

            continue

        match_results = match_results[np.argwhere(match_results[:,3] == match_results[:,3].min())][0][0]

        if match_results[3] > 2:

                # best match diverged too much

            # best match drifted too much

            continue

            measure = measures.loc[match_results[0]]

        # apply offsets

        sample = (match_results[1]/dest_to_base_scale) + dest_sample

        line = (match_results[2]/dest_to_base_scale) + dest_line

            pointpvl = isis.point_info(row["path"], sample, line, point_type="image")

        pointpvl = isis.point_info(image["path"], sample, line, point_type="image")

        groundx, groundy, groundz = pointpvl["GroundPoint"]["BodyFixedCoordinate"].value

        groundx, groundy, groundz = groundx*1000, groundy*1000, groundz*1000

            images.append(row["path"])

            constrained_net.append({

                    'pointid' : cpoint,

                    'imageid' : row['id'],

                    'serial' : row.serial,

        new_measures.append({

                'pointid' : pointid,

                'imageid' : image['id'],

                'serial' : image['serial'],

                'line' : line,

                'sample' : sample,

                'point_latlon' : p,

                'point_ground' : Point(groundx, groundy, groundz)

        })

    return new_measures

def cluster_propagate_control_network(base_cnet, walltime='00:20:00', chunksize=1000, exclude=None):

    warnings.warn('This function is not well tested. No tests currently exists \

    in the test suite for this version of the function.')

    # Setup the redis queue

    rqueue = StrictRedis(host=config['redis']['host'],

                         port=config['redis']['port'],

                         db=0)

    # Push the job messages onto the queue

    queuename = config['redis']['processing_queue']

    groups = base_cnet.groupby('pointid').groups

    for cpoint, indices in groups.items():

        measures = base_cnet.loc[indices]

        measure = measures.iloc[0]

        p = measure.point

        # get image in the destination that overlap

        lon, lat = measures["point"].iloc[0].xy

        msg = {'lon' : lon[0],

               'lat' : lat[0],

               'pointid' : cpoint,

               'paths' : measures['path'].tolist(),

               'lines' : measures['line'].tolist(),

               'samples' : measures['sample'].tolist(),

               'resolutions' : measures['resolution'].tolist(),

               'walltime' : walltime}

        rqueue.rpush(queuename, json.dumps(msg, cls=JsonEncoder))

    # Submit the jobs

    submitter = Slurm('acn_propagate',

                 job_name='cross_instrument_matcher',

                 mem_per_cpu=config['cluster']['processing_memory'],

                 time=walltime,

                 partition=config['cluster']['queue'],

                 output=config['cluster']['cluster_log_dir']+'/autocnet.cim-%j')

    job_counter = len(groups.items())