# We don't technically need this but it avoids a warning when importing pysis
import os
os.environ['ISISROOT'] = '/usgs/cpkgs/anaconda3_linux/envs/isis3.9.0'

# AutoCNet Intro As mentioned earlier AutoCNet is a method for storing control networks and has outlier detection functionality. AutoCNet also contains a suite of functions that parallelize network generation that leverages and compliments ISIS processing. The advantage of AutoCNet network generation is it takes advantage of elementwise cluster processing (these elements can be images, points, measures, etc.) and postgresql for data storage and quick relational querying.

In this notebook we are going to step through the network generation process in AutoCNet!

For Quick Access: - Load and apply configuration file - Ingest images and calculate overlaps - Distribute points in overlaps - Subpixel register points

Grab the Image Data

We are going to process Kaguya Terrain Camera (TC) images surrounding the Reiner Gamma Lunar Swirl (4.9° - 9.9° N Planetocentric Latitude and 61.3° - 56.3° W Longitude). The data is located in ‘/scratch/ladoramkershner/moon/kaguya/workshop/original/’, please use the cell below to copy the data into a directory of your choosing.

import getpass
uid = getpass.getuser()

output_directory = f'/scratch/ladoramkershner/FY21_autocnet_workshop/workshop_scratch/{uid}' # put output directory path as string here
print(output_directory)

# copy over the data to the 'lvl1' subdirectory
!mkdir -p $output_directory/lvl1/
!cp -p /scratch/ladoramkershner/moon/kaguya/workshop/original/*cub $output_directory/lvl1/

We need to create a list of the cubes, to feed into AutoCNet. It is important that the cube list handed to AutoCNet contain absolute paths, as they will serve as an accessor for loading information from the cubes later.

!ls $output_directory/lvl1/*cub > $output_directory/cubes.lis
!head $output_directory/cubes.lis

# Parse the Configuration File Return To Top

The configuration parameters are typically held in a configuration yaml file. A configuration file has been compiled for use internal to the USGS ASC facilities leveraging a shared cluster and database. Use AutoCNet’s function ‘parse_config’ to read in the yaml file and output a dictionary variable.

from autocnet.config_parser import parse_config

config_path = '/scratch/ladoramkershner/FY21_autocnet_workshop/config_moon.yml'
config = parse_config(config_path)

The config is a nested dictionary, meaning it has a larger dictionary structure defining sections for the services above and then each service section is a dictionary defining the particular configuration parameters.

import numpy as np

print('configuration dictionary keys: ')
print(np.vstack(list(config.keys())), '\n')

print('cluster configuration dictionary keys: ')
print(np.vstack(list(config['cluster'].keys())))

Although the configuration file is set up for internal use, some fields need to be altered to point to user specific areas or unique strings.

config['cluster']['cluster_log_dir'] = f'/scratch/ladoramkershner/FY21_autocnet_workshop/workshop_scratch/{uid}/logs'
config['database']['name'] = f'workshop_{uid}_kaguyatc_reinergamma'
config['redis']['basename'] = f'{uid}_queue'
config['redis']['completed_queue'] = f'{uid}_queue:comp'
config['redis']['processing_queue'] = f'{uid}_queue:proc'
config['redis']['working_queue'] = f'{uid}_queue:work'

default_log = config['cluster']['cluster_log_dir']

print(f'your log directory: {default_log}')
print('your database name:', config['database']['name'])

Create the NetworkCandidateGraph

The NetworkCandidateGraph (NCG) class can be instantiated to an object without any arguments. However, this NCG object requires configuration before it can be used for any meaningful work, so we have to run ‘config_from_dict’.

from autocnet.graph.network import NetworkCandidateGraph

ncg = NetworkCandidateGraph()
ncg.config_from_dict(config)
ncg.from_database()

# Ingest Image Data and Calculate Overlaps Return To Top

At this point our ncg variable is empty, so if we try to plot the contents we will get an empty plot.

ncg.plot()

We need to load the images into the ncg using ‘add_from_filelist’, which loads the images from a passed in list and then calculates the overlaps.

filelist = f'{output_directory}/cubes.lis' # this should contain absolute paths
ncg.add_from_filelist(filelist)

Now when we plot the ncg, we see the undirected graph, where the circles are the nodes/images and the lines are the edges/overlaps. The Kaguya TC data has a very regular overlap pattern in this area, seen by the large number of edges shared between nodes.

ncg.plot()

We have access to the image data through the ncg, but the ncg does not persist after the notebook is shut down. To persist the network, AutoCNet leverages a database for the storage of the networks images, points, and measures. The ncg has access to this database through the ncg’s ‘session_scope’. Through the ncg.session_scope() you can interact and execute queries on your database in pure SQL.

with ncg.session_scope() as session:
    img_count = session.execute("SELECT COUNT(*) FROM images").fetchall()
    overlap_count = session.execute("SELECT COUNT(*) FROM overlay").fetchall()
    print('  Number of images in database: ', img_count[0][0])
    print('Number of overlaps in database: ', overlap_count)

session.execute() is equivalent to running the input string in the database directly. It is a convenient if you are already familiar with pure sql commands, however, the return values are messy. The session.query() leverages a python module called sqlalchemy which allow pythonic calls to your database with clean output.

from autocnet.io.db.model import Images, Overlay

with ncg.session_scope() as session:
    img_count = session.query(Images).count()
    overlap_count = session.query(Overlay).count()
    print('  Number of images in database: ', img_count)
    print('Number of overlaps in database: ', overlap_count)

Additionally, session.execute() can be inconvenient if working with the actual data contained within the tables. For example, to access certain information you need to know the index where that information exists.

with ncg.session_scope() as session:
    img = session.execute("SELECT * FROM images LIMIT 1").fetchall()
    print('image index: ', img[0][0])
    print('product id: ', img[0][1])
    print('image path: ', img[0][2])
    print('image serial number: ', img[0][3])
    print('image ignore flag: ', img[0][4])
#     print('image geom: ', img[0][5]) # only uncomment after looking at other output
    print('image camera type: ', img[0][7])

However, if the structure of the database changes (order of the columns or a column is added/removed) or you cannot remember the order of the columns, working with the database data in this way is be very inconvenient. So AutoCNet built models for each table of the database tables to help interface with them.

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

with ncg.session_scope() as session:
    img = session.query(Images).first()
    print('image index: ', img.id)
    print('product id: ', img.name)
    print('image path: ', img.path)
    print('image serial number: ', img.serial)
    print('image ignore flag: ', img.ignore)
#     print('image geometry: ', img.geom) # only uncomment after looking at other output
    print('image camera type: ', img.cam_type)

Accessing the information off of the img object is more intuitive as it is property based instead of index based.

img[0][0] –> img.id img[0][1] –> img.name img[0][2] –> img.path and so on..

Additionally, if you cannot remember the exact names of the properties you want to access, you can dir() the model.

print(dir(Images))

Finally, if you uncomment the prints of the geometry in the two previous cells you see that the raw database geometry (given by session.execute()) is stored as a hexadecimal string while the Images.geom property is a shapely Multipolygon with more intuitive longitude and latitude values. The MultiPolygon also has helpful functions which allows direct access to the latitude, longitude information. To plot the geometry all we have to do is…

import matplotlib.pyplot as plt

n = 25
with ncg.session_scope() as session:
    imgs = session.query(Images).limit(n)

    fig, axs = plt.subplots(1, 1, figsize=(5,10))
    axs.set_title(f'Footprints of First {n} Images in Database')
    for img in imgs:
        x,y = img.geom.envelope.boundary.xy # this call!
        axs.plot(x,y)

# Place Points in Overlap Return To Top

The next step in the network generation process is to lay down points in the image overlaps. Before dispatching the function to the cluster, we need to make the log directory from our configuration file. If a SLURM job is submitted with a log directory argument that does not exist, the job will fail.

ppio_log_dir = default_log.replace('logs', 'ppio_logs')

print('creating directory: ', ppio_log_dir)
if not os.path.exists(ppio_log_dir):
    os.mkdir(ppio_log_dir)

We are going to use the ‘place_points_in_overlap’ function to lay the points down. For now, we will use the default size and distribution arguments, which is easily accomplished by not handing in values for these arguments. However, we need to change our camera type from the default ‘csm’ to ‘isis’.

from autocnet.spatial.overlap import place_points_in_overlap

njobs = ncg.apply('spatial.overlap.place_points_in_overlap',
                  on='overlaps', # start of function kwargs
                  cam_type='isis',
                  walltime='00:30:00', # start of apply kwargs
                  log_dir=ppio_log_dir,
                  arraychunk=50)

print(njobs)

!squeue -u $uid | head # helpful to grab job array id

The ‘place_points_in_overlaps’ function first evenly distributes points spatially into a given overlap, then it back-projects the points into the ‘top’ image. Once in image space, the function searches the area surrounding the measures to find interesting features to shift the measures to (this increases the chance of subpixel registration passing). The shifted measures are projected back to the ground and these updated longitudes and latitudes are used to propagate the points into all images associated with the overlap. So, this function requires: - An overlap (to evenly distribute points into) - Distribution kwargs (to decide how points are distributed into the overlap) - Size of the area around the measure (to search for the interesting feature) - Camera type (so it knows what to expect as inputs/output for the camera model)

Since this function operates independently on each overlap, it is ideal for parallelization with the cluster. Notice that we are not passing in a single overlap to the apply call, instead we pass “on = ‘overlaps’”. The ‘on’ argument indicates which element (image, overlap, point, measure) to apply the function.

with ncg.session_scope() as session:
    noverlay = session.query(Overlay).count()
    print(noverlay)

Multiple Ways to Check Job Array Process

Log Files

As jobs are put on the cluster, their corresponding log files are created. You can check how many jobs have been/ are being processed on the cluster by looking in the log directory.

!ls $ppio_log_dir | head -5

As more logs are placed in the log directory, you will have to specify which array job’s logs you are checking on. The naming convention of the log files generated by AutoCNet are ‘path.to.function.function_name-jobid.arrayid_taskid.out’

jobid = '' # put jobid int here
!ls $ppio_log_dir/*${jobid}_*.out | wc -l

Slurm Account

Using ‘sacct’ allows you to check the exit status of the tasks from your job array.

!sacct -j $jobid -s 'completed' | wc -l
!sacct -j $jobid -s 'failed'    | wc -l
!sacct -j $jobid -s 'timeout'   | wc -l
!sacct -j $jobid -s 'cancelled' | wc -l

The return of ‘2’ from the word count on the ‘failed’, ‘timeout’, and ‘cancelled’ job accounts are the header lines.

!sacct -j $jobid -s 'failed' | head

NCG Queue Length

The queue holds the job packages in json files called ‘queue messages’ until the cluster is ready for the job. You can view how many messages are left on the queue with the ‘queue_length’ NCG property.

print("jobs left on the queue: ", ncg.queue_length)

Reapply to the Cluster?

Sometimes jobs fail to submit to the cluster, it is prudent to check the ncg.queue_length AFTER your squeue is empty.

!squeue -u $uid

print("jobs left on the queue: ", ncg.queue_length)

When reapplying a function to the cluster, you do not need to resubmit the function arguments, because those were already serialized into the queue message. However, the cluster submission arguments can be reformatted and the ‘reapply’ argument should be set to ‘True’.

# njobs = ncg.apply('spatial.overlap.place_points_in_overlap',
#                         chunksize=redis_orphans,
#                         arraychunk=None,
#                         walltime='00:20:00',
#                         log_dir=ppio_log_dir,
#                         reapply=True)
# print(njobs)

One advantage of using of a postgresql database for data storage is that it allows for storage of geometries. You can then use relational queries to view how different elements’ geometries relate with one another.

from geoalchemy2 import functions
from geoalchemy2.shape import to_shape

with ncg.session_scope() as session:
    results = (
        session.query(
        Overlay.id,
        Overlay.geom.label('ogeom'),
        Points.geom.label('pgeom')
        )
        .join(Points, functions.ST_Contains(Overlay.geom, Points.geom)=='True')
        .filter(Overlay.id < 10) # Just view first 10 overlaps
        .all()
    )
    print('number of points: ', len(results))

    fig, axs = plt.subplots(1, 1, figsize=(10,10))
    axs.grid()

    oid = []
    for res in results:
        if res.id not in oid:
            oid.append(res.id)
            ogeom = to_shape(res.ogeom)
            ox, oy = ogeom.envelope.boundary.xy
            axs.plot(ox, oy, c='k')
        pgeom = to_shape(res.pgeom)
        px, py = pgeom.xy
        axs.scatter(px, py, c='grey')

Notice that the points are not in straight lines, this is because of the shifting place_points_in_overlap does to find interesting measure locations.

However, the default distribution of points in the overlaps looks sparse, so let’s rerun place_points_in_overlap with new distribution kwargs. Before rerunning place_points_in_overlap, the points and measures tables need to be cleared using ncg’s ‘clear_db’ method.

with ncg.session_scope() as session:
    npoints = session.query(Points).count()
    print('number of points: ', npoints)

    nmeas = session.query(Measures).count()
    print('number of measures: ', nmeas)

ncg.clear_db(tables=['points', 'measures']) # clear the 'points' and 'measures' database tables

with ncg.session_scope() as session:
    npoints = session.query(Points).count()
    print('number of points: ', npoints)

    nmeas = session.query(Measures).count()
    print('number of measures: ', nmeas)

The distribution argument for place_points_in_overlap requires two function inputs. Since overlaps are variable shapes and sizes, one integer is not sufficient to determine effective gridding along every overlaps sides. Instead, the distribution of points along the N to S edge of the overlap and the E to W edge of the overlap are determined based on a function.

The default distribution functions are: nspts_func=lambda x: ceil(round(x,1)*10) ewpts_func=lambda x: ceil(round(x,1)*5)

Where x in nspts_func is the length of the overlap’s longer edge (in km) and x in ewpts_func is the length of the overlap’s shorter edge (in km). This way a shorter edge will receive less points and a longer side will receive more points. Change the multipliers in the ‘ns’ and ‘ew’ functions below to find a satisfying distribution.

from autocnet.cg.cg import distribute_points_in_geom

def ns(x):
    from math import ceil
    return ceil(round(x,1)*15)

def ew(x):
    from math import ceil
    return ceil(round(x,1)*10)

total=0
with ncg.session_scope() as session:
    srid = config['spatial']['latitudinal_srid']
    overlaps = session.query(Overlay).filter(Overlay.geom.intersects(functions.ST_GeomFromText('LINESTRING(301.2 7.4, 303.7 7.4, 303.7 9.9, 301.2 9.9, 301.2 7.4)', srid))).all()

    print('overlaps in selected area: ', len(overlaps))
    for overlap in overlaps:
        ox, oy = overlap.geom.exterior.xy
        plt.plot(ox,oy)

        valid = distribute_points_in_geom(overlap.geom, method='classic', nspts_func=ns, ewpts_func=ew, Session=session)
        if valid:
            total += len(valid)
            px, py = list(zip(*valid))
            plt.scatter(px, py, s=1)

    print('  points in selected area: ', total)

Then rerun the apply function, setting the ‘distribute_points_kwargs’ arguments.

distribute_points_kwargs = {'nspts_func':ns, 'ewpts_func':ew, 'method':'classic'}

njobs = ncg.apply('spatial.overlap.place_points_in_overlap',
                  on='overlaps', # start of function kwargs
                  distribute_points_kwargs=distribute_points_kwargs, # NEW LINE
                  cam_type='isis',
                  walltime='00:30:00', # start of apply kwargs
                  log_dir=ppio_log_dir,
                  arraychunk=100)
print(njobs)

Check the progress of your jobs.

!squeue -u $uid | wc -l
!squeue -u $uid | head

Count the number of jobs started by looking for generated logs.

jobid = '' # put jobid int here
! ls  $ppio_log_dir/*$jobid* | wc -l

Monitor how many jobs have completed or failed.

!sacct -j $jobid -s 'completed' | wc -l
!sacct -j $jobid -s 'failed'    | wc -l

Check to see if the ncg redis queue is clear.

redis_orphans = ncg.queue_length
print("jobs left on the queue: ", redis_orphans)

Reapply cluster job if there are still jobs left on the queue

# njobs = ncg.apply('spatial.overlap.place_points_in_overlap',
#                         chunksize=redis_orphans,
#                         arraychunk=None,
#                         walltime='00:20:00',
#                         log_dir=log_dir,
#                         reapply=True)
# print(njobs)

Visualize the new distribution

with ncg.session_scope() as session:
    results = (
        session.query(
        Overlay.id,
        Overlay.geom.label('ogeom'),
        Points.geom.label('pgeom')
        )
        .join(Points, functions.ST_Contains(Overlay.geom, Points.geom)=='True')
        .filter(Overlay.id < 10)
        .all()
    )
    print('number of points: ', len(results))

    fig, axs = plt.subplots(1, 1, figsize=(10,10))
    axs.grid()

    oid = []
    for res in results:
        if res.id not in oid:
            oid.append(res.id)
            ogeom = to_shape(res.ogeom)
            ox, oy = ogeom.envelope.boundary.xy
            axs.plot(ox, oy, c='k')
        pgeom = to_shape(res.pgeom)
        px, py = pgeom.xy
        axs.scatter(px, py, c='grey')

# Subpixel Registration Return To Top

The next step is to subpixel register the measures on the newly laid points, to do this we are going to use the ‘subpixel_register_point’ function. As the name suggests, ‘subpixel_register_point’ registers the measures on a single point, which makes it parallelizable on a network’s points. Before we fire off the cluster jobs, let’s create a new subpixel registration log directory.

subpix_log_dir = default_log.replace('logs', 'subpix_logs')

print('creating directory: ', subpix_log_dir)
if not os.path.exists(subpix_log_dir):
    os.mkdir(subpix_log_dir)

from autocnet.matcher.subpixel import subpixel_register_point

?subpixel_register_point
# ncg.apply?

subpixel_template_kwargs = {'image_size':(81,81), 'template_size':(51,51)}

njobs = ncg.apply('matcher.subpixel.subpixel_register_point',
                  on='points', # start of function kwargs
                  match_kwargs=subpixel_template_kwargs,
                  geom_func='simple',
                  match_func='classic',
                  cost_func=lambda x,y:y,
                  threshold=0.6,
                  verbose=False,
                  walltime="00:30:00", # start of apply kwargs
                  log_dir=subpix_log_dir,
                  arraychunk=200,
                  chunksize=20000) # maximum chunksize = 20,000

print(njobs)

Check the progress of your jobs.

!squeue -u $uid | head

This function chooses a reference measure, affinely transforms the other images to the reference image, and clips an ‘image’ chip out of the reference image and a ‘template’ chip out of the transformed images. The template chips are marched across the image chip and the maximum correlation value and location is saved.

The solution is then evaluated to see if the maximum correlation solution is acceptable. The evaluation is done using the ‘cost_func’ and ‘threshold’ arguments. The cost_func is dependent two independent variables, the first is the distance that a point has shifted from the starting location and the second is the correlation coefficient coming out of the template matcher. The order that these variables are passed in matters. We are not going to consider the distance the measures were moved and just look at the maximum correlation value returned by the matcher. So, our function is simply: \(y\).

If the cost_func solution is greater than the threshold value, the registration is successful and the point is updated. If not, the registration is considered unsuccessful, the point is not updated, and is set to ignore.

So, ‘subpixel_register_point’ requires the following arguments: - pointid - match_kwargs (image size, template size) - cost_func - threshold

Count the number of jobs started by looking for generated logs.

jobid = '' # put jobid int here
! ls $subpix_log_dir/*$jobid* | wc -l

Monitor how many jobs have completed or failed.

!sacct -j $jobid -s 'completed' | wc -l
!sacct -j $jobid -s 'failed'    | wc -l

Check to see if the ncg redis queue is clear once squeue is empty.

redis_orphans = ncg.queue_length
print("jobs left on the queue: ", redis_orphans)

Reapply cluster job if there are still jobs left on the queue.

# job_array = ncg.apply('matcher.subpixel.subpixel_register_point',
#                       reapply=True,
#                       chunksize=redis_orphans,
#                       arraychunk=None,
#                       walltime="00:30:00",
#                       log_dir=subpix1_log_dir)
# print(job_array)

Visualize Point Registration

from plio.io.io_gdal import GeoDataset
from autocnet.transformation import roi
from autocnet.utils.utils import bytescale

roi_size = 25

with ncg.session_scope() as session:
    measures = session.query(Measures).filter(Measures.template_metric < 0.8, Measures.template_metric!=1).limit(15)
    for meas in measures:
        pid = meas.pointid
        source = session.query(Measures, Images).join(Images, Measures.imageid==Images.id).filter(Measures.pointid==pid, Measures.template_metric==1).all()
        sx = source[0][0].sample
        sy = source[0][0].line
        s_roi = roi.Roi(GeoDataset(source[0][1].path), sx, sy, size_x=roi_size, size_y=roi_size)
        s_image = bytescale(s_roi.clip())

        destination = session.query(Measures, Images).join(Images, Measures.imageid==Images.id).filter(Measures.pointid==pid, Measures.template_metric!=1).limit(1).all()
        dx = destination[0][0].sample
        dy = destination[0][0].line
        d_roi = roi.Roi(GeoDataset(destination[0][1].path), dx, dy, size_x=roi_size, size_y=roi_size)
        d_template = bytescale(d_roi.clip())

        fig, axs = plt.subplots(1, 2, figsize=(10,10));
        axs[0].imshow(s_image, cmap='Greys');
        axs[0].scatter(image_size[0], image_size[1], c='r')
        axs[0].set_title('Reference');
        axs[1].imshow(d_template, cmap='Greys');
        axs[1].scatter(image_size[0], image_size[1], c='r')
        axs[1].set_title('Template');

We are going to rerun the subpixel registration with larger chips to attempt to register the measures that failed first run.

subpixel_template_kwargs = {'image_size':(221,221), 'template_size':(81,81)}

Additionally, ‘subpixel_register_point’ can be run on a subset of points, using either the ‘filters’ or the ‘query_string’ arguments.

The ‘filters’ argument does a equals comparison on point properties and filters out points with a certain property value (e.g.: points where ignore=true). While the ‘query_string’ argument can perform inequalities and applies on the selected values. Some examples of possible filters and query_string values are

filters = {'ignore': 'true'} # filters out points where ignore=true
query_string = """
        SELECT DISTINCT pointid FROM measures
        WHERE "templateMetric" < 0.65
        """ # only grabs points with template metrics less than 0.65

filters and query_string cannot be applied at the same time. So choose one, comment out the other argument’s line and rerun the subpixel registration apply.

njobs = ncg.apply('matcher.subpixel.subpixel_register_point',
                  on='points', # start of function kwargs
#                   filters=filters,  ##### NEW LINE
                  query_string=query_string,
                  match_kwargs=subpixel_template_kwargs,
                  geom_func='simple',
                  match_func='classic',
                  cost_func=lambda x,y:y,
                  threshold=0.6,
                  verbose=False,
                  walltime="00:30:00", # start of apply kwargs
                  log_dir=subpix_log_dir,
                  arraychunk=50,
                  chunksize=20000) # maximum chunksize = 20,000

print(njobs)

Subsequent runs of ‘subpixel_register_point’ can be run on all points, if this is done AutoCNet checks for a previous subpixel registration result in the database and only updates the geometry if the new result is better. Additionally, the apriori geometry (original camera pointing) is stored in the database and subpixel registration is always done on the apriori geometry to avoid ‘measure walking’. ‘Measure walking’ refers to a measure’s geometry moving further and further away from its original location due to multiple runs of subpixel registration. In traditional software the apriori geometry may be stored, but subpixel registration is typically run off of the current geometry.

Count the number of jobs started by looking for generated logs.

! squeue -u $uid | wc -l
! squeue -u $uid | head

Count the number of jobs started by looking for generated logs.

jobid = '' # put jobid int here

! ls $log_dir/*$jobid* | wc -l

Check to see if the ncg redis queue is clear once squeue is empty.

redis_orphans = ncg.queue_length
print("jobs left on the queue: ", redis_orphans)

Reapply cluster job if there are still jobs left on the queue.

# njobs = ncg.apply('matcher.subpixel.subpixel_register_point',
#                   reapply = True,
#                   walltime="00:30:00",
#                   log_dir='/scratch/ladoramkershner/mars_quads/oxia_palus/subpix2_logs/',
#                   arraychunk=50,
#                   chunksize=20000) # maximum chunksize = 20,000

# print(njobs)

subpix2: Write out Network

Once you are finished leverage AutoCNet tools and want to move onto ISIS based analysis (qnet, jigsaw, etc.), you can use the ncg.to_isis() function to write the information in your database to an ISIS control network file.

cnet = 'reiner_gamma_morning_ns7_ew5_t121x61_t221x81.net'
ncg.to_isis(os.path.join(output_directory,cnet))