# -*- coding: utf-8 -*-
"""
Created on Wed Apr 29 15:29:34 2020

@author: neil
"""

# this is a short code snippet that loads a DEM, and shows it as a 3D surface
# this is mainly to introduce several new commands and methods
# Warning, this needs a DEM file called 'manaslu.txt' in the local directory

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

f  = open('manaslu.txt')  # open a file 'descriptor', call it 'f', we can use f to read/write to a file
                          # NOTE! the file has to exist in the same directory as the called program
Z = np.loadtxt(f)         # read in an ASCII array as a DEM, np.loadtext() reads a file of numbers 
                          # and places them into an numpy array for later use.'Z' now has the DEM in a grid
f.close()                 # always good to close files when done

rows = np.size(Z,0)       # look at Z to get the shape of the DEM
cols = np.size(Z,1)       # np.size returns the length of the requested dimension of the array

fig = plt.figure()        # this next section allows us to look at the DEM
ax = fig.add_subplot(111, projection='3d')  # add a set of 3D axes for plotting a surface

# the following is a little obscure, but a common problem.  We have the grid of elevations, but we don't
# have the x,y positions to plot the z values.  First we make vectors of correct length in x and y
Y = np.linspace(0, 19800, rows)  # this produces a vector of x values, placed 100m apart
X = np.linspace(0, 19000, cols)  # the DEM is on a 100m spacing, this gives the real world DEM size

Xv,Yv = np.meshgrid(X, Y)        # meshgrid takes the X,Y vectors and makes 2 arrays of x and y values
                                 # so we can send x,y,z values are every grid point to a plotting routine
ax.plot_surface(-Xv,Yv,Z,cmap='terrain')  #  X is east west, Y is north south, this plots a 3D surface

plt.show()