Visualization, topography, and DEMs

1
Visualization, topography, and DEMs

• How to make pretty maps

2
How high is that mountain?

• What is elevation?
• Measured from what?
• Sea Level

3
Sea level, the ellipsoid, and the geoid

• from sea level (high tide or low tide?)
• okay, mean sea level (an average)
• suppose the earth was completely covered in water
  -
• what shape would it be ?
• A. A sphere (rough approximation).
• B. An ellipsoid (the ellipsoid)(a better
  approximation).
• C. A very slightly lumpy ellipsoid (the
  geoid)(the best so far).
• So now we have three surfaces
• Geoid where the mean sea surface would be.
• Ellipsoid a geometric shape defined by an
  mathematical equation that approximates the
  geoid.
• Topography the actual surface of the earth.

4
Reference ellipsoid

• Geometric shape that closely matches the shape of
  the earth.
• Several different ones have been developed and
  optimally fit only one area of the world.
• An offset (a datum) is often applied for specific
  areas.
• Some common ellipsoids and datums are listed
  below.

Name (date) Semi-major axisflattening datum W
orld Geodetic System 1984 (WGS84) 6378137
298.257223563 Geodetic Reference System 1980
(GRS80) 6378137 298.257222101 World Geodetic
System 1984 World Geodetic System 1972
(WGS72) 6378135 298.26 World Geodetic System
1972 Clark (1866) 6378206 294.98 North
American Datum 1927
5
Geoid and gravity

• Problem we dont know this shape exactly.
• Pretty close to the ellipsoid and defined by
  deviation from the ellipsoid.
• Depends on gravity variations and hence density
  variations.
• Surface perpendicular to a vertical line (plumb
  line).

The equipotential surface of the Earth's gravity
field which best fits, in a least squares sense,
global mean sea level (U.S. National Geodetic
Survey definition).
6
So what?

• If we want to compare measurements, we need to
  use the same system.
• Ignoring it will lead to mostly small errors.
• Most GIS and software programs ask for the
  reference ellipsoid or datum.
• Handhelp GPS units usually use WGS84.
• Many topo maps we use are based on NAD27 (and the
  Clarke ellipsoid).
• Can be off by tens of m.
• Most pre-GPS topo maps are referenced to local
  sea level so are essentially referenced to the
  geoid.

7
Elevation (and other data) maps

• Want to represent 3D data on a 2D surface.
• Contour
• Shaded
• Elevation
• Slope (or gradient)
• 3D visualizations.

8
Contours

• Lines that connect points of equal elevation.
• Do not split or cross other contour lines.
• Closely spaced lines indicate slope widely
  spaced lines indicate flatter areas.
• Topography maps generally have smooth curves
  reflecting smooth changes.
• Faults in subsurface maps often have
  discontinuities.

9
Creating contour maps

• An interpretation of data.
• Does not require evenly spaced data.
• By hand
• Can be better easily to include a priori
  information
• Impractical for large datasets
• By computer
• Often grid the data first
• May create artifacts
• Edge effects
• Misses control points
• Bulls eye
• Need to understand data and method

10
Hand contouring

• Various methods
• Place contours divided linearly between points
  (basically linear interpolation)
• Parallel draw contours parallel to each other.
• Equal-space assume uniform slope over all areas.
• Interpretative do anything you want as long as
  the data is honored.

11
Computer contouring

• Delaunay triangulation
• Creates a regular grid (or surface) from
  irregularly spaced data.
• Various methods
• Linear interpolation
• Weighted interpolation
• Kriging and geostatistics (may be best for
  geology)
• Advantages of grids
• Even data distribution
• Allows easy application of mathematical
  operations
• Filtering
• Smoothing
• Display

12
Example a set of numbers

• 0 0 0 2 10 22
• 0 0 0 3 15 26
• 0 0 3 10 21 32
• 0 0 4 13 22 29
• 0 0 0 7 17 26

hand
computer
13
Color shading
14
3D perspective
15
Perspective with shading
16
comparison
map good for reference, bad for visualization
perspective good for visualization, bad
for reference
17
More fancy
18
Volume visualization
19
DEM

• Digital elevation models
• Gridded elevations
• Not always accurate (thats why it is a model)
• Getting more widely available.
• Useful for
• Shaded and perspective maps
• Hydrological measurements

20
Some useful DEM bathymetry

• National Elevation Data set (NED)
• US 30 m resolution (except Alaska)
• Seamless http//seamless.usgs.gov/
• Shuttle Radar Topography (SRTM)
• US 30 m http//seamless.usgs.gov/
• Worldwide 90m (between 60 N and 60S)
• Some gaps and problems in high slopes and lakes.
• World topo bathymetry etopo2http//www.ngdc.no
  aa.gov/mgg/

21
Miscellaneous DEM DMM

• Mars (at 1 km resolution)
• High resolution offshore CA bathymetry.http//wrgi
  s.wr.usgs.gov/dds/dds-55/pacmaps/site.htm

22
Matlab
23
Matlab
Software package to do math Mainly discrete
(numbers) but also includes symbolic (as in
algebra) Includes scripting capability Fair
graphics capability Currently has lots of
toolboxes to do almost anything Currently is
getting expensive Student license
available Another similar package is Octave
free, slower, not as many functions but most of
the widely used ones
24
Command window
Simple operations addition - subtraction / divi
sion multiplication power Also try the
help command
25
Can also assign variable with an equal sign
variables
x 1 means put the value 1 into the variable
x y 2 puts the value 2 into the variable
y x y now yields 3 Adding a semicolon
stops matlab from echoing the answer or
assignment. Use up arrow to go to previous
commands
history
26
Variables can hold more than 1 number either a
vector or a matrix. The operations can apply to
matrices as well. This saves a lot of time. We
could also assign a vector gtgt x
00.110 Count to 10 by 0.1 Put it all in x
27
Graphics 2D and 3D
28
Radioactive decay

• Radioactive elements spontaneously and randomly
  emit particles.
• The rate of decay differs for different elements
  and isotopes.
• Carbon 14 decays relatively quickly potassium
  slowly.
• Parent elements decay to daughter elements.

29
Decay process

• Exponential process so that half the parent
  element decays in a certain time span.
• Controlled by several equations.
• N/N0exp(-kt)
• N is number of atoms after a time t.
• N0 is original number of atoms.
• K is the decay constant (different for each
  element/isotope)
• Exp is exponent using e
• k ln(2)/t1/2
• ln is natural log
• The number of disintegrations depend on the
  number of parents.
• Carbon 14 has a large decay constant potassium
  has a small one.

30
Where did these equations come from?

• Suppose we let N(t) denote the amount of material
  at a given time t.
• The change in N (dN) over a time (dt) depends on
  N.
• dN/dt -kN
• This is an ordinary differential equation (ODE).
• A couple of ways to solve this
• Take a class in ODE (may be bad for the GPA and
  you will forget it anyway after a while)
• Get a book or look it up on google
• Use matlab

31
Using Matlab symbolic calculator

• gtgt A dsolve('DN -kN')
• A
• C1exp(-kt)
• Which is pretty much what we said before.

32
Now for k

• K is the time for half the element to decay
• N/N0 exp(-kt)
• Half gone means N/N0 ½ so ½ exp(-kt)
• Take the natural log of both sides
• ln(1/2) -kt
• K ln(2)/t remember ln(1/2) - ln(2)
• In Matlab
• gtgt syms k t
• gtgt k solve('1/2 exp(-1kt)')
• k log(2)/k

33
Other symbolic tricks

• solve solves normal equationsdsolve does
  differential equations
• x solve('ax2 bx c 0')
• Simplify makes things simpler
• Try f cos(x)2 sin(x)2
• f simple(f)
• What should f equal?

34
Linear algebra

• Handles matrices
• gtgt A 1,2,33,4,68,6,2
• A 1 2 3
• 3 4 6
• 8 6 2
• gtgt B 3,2,83,1,69,6,2
• B 3 2 8
• 3 1 6
• 9 6 2
• gtgt AB
• ans
• 36 22 26
• 75 46 60
• 60 34 104
• gt inv(A)C
• ans
• 3.0000 2.0000 8.0000
• 3.0000 1.0000 6.0000
• 9.0000 6.0000 2.0000

35
Programming constructs

• Loops
• gtgt for i 13
• i
• end
• If - then statements
• if I J
• A(I,J) 2
• else
• A(I,J) 0
• end
• Functions
• -------------------------
• function mean avg(x,n)
• MEAN subfunction
• mean sum(x)/n

36
The Theis equation

• Extracting water from a well will lower the
  aquifer surface.
• It does not lower evenly (like in a bathtub)
• It forms a cone of depression around the well.
• The depth and radius of this depressed area
  depends on the rate of extraction and the
  properties of the aquifer.

well
37
Want to predict drawdown
original surface
H0-h
b
h
H0
r
Q pumping rate H hydraulic head H0 initial
hydraulic head H0-h drawdown T aquifer
transmissivity t time since pumping
began r radial distance from well S aquifer
storativity
H0-h (Q/rpT)W(u) U (r2S)/(4Tt) W(u)
-0.5772-ln(u) u u2/22u3/33-u4/44
38
strategy

• Understand the problem
• Start simple
• Proceed in steps
• Do 1 part test
• Do 2 part test
• And so on
• If we try to do it all at once there are a lot of
  places that can go wrong.

39
functions

• We will make a new function

---------------------------------- the well
function function w well(u) the well
function calculate out to 8 log is the
natural log w -0.5772 - log(u) u -
(u2)/(22) (u3)/(33) - (u4)/(44)
(u5)/(55) - (u6)/(66) (u7)/(77) -
(u8)/(88) -----------------------------------
40
Test with number from book
gtgt well(10(-9)) ans 20.1461 gtgt
Fetter (2001) says 20.15 so it looks like it is
working okay (usually try a few number to make
sure)
41
Now do the rest
The Theis function Q pumping rate drawdown
is drawdown (h0-h) r radial distance from
well S storativity trans transmissivity
time time since pumping began function
u,drawdown theis(r,Q,S,trans,time) u
(r2)S/(4transtime) drawdown
(Q/(4.03.14159trans))well(u) ----------------
------------------
Uses our previously defined well function
42
Final check
gtgt r 7 gtgt Q 2725 gtgt S 0.0051 gtgt trans
299 gtgt time 1 gtgt u,drawdown
theis(r,Q,S,trans,time) u 2.0895e-004 drawdown
5.7269 gtgt
Which matches Fetters example problem
43
Now calculate over a grid and vary r (the radial
distance)
calculate over a 200 m grid with 2 meter
spacing well is at grid point (100,100) r
7 Q 2725 S 0.0051 trans 299 time 1 for i
1101 for j 1101 x 2i
y 2j r sqrt((100-x)2
(100-y)2) u,drawdown
theis(r,Q,S,trans,time) ddgrid(i,j)
-1drawdown end end
44
We get