diff --git a/knoten/csm.py b/knoten/csm.py
index 31df6b5f590f22b9f5f91ce37cfeb7f2dc7f0bd0..44fdbb458eae4f2d677624c0ff5b2a8be64f442e 100644
--- a/knoten/csm.py
+++ b/knoten/csm.py
@@ -9,6 +9,9 @@ import numpy as np
 import pyproj
 import requests
 import scipy.stats
+from singledispatch import singledispatch
+
+from plio.io.io_gdal import GeoDataset
 
 class NumpyEncoder(json.JSONEncoder):
     def default(self, obj):
@@ -94,6 +97,68 @@ def create_csm(image):
             if plugin.canModelBeConstructedFromISD(isd, model_name):
                 return plugin.constructModelFromISD(isd, model_name)
 
+@singledispatch
+def generate_ground_point(dem, image_pt, camera):
+    '''
+    Generates a longitude, latitude, and altitude coordinate for a given x, y
+    pixel coordinate
+
+    Parameters
+    ----------
+    dem : float or object
+          Either a float that represents the height above the datum or a
+          GeoDataset object generated from Plio off of a Digital Elevation
+          Model (DEM)
+    image_pt : tuple
+               Pair of x, y coordinates in pixel space
+    camera : object
+             USGSCSM camera model object
+    max_its : int, optional
+              Maximum number of iterations to go through if the height does not
+              converge
+    tolerance : float, optional
+                Number of decimal places to solve to when solving for height
+
+    Returns
+    -------
+    : object
+      CSM EcefCoord containing the newly computed lon, lat, and alt values
+      corresponding to the original image_pt coordinates
+    '''
+    if not isinstance(image_pt, csmapi.ImageCoord):
+        # Support a call where image_pt is in the form (x,y)
+        image_pt = csmapi.ImageCoord(*image_pt)
+
+    return camera.imageToGround(image_pt, dem)
+
+@generate_ground_point.register(GeoDataset)
+def _(dem, image_pt, camera, max_its = 20, tolerance = 0.001):
+    if not isinstance(image_pt, csmapi.ImageCoord):
+        # Support a call where image_pt is in the form (x,y)
+        image_pt = csmapi.ImageCoord(*image_pt)
+
+    intersection = generate_ground_point(0.0, image_pt, camera)
+    iterations = 0
+    semi_major, semi_minor = get_radii(camera)
+    ecef = pyproj.Proj(proj='geocent', a=semi_major, b=semi_minor)
+    lla = pyproj.Proj(proj='latlon', a=semi_major, b=semi_minor)
+    while iterations != max_its:
+        lon, lat, alt = pyproj.transform(ecef, lla, intersection.x, intersection.y, intersection.z)
+
+        px, py = dem.latlon_to_pixel(lon, lat)
+        height = dem.read_array(1, [px, py, 1, 1])[0][0]
+
+        next_intersection = camera.imageToGround(image_pt, float(height))
+        dist = max(abs(intersection.x - next_intersection.x),
+            abs(intersection.y - next_intersection.y),
+            abs(intersection.z - next_intersection.z))
+
+        intersection = next_intersection
+        iterations += 1
+        if dist < tolerance:
+            break
+    return intersection
+
 def generate_boundary(isize, npoints=10):
     '''
     Generates a bounding box given a camera model
@@ -116,7 +181,7 @@ def generate_boundary(isize, npoints=10):
 
     return boundary
 
-def generate_latlon_boundary(camera, boundary, radii=None):
+def generate_latlon_boundary(camera, boundary, dem=0.0, radii=None, **kwargs):
     '''
     Generates a latlon bounding box given a camera model
 
@@ -129,7 +194,7 @@ def generate_latlon_boundary(camera, boundary, radii=None):
     radii : tuple
             in the form (semimajor, semiminor) axes in meters. The default
             None, will attempt to get the radii from the camera object.
- 
+
     Returns
     -------
     lons : list
@@ -139,6 +204,7 @@ def generate_latlon_boundary(camera, boundary, radii=None):
     alts : list
            List of altitude values
     '''
+
     if radii is None:
         semi_major, semi_minor = get_radii(camera)
     else:
@@ -152,7 +218,7 @@ def generate_latlon_boundary(camera, boundary, radii=None):
     for i, b in enumerate(boundary):
         # Could be potential errors or warnings from imageToGround
         try:
-            gnd = camera.imageToGround(csmapi.ImageCoord(*b), 0)
+            gnd = generate_ground_point(dem, b, camera, **kwargs)
         except:
             pass
 
@@ -171,15 +237,6 @@ def generate_gcps(camera, boundary, radii=None):
              csmapi generated camera model
     boundary : list
                of image boundary coordinates
-    nnodes : int
-             Not sure
-    semi_major : int
-                 Semimajor axis of the target body
-    semi_minor : int
-                 Semiminor axis of the target body
-    n_points : int
-               Number of points to generate between the corners of the bounding
-               box per side.
     Returns
     -------
     gcps : list
@@ -203,22 +260,13 @@ def generate_gcps(camera, boundary, radii=None):
 
     return gcps
 
-def generate_latlon_footprint(camera, boundary, radii=None):
+def generate_latlon_footprint(camera, boundary, dem=0.0, radii=None, **kwargs):
     '''
     Generates a latlon footprint from a csmapi generated camera model
     Parameters
     ----------
     camera : object
              csmapi generated camera model
-    nnodes : int
-             Not sure
-    semi_major : int
-                 Semimajor axis of the target body
-    semi_minor : int
-                 Semiminor axis of the target body
-    n_points : int
-               Number of points to generate between the corners of the bounding
-               box per side.
     Returns
     -------
     : object
@@ -229,7 +277,7 @@ def generate_latlon_footprint(camera, boundary, radii=None):
     else:
         semi_major, semi_minor = radii
 
-    lons, lats, _ = generate_latlon_boundary(camera, boundary, radii=radii)
+    lons, lats, _ = generate_latlon_boundary(camera, boundary, dem=dem, radii=radii, **kwargs)
 
     # Transform coords from -180, 180 to 0, 360
     # Makes crossing the meridian easier to identify
@@ -305,7 +353,7 @@ def generate_bodyfixed_footprint(camera, boundary, radii=None):
         semi_major, semi_minor = get_radii(camera)
     else:
         semi_major, semi_minor = radii
-        
+
     latlon_fp = generate_latlon_footprint(camera, boundary, radii=radii)
 
     ecef = pyproj.Proj(proj='geocent', a=semi_major, b=semi_minor)