map elevation data

import gdal
d = gdal.Open('10n060e_20101117_gmted_mea300.tif')
arr = d.ReadAsArray()

def get_elevation_at_point(x, y, elevation_array):
  # skip interpolation for now
  return elevation_array[int(x), int(y)]

def get_elevation_profile_for_path(x, y, elevation_array):
  e = [get_elevation_at_point(x[i], y[i], elevation_array) for i in range(len(x))]
  return np.array(e)
  
def get_slope_profile_for_path(x,y,elevation_array):
  gradient = np.gradient(elevation_array)
  sx = [gradient[0][x[i],y[i]] for i in range(len(x))]
  sy = [gradient[1][x[i],y[i]] for i in range(len(x))]
  return np.array(sx) + np.array(sy) # not quite right - should separate the gradient in the direction of travel from the traverse angle
  
# find countours and ridges, link contours and ridges together to minimize total path length
# ridges are basically edges - try Canny edge detector
# contours - level sets - try an iterated findContours
# then connect starting and finishing point using contours and ridges
#    subject to additional constraints: e.g. max slope when travelling on a contour
# these contours / ridges are essentially like a road network. learn about the algorithms used to find efficient paths on a road network
# connections between contours and ridges are nodes
# use shortest path algorithms to find possible routes, e.g. A*
# http://www.redblobgames.com/pathfinding/a-star/introduction.html
# http://www.redblobgames.com/pathfinding/a-star/implementation.html

# could simply use A* with changes in elevation as a cost, but this would be very expensive for large high resolution maps. 
# Finding the contours / ridges simplifies the search problem

# A* package in networkx:
# https://networkx.github.io/documentation/networkx-1.10/reference/algorithms.shortest_paths.html