Unverified Commit 193056a6 authored by Jesse Mapel's avatar Jesse Mapel Committed by GitHub
Browse files

Moves DEM into an interface and adds the ability to use an ellipsoid (#576) 



* Moved DEM to an interface

* Actually added surface

* Added new options to demo config

* Updated Surface doc strings

Co-authored-by: default avatarjlaura <jlaura@usgs.gov>
parent 25427b05

autocnet/graph/network.py

+29 −25
Original line number Diff line number Diff line
@@ -55,6 +55,7 @@ from autocnet.vis.graph_view import plot_graph, cluster_plot 
from autocnet.control import control

from autocnet.spatial.overlap import compute_overlaps_sql

from autocnet.spatial.isis import point_info

from autocnet.spatial.surface import GdalDem, EllipsoidDem

from autocnet.transformation.spatial import reproject, og2oc



#np.warnings.filterwarnings('ignore')
@@ -1430,11 +1431,14 @@ class NetworkCandidateGraph(CandidateGraph): 


    def _setup_dem(self):

        spatial = self.config['spatial']

        semi_major = spatial.get('semimajor_rad')

        semi_minor = spatial.get('semiminor_rad')

        dem_type = spatial.get('dem_type')

        dem = spatial.get('dem', False)

        if dem:

            self.dem = GeoDataset(dem)

            self.dem = GdalDem(dem, semi_major, semi_minor, dem_type)

        else:

            self.dem = None

            self.dem = EllipsoidDem(semi_major, semi_minor)



    @property

    def Session(self):

autocnet/matcher/cross_instrument_matcher.py

+7 −8
Original line number Diff line number Diff line
@@ -163,8 +163,8 @@ def propagate_point(Session, 
                              'pgbouncer_port':6543,

                              'name':'somename'}



    dem       : plio.io.io_gdal.GeoDataset

                digital elevation model of target body

    dem       : surface

                surface model of target body



    lon       : np.float

                longitude of point you want to project
@@ -291,8 +291,7 @@ def propagate_point(Session, 
    if len(best_results[:,3])==1 and best_results[:,3][0] is None:

        return new_measures



    px, py = dem.latlon_to_pixel(lat, lon)

    height = dem.read_array(1, [px, py, 1, 1])[0][0]

    height = dem.get_height(lat, lon)



    semi_major = config['spatial']['semimajor_rad']

    semi_minor = config['spatial']['semiminor_rad']
@@ -350,8 +349,8 @@ def propagate_control_network(Session, 
                              'pgbouncer_port':6543,

                              'name':'somename'}



    dem       : plio.io.io_gdal.GeoDataset

                Digital elevation model of target body

    dem       : surface

                surface model of target body



    base_cnet : pd.DataFrame

                Dataframe representing the points you want to propagate. Must contain 'line', 'sample' location of

autocnet/matcher/ground.py

+35 −36
Original line number Diff line number Diff line
@@ -118,8 +118,7 @@ def propagate_ground_point(point, 
    dem = ncg.dem

    config = ncg.config



    px, py = dem.latlon_to_pixel(lat, lon)

    height = dem.read_array(1, [px, py, 1, 1])[0][0]

    height = dem.get_height(lat, lon)



    semi_major = config['spatial']['semimajor_rad']

    semi_minor = config['spatial']['semiminor_rad']

autocnet/spatial/overlap.py

+14 −18
Original line number Diff line number Diff line
@@ -185,8 +185,7 @@ def place_points_in_overlap(overlap, 


        # Calculate the height, the distance (in meters) above or

        # below the aeroid (meters above or below the BCBF spheroid).

        px, py = ncg.dem.latlon_to_pixel(lat, lon)

        height = ncg.dem.read_array(1, [px, py, 1, 1])[0][0]

        height = ncg.dem.get_height(lat, lon)



        # Need to get the first node and then convert from lat/lon to image space

        for reference_index, node in enumerate(nodes):
@@ -260,8 +259,7 @@ def place_points_in_overlap(overlap, 
                                                           semi_major, semi_minor, 'geocent', 'latlon')

            updated_lon, updated_lat = og2oc(updated_lon_og, updated_lat_og, semi_major, semi_minor)



            px, py = ncg.dem.latlon_to_pixel(updated_lat, updated_lon)

            updated_height = ncg.dem.read_array(1, [px, py, 1, 1])[0][0]

            updated_height = ncg.dem.get_height(updated_lat, updated_lon)





            # Get the BCEF coordinate from the lon, lat
@@ -403,8 +401,7 @@ def place_points_in_image(image, 


        # Calculate the height, the distance (in meters) above or

        # below the aeroid (meters above or below the BCBF spheroid).

        px, py = ncg.dem.latlon_to_pixel(lat, lon)

        height = ncg.dem.read_array(1, [px, py, 1, 1])[0][0]

        height = ncg.dem.get_height(lat, lon)



        with ncg.session_scope() as session:

            intersecting_images = session.query(Images.id, Images.path).filter(Images.geom.ST_Intersects(point_geometry)).all()
@@ -477,8 +474,7 @@ def place_points_in_image(image, 
                                                           semi_major, semi_minor, 'geocent', 'latlon')

            updated_lon, updated_lat = og2oc(updated_lon_og, updated_lat_og, semi_major, semi_minor)



            px, py = ncg.dem.latlon_to_pixel(updated_lat, updated_lon)

            updated_height = ncg.dem.read_array(1, [px, py, 1, 1])[0][0]

            updated_height = ncg.dem.get_height(updated_lat, updated_lon)





            # Get the BCEF coordinate from the lon, lat
@@ -583,4 +579,4 @@ def add_measures_to_point(pointid, cam_type='isis', ncg=None, Session=None): 
                    i += 1

            if i >= 2:

                point.ignore = False

    return

autocnet/spatial/surface.py

0 → 100644
+142 −0
Original line number Diff line number Diff line 
"""

A set of classes that represent the target surface. Each class implements the

get_height and get_radius functions for computing the height and radius respectively

at a given ground location (geocentric latitude and longitude).

"""



import numpy as np

from plio.io.io_gdal import GeoDataset



class EllipsoidDem:

    """

    A biaxial ellipsoid surface model.

    """



    def __init__(self, semi_major, semi_minor = None):

        """

        Create an ellipsoid DEM from a set of radii



        Parameters

        ----------

        semi_major : float

                     The equatorial semi-major radius of the ellipsoid.

        semi_minor : float

                     The polar semi-minor radius of the ellipsoid.

        """

        self.a = semi_major

        self.b = semi_major

        self.c = semi_major



        if semi_minor is not None:

            self.c = semi_minor



    def get_height(self, lat, lon):

        """

        Get the height above the ellipsoid at a ground location



        Parameters

        ----------

        lat : float

              The geocentric latitude in degrees

        lon : float

              The longitude in degrees

        """

        return 0



    def get_radius(self, lat, lon):

        """

        Get the radius at a ground location



        Parameters

        ----------

        lat : float

              The geocentric latitude in degrees

        lon : float

              The longitude in degrees

        """

        cos_lon = np.cos(np.deg2rad(lon))

        sin_lon = np.sin(np.deg2rad(lon))

        cos_lat = np.cos(np.deg2rad(lat))

        sin_lat = np.sin(np.deg2rad(lat))



        denom = self.b * self.b * cos_lon * cos_lon

        denom += self.a * self.a * sin_lon * sin_lon

        denom *= self.c * self.c * cos_lat * cos_lat

        denom += self.a * self.a * self.b * self.b * sin_lat * sin_lat



        return (self.a * self.b * self.c) / np.sqrt(denom)



class GdalDem(EllipsoidDem):

    """

    A raster DEM surface model.

    """



    def __init__(self, dem, semi_major, semi_minor = None, dem_type=None):

        """

        Create a GDAL dem from a dem file



        Parameters

        ----------

        dem : str

              The DEM file

        semi_major : float

                     The equatorial semi-major radius of the reference ellipsoid.

        semi_minor : float

                     The polar semi-minor radius of the reference ellipsoid.

        dem_type : str

                   The type of DEM, either height above reference ellipsoid or radius.

        """

        super().__init__(semi_major, semi_minor)

        dem_types = ('height', 'radius')

        if dem_type is None:

            dem_type = dem_types[0]

        if dem_type not in dem_types:

            raise ValueError(f'DEM type {dem_type} is not a valid option.')

        self.dem = GeoDataset(dem)

        self.dem_type = dem_type



    def get_raster_value(self, lat, lon):

        """

        Get the value of the dem raster at a ground location



        Parameters

        ----------

        lat : float

              The geocentric latitude in degrees

        lon : float

              The longitude in degrees

        """

        px, py = self.dem.latlon_to_pixel(lat, lon)

        return self.dem.read_array(1, [px, py, 1, 1])[0][0]



    def get_height(self, lat, lon):

        """

        Get the height above the ellipsoid at a ground location



        Parameters

        ----------

        lat : float

              The geocentric latitude in degrees

        lon : float

              The longitude in degrees

        """

        height = self.get_raster_value(lat, lon)

        if self.dem_type == 'radius':

            height -= super().get_radius(lat, lon)

        return height



    def get_radius(self, lat, lon):

        """

        Get the radius at a ground location



        Parameters

        ----------

        lat : float

              The geocentric latitude in degrees

        lon : float

              The longitude in degrees

        """

        radius = self.get_raster_value(lat, lon)

        if self.dem_type == 'height':

            radius += super().get_radius(lat, lon)

        return radius