diff --git a/environment.yml b/environment.yml index 8058563b69820a46285524b55a53d55d6afc211e..62986d87445aceef07f429afa33e9e4ce4806d6c 100644 --- a/environment.yml +++ b/environment.yml @@ -3,9 +3,11 @@ channels: - conda-forge - usgs-astrogeology dependencies: + - bokeh - coveralls - csmapi - gdal + - holoviews - jinja2 - jupyter - matplotlib diff --git a/knoten/vis.py b/knoten/vis.py new file mode 100644 index 0000000000000000000000000000000000000000..e9d07210ccf4f528d68a3243c081fad8ea3d02c2 --- /dev/null +++ b/knoten/vis.py @@ -0,0 +1,205 @@ +import tempfile + +import json +import pvl +import pyproj +import csmapi + +from knoten import csm + +from numbers import Number + +import numpy as np +import pandas as pd + +from pysis import isis +from pysis.exceptions import ProcessError + +import holoviews as hv +from holoviews import opts, dim +from bokeh.models import HoverTool + + +def reproject(record, semi_major, semi_minor, source_proj, dest_proj, **kwargs): + """ + Thin wrapper around PyProj's Transform() function to transform 1 or more three-dimensional + point from one coordinate system to another. If converting between Cartesian + body-centered body-fixed (BCBF) coordinates and Longitude/Latitude/Altitude coordinates, + the values input for semi-major and semi-minor axes determine whether latitudes are + planetographic or planetocentric and determine the shape of the datum for altitudes. + If semi_major == semi_minor, then latitudes are interpreted/created as planetocentric + and altitudes are interpreted/created as referenced to a spherical datum. + If semi_major != semi_minor, then latitudes are interpreted/created as planetographic + and altitudes are interpreted/created as referenced to an ellipsoidal datum. + + Parameters + ---------- + record : object + Pandas series object + + semi_major : float + Radius from the center of the body to the equater + + semi_minor : float + Radius from the pole to the center of mass + + source_proj : str + Pyproj string that defines a projection space ie. 'geocent' + + dest_proj : str + Pyproj string that defines a project space ie. 'latlon' + + Returns + ------- + : list + Transformed coordinates as y, x, z + + """ + source_pyproj = pyproj.Proj(proj = source_proj, a = semi_major, b = semi_minor) + dest_pyproj = pyproj.Proj(proj = dest_proj, a = semi_major, b = semi_minor) + + y, x, z = pyproj.transform(source_pyproj, dest_pyproj, record[0], record[1], record[2], **kwargs) + + return y, x, z + + +def point_info(cube_path, x, y, point_type, allow_outside=False): + """ + Use Isis's campt to get image/ground point info from an image + + Parameters + ---------- + cube_path : str + path to the input cube + + x : float + point in the x direction. Either a sample or a longitude value + depending on the point_type flag + + y : float + point in the y direction. Either a line or a latitude value + depending on the point_type flag + + point_type : str + Options: {"image", "ground"} + Pass "image" if x,y are in image space (sample, line) or + "ground" if in ground space (longitude, lattiude) + + Returns + ------- + : PvlObject + Pvl object containing campt returns + """ + point_type = point_type.lower() + + if point_type not in {"image", "ground"}: + raise Exception(f'{point_type} is not a valid point type, valid types are ["image", "ground"]') + + + if isinstance(x, Number) and isinstance(y, Number): + x, y = [x], [y] + + with tempfile.NamedTemporaryFile("w+") as f: + # ISIS wants points in a file, so write to a temp file + if point_type == "ground": + # campt uses lat, lon for ground but sample, line for image. + # So swap x,y for ground-to-image calls + x,y = y,x + elif point_type == "image": + # convert to ISIS pixels + x = np.add(x, .5) + y = np.add(y, .5) + + f.write("\n".join(["{}, {}".format(xval,yval) for xval,yval in zip(x, y)])) + f.flush() + try: + pvlres = isis.campt(from_=cube_path, coordlist=f.name, allowoutside=allow_outside, usecoordlist=True, coordtype=point_type) + except ProcessError as e: + warn(f"CAMPT call failed, image: {cube_path}\n{e.stderr}") + return + + pvlres = pvl.loads(pvlres) + if len(x) > 1 and len(y) > 1: + for r in pvlres: + # convert all pixels to PLIO pixels from ISIS + r[1]["Sample"] -= .5 + r[1]["Line"] -= .5 + else: + pvlres["GroundPoint"]["Sample"] -= .5 + pvlres["GroundPoint"]["Line"] -= .5 + + return pvlres + + +def reprojection_diff(isd, cube, nx=4, ny=8): + """ + """ + hv.extension('bokeh') + + isdjson = json.load(open(isd)) + + nlines = isdjson['image_lines'] + nsamples = isdjson['image_samples'] + + # generate meshgrid + xs, ys = np.mgrid[0:nsamples:nsamples/nx, 0:nlines:nlines/ny] + xs, ys = xs.flatten(), ys.flatten() + + csmcam = csm.create_csm(isd) + + # CS101 C++ programming class style dividers + ############################## + isis_pts = point_info(cube, xs, ys, 'image') + isisgnds = np.asarray([g[1]['BodyFixedCoordinate'].value for g in isis_pts]) + csm_pts = np.asarray([[p.samp, p.line] for p in [csmcam.groundToImage(csmapi.EcefCoord(*(np.asarray(bf)*1000))) for bf in isisgnds]]) + isis2csm_diff = csm_pts - np.asarray([xs,ys]).T + isis2csm_diffmag = np.linalg.norm(isis2csm_diff, axis=1) + isis2csm_angles = np.arctan2(*isis2csm_diff.T[::-1]) + + data = np.asarray([csm_pts.T[0], csm_pts.T[1], xs, ys, isis2csm_diff.T[0], isis2csm_diff.T[1], isis2csm_diffmag, isis2csm_angles]).T + data = pd.DataFrame(data, columns=['x', 'y', 'isisx','isisy', 'diffx', 'diffy', 'magnitude', 'angles']) + + isis2ground2csm_plot = hv.VectorField((data['x'], data['y'], data['angles'], data['magnitude']), group='isis2ground2csmimage' ).opts(opts.VectorField(pivot='tail', colorbar=True, cmap='coolwarm', title='ISIS2Ground->CSM2Image Pixel Diff', arrow_heads=True, magnitude='Magnitude', color=dim('Magnitude'))) + isis2ground2csm_plot = isis2ground2csm_plot.redim(x='sample', y='line') + isis2ground2csm_plot = isis2ground2csm_plot.opts(plot=dict(width=500, height=1000)) + isis2ground2csm_plot = isis2ground2csm_plot*hv.Points(data, group='isis2ground2csmimage').opts(size=5, tools=['hover'], invert_yaxis=True) + isis2ground2csm_plot = isis2ground2csm_plot.redim.range(a=(data['magnitude'].min(), data['magnitude'].max())) + ############################## + + ############################## + csmgnds = np.asarray([[p.x, p.y, p.z] for p in [csmcam.imageToGround(csmapi.ImageCoord(y,x), 0) for x,y in zip(xs,ys)]]) + csmlon, csmlat, _ = reproject(csmgnds.T, isdjson['radii']['semimajor'], isdjson['radii']['semiminor'], 'geocent', 'latlong') + isis_imgpts = point_info(cube, csmlon, csmlat, 'ground') + isis_imgpts = np.asarray([(p[1]['Sample'], p[1]['Line']) for p in isis_imgpts]) + + csm2isis_diff = isis_imgpts - np.asarray([xs,ys]).T + csm2isis_diffmag = np.linalg.norm(csm2isis_diff, axis=1) + csm2isis_angles = np.arctan2(*csm2isis_diff.T[::-1]) + + csm2isis_data = np.asarray([xs, ys, isis_imgpts.T[0], isis_imgpts.T[1], csm2isis_diff.T[0], csm2isis_diff.T[1], csm2isis_diffmag, csm2isis_angles]).T + csm2isis_data = pd.DataFrame(csm2isis_data, columns=['x', 'y', 'csmx','csmy', 'diffx', 'diffy', 'magnitude', 'angles']) + + csm2ground2isis_plot = hv.VectorField((csm2isis_data['x'], csm2isis_data['y'], csm2isis_data['angles'], csm2isis_data['magnitude']), group='csmground2image2isis').opts(opts.VectorField(pivot='tail', colorbar=True, cmap='coolwarm', title='CSM2Ground->ISIS2Image Pixel Diff', arrow_heads=True, magnitude='Magnitude', color=dim('Magnitude'))) + csm2ground2isis_plot = csm2ground2isis_plot.redim(x='sample', y='line') + csm2ground2isis_plot = csm2ground2isis_plot.opts(plot=dict(width=500, height=1000)) + csm2ground2isis_plot = csm2ground2isis_plot*hv.Points(csm2isis_data, group='csmground2image2isis').opts(size=5, tools=['hover'], invert_yaxis=True) + ############################### + + ############################### + isis_lonlat = np.asarray([[p[1]['PositiveEast360Longitude'].value, p[1]['PlanetocentricLatitude'].value] for p in isis_pts]) + csm_lonlat = np.asarray([csmlon+360, csmlat]).T + + isiscsm_difflatlon = isis_lonlat - csm_lonlat + isiscsm_difflatlonmag = np.linalg.norm(isiscsm_difflatlon, axis=1) + isiscsm_angleslatlon = np.arctan2(*isiscsm_difflatlon.T[::-1]) + + isiscsm_latlondata = np.asarray([isis_lonlat.T[0], isis_lonlat.T[1], csm_lonlat.T[0], csm_lonlat.T[1], isiscsm_difflatlon.T[0], isiscsm_difflatlon.T[1], isiscsm_difflatlonmag, isiscsm_angleslatlon]).T + isiscsm_latlondata = pd.DataFrame(isiscsm_latlondata, columns=['isislon', 'isislat', 'csmlon','csmlat', 'difflon', 'difflat', 'magnitude', 'angles']) + + isiscsm_plotlatlon = hv.VectorField((isiscsm_latlondata['isislon'], isiscsm_latlondata['isislat'], isiscsm_latlondata['angles'], isiscsm_latlondata['magnitude']), group='isisvscsmlatlon').opts(opts.VectorField(pivot='tail', colorbar=True, cmap='coolwarm', title='Image2Ground latlon Diff', arrow_heads=True, magnitude='Magnitude', color=dim('Magnitude'))) + isiscsm_plotlatlon = isiscsm_plotlatlon.redim(x='longitude', y='latitude') + isiscsm_plotlatlon = isiscsm_plotlatlon.opts(plot=dict(width=500, height=1000)) + isiscsm_plotlatlon = isiscsm_plotlatlon*hv.Points(isiscsm_latlondata, ['isislon', 'isislat'], group='isisvscsmlatlon').opts(size=5, tools=['hover'], invert_yaxis=True) + ############################### + + return isis2ground2csm_plot, csm2ground2isis_plot, isiscsm_plotlatlon, data, csm2isis_data, isiscsm_latlondata