Transcript Spatial Analysis Using Grids - UNL School of Natural Resources

Handling Gridded Data:
Topography and Projections
GIS in Water Resources
Fall 2014
by Ayse Kilic with materials from
David G. Tarboton, Utah State University and from ESRI
software
Learning Objectives
• Calculation of slope on a raster using
– ArcGIS method based in finite differences
– D8 steepest single flow direction
– D steepest outward slope on grid centered triangular facets
• Map Projections
– State Plane Coordinate System
– UTM (Universal Transverse Mercator Coordinate System)
Spatial Surfaces used in Hydrology
Elevation Surface — the ground surface elevation at
each point -- Expressed as a Digital Elevation Model
for Gridded Data
Types of Elevation Data available
Data
Spatial reference
Pixel size
•GCS_WGS_1984
•Decimal degrees
(Global
Topography) •WGS 1984
30 arcsec
(1 km)
GTOPO
SRTM
(Shuttle Radar
Topography
Mission)
NED 30
(National
Elevation
Data)
NED 10
Z
Bit Depth
units
16-bit
signed/un
m
signed
Integer
•GCS_WGS_1984
•Decimal degrees
•WGS 1984
90 m
m
16-bit
signed
Integer
•GDC_North_America_1983
•Decimal degrees
•NAD 1983
1 arcsec
(30 m)
m
Float
•GDC_North_America_1983
•Decimal degrees
•NAD 1983
1/3 arcsec
m
(10 m)
Float
3 ft
Float
•NAD83_HARN_StatePlane_Oregon_North
Lidar
•Foot
(DEM/DSM)
•NAD 1983 HARN
ft
Slope Handout
http://snr.unl.edu/kilic/giswr/2014/Slope.pdf
Determine the length, slope and azimuth of the line AB.
3-D detail of the Tongue river at the WY/Montana border from LIDAR.
LIDAR from aircraft or from the
ground can provide amazing detail
on elevation, including individual
tree heights and hydraulic channels
Roberto Gutierrez
University of Texas at Austin
Topographic Slope
• Used to determine how (quickly) water flows downhill
and concentrates into streams
• Topographic slope can be determined from a DEM
Topographic Slope
• There are three alternative sets of inputs (choose one)
– Surface derivative z (dz/dx, dz/dy)
– Vector with x and y components (Sx, Sy). Slope in x and y direction.
– Vector with magnitude (slope) and direction (aspect) (S, )
ArcGIS “Slope” tool
• Calculates the maximum rate of change in
value from that cell to its neighbors
• Calculates for each cell
• Represents the rate of change of elevation for
each digital elevation model (DEM) cell
(slope). Slope is the first derivative of a DEM
• The lower the slope value, the flatter the
terrain; the higher the slope value, the steeper
the terrain.
Definition of X, Y, and Z in 3D space
Z axis is the direction that elevation changes (up or down)
Origin is the location of the point of interest (pixel or grid cell)
X axis is the direction that
X has a changing value
(East-West in ArcGIS)
Y axis is the direction that Y
has a changing value
(North-South in ArcGIS)
X, and Y are horizontal distances
Z is the vertical distance
The X, Y, Z axes are at right angles to one another
Definition of Slope
𝑆𝑙𝑜𝑝𝑒 𝑖𝑛 𝑑𝑒𝑔𝑟𝑒𝑒𝑠(𝛳) =
𝑅𝑖𝑠𝑒
𝑇𝑎𝑛 (𝛳) =
𝑅𝑢𝑛
𝑃𝑒𝑟𝑐𝑒𝑛𝑡 𝑜𝑓 𝑆𝑙𝑜𝑝𝑒(𝛳) =
𝑅𝑖𝑠𝑒
𝑅𝑢𝑛
𝑅𝑖𝑠𝑒
𝑅𝑢𝑛
(0)
* 100 (%)
Run is the horizontal distance calculated using X and Y
Rise is the vertical distance calculated using Z (elevation)
Slope ranges (-900, +900) or (-infinity %, +infinity %)
How Slope works
𝑆𝑙𝑜𝑝𝑒(𝛳) =
𝑆𝑙𝑜𝑝𝑒(𝛳) =
𝑅𝑖𝑠𝑒
𝑅𝑢𝑛
𝑅𝑖𝑠𝑒
𝑅𝑢𝑛
(0)
* 100 (%)
When the angle (𝛳) is 45 degrees, the rise is equal to the run,
and the percent rise is 100 percent
When the slope angle (𝛳 ) approaches vertical (90
degrees), the percent rise (slope) begins to approach
infinity.
Run is the horizontal distance calculated using X and Y
Rise is the vertical distance calculated using Z (elevation)
Pythagorean theorem
Used to calculate Run where a =ΔY and b = ΔX
The calculated “c” is the “Run”
ArcGIS “Slope” tool
y
a
b
c
d
e
f
g
h
i
Calculates slope for each cell. In this illustration, it is for Cell “e”
For each cell, the Slope tool calculates the maximum of the rate
of change in value from that cell to each of its eight neighbors
The rate of change in the x direction for cell e is
calculated with the following algorithm
x
dz
a + 2d + g − c + 2f + i
=−
dx
8∆
The rate of change in the y direction
for cell e is calculated with the
following algorithm
∆
b
a
d
c
e
g
y
f
h
i
2∆
dz
g + 2h + i − a + 2b + c
=−
dy
8∆
The negative sign in front of the
equations is because x increases to
the right (east) and y increases to the
x
north. Now dz/dx is + if z increases
with increasing x.
ArcGIS “Slope” tool
The two equations for dz/dx and dz/dy are simplified from the first equation below.
The basis for that equation is illustrated in the Figure and represents an average of
central finite differences over each of the three rows of cells, with the middle row
counting twice as it appears in averages on each side.
𝑎−𝑐
2∆
𝑑−𝑓
2∆
c
∆
b
a
f
e
d
=-
dz
a + 2d + g − c + 2f + i
=−
dx
8∆
𝑔−𝑖
2∆
i
dz
g + 2h + i − a + 2b + c
=−
dy
8∆
h
g
y
dz
dx
2∆
𝑎−𝑐 𝑑−𝑓 𝑑−𝑓 𝑔−𝑖
+
+
2∆
2∆ + 2∆
2∆
2
2
2
x
m
Slope  (dz / dx)  (dz / dy ) 
m
2
2
The negative sign in front of the equations is because we are computing uphill slope
Definition of Azimuth
𝑦
Y axis is the
direction that Y has
a changing value Δ𝑦
(South to North in
ArcGIS)
Δ𝑥
𝛼
𝛼 = Azimuth, angle
defined as degrees
clockwise from North
𝑥
This is my grid cell location
X axis is the direction that X
has a changing value (West to
East in ArcGIS)
Definition of Azimuth
𝑦
Δ𝑦
Δ𝑥
𝛼
∆𝑋
𝑇𝑎𝑛 (𝛼) =
∆𝑌
Solve for α by Inverting the Tangent
Function (ArcTan)
𝛼 = Azimuth, angle defined as degrees
clockwise from North
𝑥
∆𝑋
𝛼 = 𝐴𝑟𝑐𝑇𝑎𝑛
= Azimuth
∆𝑌
The other way to write ArcTan is Tan-1
Azimuth=
Convert from radians to degrees (180/π)
Azimuth is the angle between North and any desired direction you want to travel
ArcGIS Aspect – the steepest downslope direction
If I pour water on the ground, which direction does it flow?
Aspect is the azimuth associated with the steepest downhill
slope.
Therefore, we use slopes instead of distances in the tangent
function.
In Arc, with grid cells it is easiest to calculate Aspect using the
ratio of slopes (dz/dx) and (dz/dy).
dz
dy
 dz / dx 
atan

 dz / dy 
𝛼
dz
dx
= Aspect
Example for topographic slope
Mesh spacing=30 m
Slope/Aspect at cell e?
30
a
b
80
d
c
74
e
69
g
67
h
60
52
63
dz
(a  2d  g) - (c  2f  i)

dx
8 * x_mesh_spacing
o
f 145.2
(80  2 * 69  60)  (63  2 * 56  48)


 0.229 m/m
56
8 * 30
dz
(g  2h  i) - (a  2b  c)

i
dy
8 * y_mesh_spacing
48
(60  2 * 52  48)  (80  2 * 74  63)

 0.329 m/m
8 * 30
Slope  (dz / dx) 2  (dz / dy ) 2  (0.229) 2  0.329 2
 0.401 m/m
Slope  atan (0.401)
Note that this is the slope in Uphill direction (it is a positive number)
180

 21.8o
Converts slope from m/m to degrees (180/π)
Example for Aspect
Mesh spacing=30 m
Aspect at cell e?
30
a
b
o
-34.8
80
74
d
e
69
g
63
f 145.2o
67
h
60
c
56
i
52
48
dz
(a  2d  g) - (c  2f  i)

dx
8 * x_mesh_spacing
(80  2 * 69  60)  (63  2 * 56  48)

 0.229 m/m
8 * 30
dz
(g  2h  i) - (a  2b  c)

dy
8 * y_mesh_spacing
(60  2 * 52  48)  (80  2 * 74  63)

 0.329 m/m
8 * 30
 dz / dx 
  0.229  180
  atan
Aspect  atan
 34.8o

 0.329  
 dz / dy 
One more adjustment: The above Aspect is in the direction of increasing elevation (increasing dz/dx).
We need to add 180o to this calculated aspect to get the direction of decreasing z (i.e., the steepest
downhill slope)
 34.8  180o  145.2o
The Atan function is multivalued on the full circle and only unique in a range of
180 degrees. To unambiguously determine the direction from two components
you really need the atan2 function that keeps the sign on y and x components
separately.
For example, let y = y component of a vector
x = x component of a vector
atan(x/y) gives the direction of the vector as an angle (with the ratio x/y since
angle here is measured from north).
But x/y is the same value if y is positive and x negative, or x positive and y
negative. So once you take the ratio x/y, if you get a negative number you do not
know which (y or x) was negative.
A way to resolve this is
angle = atan(x/y)
if(0 < angle < 180 and dz/dx < 0) then
aspect = angle + 180 (flip the direction because dz/dx is negative)
else
aspect = angle
endif
D8 steepest single flow direction
(Eight Direction Pour Point Model)
In a gridded system, water can only flow to one of eight adjacent cells
32
64
16
8
4
128
The direction of flow is determined
by the direction of steepest descent:
1
Maximum_drop = (change_in_z-value
/ distance) * 100
2
This is maximum percentage drop.
Defined as “Hydrologic slope” in
ArcGIS
ESRI Direction encoding (ArcGIS)
Hydrologic Slope (Flow Direction Tool)
Find Direction of Steepest Descent (ArcGIS)
30
30
80
74
63
80
74
63
69
67
56
69
67
56
60
52
48
60
52
48
67  52
 0.50
Slope:
30
67  48
Slope:
 0.45
30 2
For diagonal direction, the denominator
for slope includes square root of 2
Limitation due to 8 grid directions.
?
The true flow direction follows the red arrow. However, we can only choose one of the blue
arrows because we have to use one of eight adjacent cells.
The D Algorithm
Steepest direction
downslope
Proportion flowing to
neighboring grid cell 3
is 2/(1+2)
Proportion
flowing to
neighboring
grid cell 4 is
1/(1+2)
3
4
2
1
2
Flow
direction.
5
1
6
8
7
Tarboton, D. G., (1997), "A New Method for the Determination of Flow Directions and
Contributing Areas in Grid Digital Elevation Models," Water Resources Research,
33(2): 309-319.) (http://www.engineering.usu.edu/cee/faculty/dtarb/dinf.pdf)
The D Algorithm
Steepest direction
downslope
4
z3 3
z2
2
2
0
zo
1
z1
1
5

 z1  z 2 

1  atan
 z0  z1 
 z1  z 2   z0  z1 
S 

 
     
2
6
8
7
2
If 1 does not fit within the triangle, the angle is chosen along the steepest edge or diagonal
resulting in a slope and direction equivalent to D8
D∞ Example
30
z4
80
z5
z3
z2
74
63
zo
69
67
z6
z7
60
52
56z1
z8
48
 z 7  z8 

1  atan
 z0  z7 
 52  48 
o
 atan
  14.9
 67  52 
z z  z z 
S   1 2   0 1 
     
2
2
 52  48   67  52 
S 
 

 30   30 
 0.517
2
284.9o
14.9o
The tool is available at
http://hydrology.usu.edu/taudem/taudem5/documentation.html
2
Automating Processes using Model Builder
Using a DEM tif file as input
Elevation (m) for Upper Klamath Lake Basin, OR
Elevation Contours for Wood River Valley
Watershed of Upper Klamath Lake Basin
Slope (%) for Upper Klamath Lake Basin, OR
(-infinity, + infinity)
Slope (Degree) for Upper Klamath Lake
Basin, OR
Aspect (Degree) for Upper Klamath
Lake Basin, OR
Percentage Drop (Degree) for Upper Klamath
Flow Direction
Integer raster whose values range from 1 to 255
32
64
16
8
128
1
4
2
Hillshade
Hypothetical illumination of a surface by determining illumination values for each cell in a
raster. It does this by setting a position for a hypothetical light source and calculating the
illumination values of each cell in relation to neighboring cells.
Azimuth
The azimuth is the angular direction of the sun, measured
from north in clockwise degrees from 0 to 360. An
azimuth of 90º is east. The default azimuth is 315º (NW).
Altitude
The altitude is the slope or angle of the illumination
source above the horizon. The units are in degrees, from 0
(on the horizon) to 90 (overhead). The default is 45
degrees.
Hillshade
Map Projection Parameters
State Plane Coordinate Systems
• Like UTM but customized to minimize error
within States-1930
• NAD27 – feet
• NAD83 – feet or meters
39
State Plane Coordinate Systems
•
•
•
•
TM for North South oriented states
Lambert Conformal Conic- East West orientation
Depending on TM or LCC, zones baselines and meridians positioned differently
Baselines are placed under the zones Southern Border and measured North from
there
40
State Plane zones
State plane zones of Minnesota (1/6 of the zone width) and details of standard
parallel placement (Lambert Conformal Conic)
State Plane Coordinate System-overview
• Used by state/local governments
– Units are feet (now appearing in meters)
• Panhandle of Alaska: oblique Mercator (at angle)
• N/S states (Missouri): Transverse Mercator projection
• E/W states: Lambert Conformal Conic for an East-West
State (California, Nebraska)
• Doesn’t make sense?
– States divided into smaller zones
– Distortion is extremely small for local area
• Coordinates: SW corner is 0, 0 .
• Mapping large areas (over >1 zone) is trouble!
What are the coordinates of the origin (fo, lo) and
the corresponding (Xo, Yo) ?
(fo, lo) = (39° 50’ N, 100° W)
(lo, fo) = (100° W, 39° 50’ N)
(Xo, Yo) = (1640416.667, 0.000)
Don’t get confused.
“X” is always
associated with
Longitude, “Y” with
Latitude
Zones in the Texas State Plane
Universal Transverse Mercator Coordinate System
• Line of intersection at a central meridian.
• Distances measured with respect to the central
meridian and from the equator in meters
• 60 zones around the Earth from East to West
• False Northing and False Easting
46
Universal Transverse Mercator Coordinate System
• False Easting to avoid negative values: Central
meridian is denoted as 500,000 not “0”
• Measurements are shown as positive
• So if a point is 11,254 meters west of the
central meridian, it is shown as 500,00011,254 or 488,746.
• False northing. In the southern hemisphere,
10,000,000 is the value assigned to the
equator. So something 120,000 meters south
of the equator is depicted as 9,880,000
47
Word in UTM Grids
UTM zones – lower 48 states
Universal Transverse Mercator (UTM Zone 14)
Elevation (feet)
UTM (Location of N at Memorial football stadium)
(Xo, Yo) = (693,500 4,521,383) Meters
Base map (streets) is from ArcGIS.com
ArcGIS.Com ready to use maps
including elevation services
http://www.arcgis.com/features/maps/earth.html
Elevation
Land Cover
Soils
Elevation Services
http://elevation.arcgis.com/arcgis/services
Summary Concepts
• The elevation surface represented by a grid digital
elevation model is used to derive slope important
for surface flow
• The eight direction pour point model approximates
the surface flow using eight discrete grid
directions.
• The D vector surface flow model approximates
the surface flow as a flow vector from each grid cell
apportioned between down slope grid cells.