Using PLIO to analyze control networks

PLIO is a general purpose library for reading data from various sources. In this workshop, we will be using PLIO’s ability to read ISIS control networks into a Pandas dataframe.

# PLIO uses pysis for some other things. We don't technically need this but it avoids a warning.
import os
os.environ['ISISROOT'] = '/usgs/cpkgs/anaconda3_linux/envs/isis4.3.0'
os.environ['ISISDATA'] = '/usgs/cpkgs/isis3/isis_data'

# 3D plotting toolkit for matplotlib
from mpl_toolkits.mplot3d import Axes3D

# Numerical Python library
import numpy as np

Our networks

All of this data was generously provided by Lynn Weller and Mike Bland from their Europa control project.

The first network is a very rough starting point. The Galileo images of Europa were put through the findfeatures application and then all of the resulting networks were merged together. This network has many known issues including islands, massive residuals, and poor coverage.

The second network is the final network containing Galileo and Voyager images of Europa. The issues from the initial network have been resolved and the final point cloud covers the majority of the body.

galileo_net = '/scratch/jmapel/europa/networks/GLL_FFCombined_thin_SubReg2_Del_2.net'
final_net = '/scratch/jmapel/europa/networks/GalileoVoyager_Europa_Merged_2020_CilixFree.net'

The control network dataframe

PLIO directly ingests the data from the control network file. Each row in the dataframe is a single control measure and each column is a field from the protobuf control network. The data for control points is stored implicitly in its measures.

# This function is what reads a control network file
from plio.io.io_controlnetwork import from_isis

galileo_df = from_isis(galileo_net)
galileo_df.describe()

Exercise: How many measures are there in the network? How many points are there in the network? How many images are there in the network?

tip: use len(dataframe) to find the number of rows in a dataframe

tip: use dataframe[“columnName”].nunique() to find the number of unique values in a column

Data types

The different columns of our dataframe store different types of data. The cell below shows all of the the data types in the dataframe. You can see all of the different possible datatypes for a dataframe in the pandas docs.

galileo_df.dtypes

Most of the data types are straightforward. For example, the line and sample are 64-bit floats. Let’s dig into the more unusual types.

pointType, measureType, aprioriSurfPointSource, and aprioriRadiusSource are 64 bit integers, but those integers correspond to enumerations. For example, a pointType of 2 means Free. See the tables below for all of the enumerations

galileo_df[['pointType', 'measureType', 'aprioriSurfPointSource']].head()

pointType

Value	Name
0	Tie (obsolete)
1	Ground (obsolete)
2	Free
3	Constrained
4	Fixed

measureType

Value	Name
0	Candidate
1	Manual
2	RegisteredPixel
3	RegisteredSubPixel

aprioriSurfPointSource & aprioriRadiusSource

Value	Name
0	None
1	User
2	AverageOfMeasures
3	Reference
4	Ellipsoid
5	DEM
6	Basemap
7	BundleSolution

Exercise: Have any measure in this network been sub-pixel registered?

tip: look at the measure types

id, pointChoosername, pointDatetime, aprioriSurfPointSourceFile, aprioriRadiusSourceFile, serialnumber, measureChoosername, and measureDatetime are all listed as objects but are simply strings.

galileo_df[['id', 'serialnumber', 'pointChoosername', 'pointDatetime', 'measureChoosername', 'measureDatetime']].head()

adjustedCovar, pointLog, and measureLog are more complicated. We will go over adjustedCovar later with the final Euroap network. pointLog is leftover from older network formats and can be ignored. measureLog contains information about the registration of the measure.

galileo_df.loc[1,'measureLog']

Data availability

Depending on how your network was generated and what processing has been done, many fields will not be set. If a numerical field has a value of 0, then it has not been set. For example, our network has not been bundle adjusted, so there are only a priori ground points.

galileo_df[['aprioriX', 'aprioriY', 'aprioriZ', 'adjustedX', 'adjustedY', 'adjustedZ']].describe()

Exercise: Can you find all of the fields that are completely unset in our control network?

tip: numerical fields default to 0, strings default to an empty string ““, and boolean values default to False.

You can also check which columns are default programmaticaly. The following cell checks if all of the values in a column are a default value.

(galileo_df==0).all() | (galileo_df=="").all() | (galileo_df==False).all()

Looking at a bundle adjusted control network

Our Galileo network is interesting but networks have significantly more useful information in them after bundle adjustment. So, let’s take a look at the final Europa network.

final_net_df = from_isis(final_net)
final_net_df.describe()

Exercise: What fields are set in the bundle adjusted network that weren’t previously?

Analyzing the measures

The data in a control network dataframe is not always in the format we want to work with. The measure residuals are broken down into the line and sample residuals. The following cell computes the full magnitude of the residuals and adds it to the dataframe under the “residualMag” column.

final_net_df['residualMag'] = np.sqrt(final_net_df['sampleResidual']**2 + final_net_df['lineResidual']**2)

Now let’s plot the residuals and see if we can form any theories. The next cell imports matplotlib for plotting tools and then plots the residuals in terms of sample and line residual. Note that the color of points is based on the residual magnitude, whcih should give a nice bullseye effect.

# This allows us to interact with our plots. This must be set before importing pyplot
%matplotlib notebook

# General plotting library
import matplotlib
import matplotlib.pyplot as plt

resid_fig = plt.figure(figsize=(6, 6))
resid_ax = resid_fig.add_subplot(111)
resid_scatter = resid_ax.scatter(final_net_df['sampleResidual'], final_net_df['lineResidual'], c=final_net_df['residualMag'], marker='+')
resid_ax.set_aspect('equal')
plt.axhline(0, color='black')
plt.axvline(0, color='black')
resid_cbar = plt.colorbar(resid_scatter)
resid_fig.suptitle('Bundle Adjusted Measure Residuals')
resid_ax.set_xlabel('Sample Residual')
resid_ax.set_ylabel('Line Residual')
resid_cbar.set_label('Residual Magnitude')
plt.show()

We can also color our points based on other properties. Let’s try and separate the measures out by mission. The serial numbers should help us so let’s look at the serial numbers for all of our images.

final_net_df['serialnumber'].unique()

Each serial number starts with the mission name, which makes separating them out easy. All we need to do is check if the beginning of the serial number matches our mission.

The pd.DataFrame.str package allows us to do this type of string comparisons quickly and easily. Here we will use the DataFrame.str.startswith method.

final_galileo_df = final_net_df[final_net_df['serialnumber'].str.startswith('Galileo')]
final_voyager1_df = final_net_df[final_net_df['serialnumber'].str.startswith('Voyager1')]
final_voyager2_df = final_net_df[final_net_df['serialnumber'].str.startswith('Voyager2')]

Now let’s plot the measures and color them based on their mission.

inst_resid_fig = plt.figure(figsize=(6, 6))
inst_resid_ax = inst_resid_fig.add_subplot(111)
inst_resid_ax.scatter(final_galileo_df['sampleResidual'], final_galileo_df['lineResidual'], color='Green', marker='+', alpha=0.25, label='Galileo')
inst_resid_ax.scatter(final_voyager1_df['sampleResidual'], final_voyager1_df['lineResidual'], color='Red', marker='+', alpha=0.25, label='Voyager1')
inst_resid_ax.scatter(final_voyager2_df['sampleResidual'], final_voyager2_df['lineResidual'], color='Blue', marker='+', alpha=0.25, label='Voyager2')
inst_resid_ax.set_aspect('equal')
plt.axhline(0, color='black')
plt.axvline(0, color='black')
plt.legend()
inst_resid_fig.suptitle('Bundle Adjusted Measure Residuals by Mission')
inst_resid_ax.set_xlabel('Sample Residual')
inst_resid_ax.set_ylabel('Line Residual')
plt.show()

What can you say about the residuals for the different missions based on our plot?

Exercise: What the descriptive statistics for the residual magnitude of the Galileo measures? What about for Voyager 1 and Voyager 2?

final_galileo_df['residualMag'].describe()

final_voyager1_df['residualMag'].describe()

final_voyager2_df['residualMag'].describe()

Do you notice anything interesting about the residual magnitudes for the different instruments? How does this compare to what you noticed with the scatter plot?

We can even test if the measure residuals are normally distributed. The following cell performs a chi-squared test to see if the residual magnitudes could reasonably come from a normal distribution. This is important because it will tell us if we have large blunders in our network or systematic error from something like a bad sensor model.

# Statistics library
from scipy import stats

alpha = 1e-3 # 99.999% confidence
_, normal_test_result = stats.normaltest(final_voyager1_df['residualMag'])
print(f'Chi-squared test statistic: {normal_test_result}')
if (normal_test_result < alpha):
    print("The residuals are normally distributed")
else:
    print("The residuals may not be normally distributed")

Analyzing the points

The information for control points is duplicated for each measure they have. So, the first step in looking at control point data is to extract only the data we want from the dataframe. This will make the dataframe easier to read and it will make things run quicker.

To do this, we’re going to first extract all of the columns with point data. Then, we’re going extract the first measure from each point. After all is said and done, we will have a dataframe with columns related to the point info and only one row for each point.

point_columns = ['id',
                 'pointType',
                 'pointChoosername',
                 'pointDatetime',
                 'pointEditLock',
                 'pointIgnore',
                 'pointJigsawRejected',
                 'aprioriSurfPointSource',
                 'aprioriSurfPointSourceFile',
                 'aprioriRadiusSource',
                 'aprioriRadiusSourceFile',
                 'latitudeConstrained',
                 'longitudeConstrained',
                 'radiusConstrained',
                 'aprioriX',
                 'aprioriY',
                 'aprioriZ',
                 'aprioriCovar',
                 'adjustedX',
                 'adjustedY',
                 'adjustedZ',
                 'adjustedCovar',
                 'pointLog']
final_points_df = final_net_df[point_columns].drop_duplicates('id')
final_points_df.describe()

Next, we’re going to transform the point data so that it’s more useful to us. This cell will take the (X, Y, Z) adjusted ground points and convert them to (lat, lon, radius) using a library called pyproj. pyproj is a very powerful projections library and can do many cartofraphic transformations and projections.

Note: This cell will generate a warning because we are using old pyproj.Proj calls which will eventually need to change. For now we can ignore the warning.

# Projection library for switching between rectangular and latitudinal
os.environ['PROJ_LIB'] = '/usgs/cpkgs/anaconda3_linux/envs/autocnet/share/proj'
import pyproj

# Compute the lat/lon/alt
europa_radii = [1562600, 1560300, 1559500]
ecef = pyproj.Proj(proj='geocent', a=europa_radii[0], b=europa_radii[1], c=europa_radii[2])
lla = pyproj.Proj(proj='latlong', a=europa_radii[0], b=europa_radii[1], c=europa_radii[2])
lon, lat, alt = pyproj.transform(ecef, lla, final_points_df['adjustedX'].values, final_points_df['adjustedY'].values, final_points_df['adjustedZ'].values, radians=True)

# Store the data in the dataframe
final_points_df['latitude'] = lat
final_points_df['longitude'] = lon
final_points_df['altitude'] = alt

# We will also want the point radii
final_points_df['radius'] = np.sqrt(final_points_df['adjustedX']**2 + final_points_df['adjustedY']**2 + final_points_df['adjustedZ']**2)

Because of how we defined our projection, the latitude and longitude values will be in radians. Also, the longitude will be in 180 postiive East. You can change this by modifying how you use pyproj but that is outside of this workshop.

final_points_df[["latitude", "longitude", "altitude", "radius"]].describe()

Exercise: Convert the latitude and longitude from radians to degrees:

Similar to how we computed the residual magnitude, we want to compute the average residual magnitude for each point. The following cell goes back to our original dataframe, computes the mean point by point, and then saves all of the results in our new dataframe.

Note: This cell can take a while to run because it has to re-access the dataframe for every point

final_points_df["averageResidual"] = 0
for point_id, group in final_net_df.groupby('id'):
    final_points_df.loc[final_points_df.id == point_id, "averageResidual"] = group['residualMag'].mean()

Exercise: What is the 95th percentile for the average residuals?

Plotting the points

Now that we have latitudes and longitudes for each point, we can generate some simple plots to look at them.

point_map = plt.figure(figsize=(10, 10))
point_ax = point_map.add_subplot(111)
point_ax.scatter(final_points_df["longitude"], final_points_df["latitude"], marker='+')
point_map.suptitle('Control Points')
point_ax.set_xlabel('Longitude')
point_ax.set_ylabel('Latitude')
plt.show()

It can also be helpful to color the points based on different values. The following cell draws the same plot but colors each point based on its average residual. Because the residuals are not uniformly distributed we also apply a lograithmic scale to the colors that you can see in the colorbar.

point_resid_map = plt.figure(figsize=(10, 10))
point_resid_ax = point_resid_map.add_subplot(111)
point_resid_norm = matplotlib.colors.LogNorm(vmax=final_points_df["averageResidual"].max())
point_resid_scatter = point_resid_ax.scatter(final_points_df["longitude"], final_points_df["latitude"], c=final_points_df["averageResidual"], alpha=0.5, norm=point_resid_norm, marker='+', cmap=plt.get_cmap('plasma'))
point_resid_cbar = plt.colorbar(point_resid_scatter)
point_resid_map.suptitle('Control Points')
point_resid_ax.set_xlabel('Longitude')
point_resid_ax.set_ylabel('Latitude')
point_resid_cbar.set_label('Average Residual Magnitude')
plt.show()

Plotting individual points can be helpful getting a general idea for the distribution of the points, but it can be hard to interpret the data in area where there are many points all ontop of each other. So, let’s combine near by points and determine the residual based on the region.

To do this, we’re going to bin the points into a regular grid across the latitude and longitude and then compute the mean within each bin.

Try changing the grid_step value and re-running the two cells

grid_step = 10

final_points_df['lonBin'] = final_points_df['longitude'].apply(lambda x: [e for e in range(-180, 180, grid_step) if e <= x][-1])
final_points_df['latBin'] = final_points_df['latitude'].apply(lambda x: [e for e in range(-90, 90, grid_step) if e <= x][-1])

avg_resid_binned = final_points_df.groupby(['lonBin', 'latBin'])['averageResidual'].mean()

filled_data = []
for lon_bin in range(-180, 180, grid_step):
    for lat_bin in range(-90, 90, grid_step):
        try:
            filled_data.append(avg_resid_binned.loc[lon_bin, lat_bin])
        except:
            filled_data.append(0)
filled_data = np.array(filled_data).reshape((int(360/grid_step), int(180/grid_step))).T

avg_gridded = plt.figure(figsize=(10, 5))
avg_gridded_ax = avg_gridded.add_subplot(111)
avg_gridded_plot = avg_gridded_ax.imshow(filled_data, origin='lower', extent= [-180, 180, -90, 90], cmap=plt.get_cmap('plasma'))
avg_gridded_ax.scatter(final_points_df["longitude"], final_points_df["latitude"], color='black', marker='+', alpha=0.1)
avg_gridded_cbar = plt.colorbar(avg_gridded_plot)
avg_gridded.suptitle('Average Residual by lat/lon grid')
avg_gridded_ax.set_xlabel('Longitude')
avg_gridded_ax.set_ylabel('Latitude')
avg_gridded_cbar.set_label('Average Residual Magnitude')
plt.show()

3D Plotting

2D plotting either requires these simple equal area projections or converting to another projection via pyproj. Instead, let’s look at our data in true 3D.

The following cell plots the same data as before but plots it in 3d instead of just a 2d projection

resid_fig_3d = plt.figure(figsize=(10, 10))
resid_ax_3d = resid_fig_3d.add_subplot(111, projection='3d')
resid_plot_3d = resid_ax_3d.scatter(final_points_df['adjustedX'], final_points_df['adjustedY'], final_points_df['adjustedZ'], c=final_points_df["averageResidual"], alpha=0.5, norm=point_resid_norm, marker='+', cmap=plt.get_cmap('plasma'))
resid_cbar_3d = plt.colorbar(resid_plot_3d)
resid_fig_3d.suptitle('3D Control Points')
resid_cbar_3d.set_label('Average Residual Magnitude (pix)')
plt.show()