Introduction to Python

Introduction

At the 2012 NCED Summer Institute in Earth Surface Dynamics (SIESD), we realized that there is a large community desire to transition to open source programming languages and tools, but that a roadblock is the lack of a coherent set of information to allow individual scientists to act on these desires without needing to invest a large amount of time. Here I will introduce programming in Python, a fairly intuitive open-source interpreted programming language, and some basic finite difference numerical methods through a simple example of a linear-diffusive hillslope.

Python

What is it?

Python is a general-purpose, high-level, interpreted programming language. This means that:

• General-purpose: You can use Python for numerical analysis, mapping and GIS, web design, interfacing with computer hardware, plotting and graphics, database management, image processing, and certainly lots of additional applications that scientists do not often have need to use. These applications are embodied in specific modules that can be imported into your Python script.
• High-level: Python uses array operations and pre-constructed functions to allow you to write your code in many fewer lines than would be required in a low-level language like C.
• Interpreted: Some languages use a compiler to convert what you write into bytecode that can be interpreted by the computer. Python is "interpreted", meaning that it calls pre-compiled chunks of code that are combined to run the program. It is easier and faster, but can be slower for computationally-intensive applications. However, your own compiled C code can be included in Python scripts for these applications, and the ease of use is a huge advantage!

Downloading and Installing

We will be using Python 2.7, the Scipy and Numpy packages that include tools for many common scientific and numerical applications and n-dimensional arrays, and the Matplotlib package that comprises a range of plotting tools.

I also highly recommend:

• IPython, a high-quality python shell (terminal window with syntax highlighting, tab-completion, etc.)
• A nice text editor that does syntax highlighting: gedit for any OS (though most commonly on Linux), notepad++ for Windows, and TextWrangler for Mac all seem good choices.

In addition, some colleagues of mine really like Spyder, which is a all-inclusive interface to edit and execute Python code, and is reminiscent of Matlab.

Here I will provide instructions to install these on the different major operating systems:

Linux

In Linux, the easiest way to install Python and its packages is with the package manager. This, as you probably already know, installs the required tools and all of their dependencies (programs that must be installed in order that the tool that you want to use works properly). I will be providing instructions for Ubuntu, but these should be easily adaptable to Debian and fairly similar to those for Red Hat / Fedora.

sudo apt-get install python2.7 python-numpy python-scipy python-matplotlib

You will have to enter your password for these programs to be installed.

I like the default gedit, but you may have your own favorite syntax-highlighting text editor (e.g., Kate).

To obtain IPython and Spyder, you type, respectively:

sudo apt-get install ipython
sudo apt-get install spyder

Nice and simple!

Helpful tools and information

• iPython
• Spyder

• Mathesaurus Matlab-Numpy

How to Write a Model in Python

Theory

Numerical models are fundamentally a computer-readable way to interpret rules that we have provided about how the natural world should work. Here, we are creating a simple model of a hillslope profile with rivers incising on either end. In our case, these rules take the form of partial differential equations. It is extremely important to have a full grasp of and—insofar as it is possible—analytical solution for these equations before beginning.

Here, we make the simple assumption that regolith discharge is a linear function of slope:

[math]\displaystyle{ Q = - k S }[/math]

Here, [math]\displaystyle{ Q }[/math] is regolith discharge [square meters per year], [math]\displaystyle{ S }[/math] is slope [unitless], and [math]\displaystyle{ k }[/math] is a constant of proportionality that scales slope to discharge. The negative sign is present because we want positive discharge where there are negative slopes (i.e. elevation decreasing as [math]\displaystyle{ x }[/math]-distance increases).

[math]\displaystyle{ Q = k S }[/math]

Source code

First, I will write the whole source code. Then I will break it down into sections where I explain each component. Then I will present a new model that uses an implicit matrix-style solution instead of a foward-time-stepping solution, and write about numerical (finite difference) methods.

#! /usr/bin/python

# Import modules
import numpy as np
from matplotlib import pyplot as plt

################################################################################
# hill.py
# 
# HILLSLOPE PROFILE
# 
# Center of hill is at 0, two incising channels at steady rate at edges
# no channel-hillslope feedback and assuming entire hill is made of mobile 
# regolith
# 
# This is a very simple example code in which variables are defined inline, 
# and the equations are solved by forward difference methods without 
# consideration for stability
# 
# Written by Andrew D. Wickert; August, 2013
################################################################################

# We will have the simple boundary condition of rivers on either end of the 
# hillslope that are incising at the same, steady rate.
zdot_channel = 2E-4 # [m], 0.5 mm/yr boundary incision

# time
dt = 100. # [years], the time step (from guessing, not from stability analysis)
t_final = 4E5 # [years], the final time for hillslope evolution
timesteps = np.arange(0,t_final + dt/10.,dt) # [years]

# Set up domain
dx = 10 # [m]
xmax = 500 # [m]
x = np.arange(-xmax, xmax + dx/10., dx) # [m], +dx/10. to make sure that edges are included
z = t_final * zdot_channel * np.ones(x.shape) # [m], elevation in meters - set such that the edges are 0 at t_final

# We're using a very simplified assumption that the rate of hillslope 
# material transport is linearly proportional to local slope, and that this 
# constant of proportionality is constant across the hill
k = 5E-3 # Slope -- soil flux scaling

# Loop through time and evolve elevation
for t in timesteps:
  # The channels set the boundaries
  z[0] -= zdot_channel * dt
  z[-1] -= zdot_channel * dt
  # for np.diff, out[n] = a[n+1] - a[n]
  # We are calculating the slopes between each of the cells, then using 
  # these to solve for discharge of material between interior cells.
  S = np.diff(z) # [-] slope
  Q = -k*S # [m**2/yr] discharge, goes downslope, hence the negative sign # POSITIVE - WHY? CANCELLED OUT AGAIN!
  # The change in internal elevation, due to conservation of mass, is equal 
  # once again to the x-defivative
  dzdt_interior = np.diff(Q)
  z[1:-1] -= dzdt_interior * dt
  # (It would have been quicker to take both derivatives at once, but would be 
  # somewhat less intuitive for descriptive purposes)

# Plot
plt.plot(x,z,'k-',linewidth=3)
plt.title('Final hillslope profile', fontsize=20, weight='bold')
plt.xlabel('Distance across hillside [m]', fontsize=16)
plt.ylabel('Elevation [m]', fontsize=16)
plt.xlim((x[0],x[-1]))
plt.show()

The "Shebang"

This script starts with:

#! /usr/bin/python

This is a "shebang" that tells the program where to look for the Python interpreter on Unix-like systems (the most popular being Mac and Linux, the former more popular for desktop applications and the latter more popular for supercomputers). You can find where Python is located by typing:

which python

at the Unix command prompt (i.e. terminal).