Docstoc

EES 303

Document Sample
EES 303 Powered By Docstoc
					                                 Short Course on DEM and Morphotectonics
                    Using DEMs for Quantitative Geomorphic Analyses
INTRODUCTION

The overall goal of this skill set is to learn how to acquire, manipulate, and extract topographic
metrics from DEMs in an effort to understand the active tectonics of a region. In our case, the
region of interest is the northern Italian Apennines. A DEM for this region with a horizontal
resolution of 250 m and nominal vertical resoultion of 1 m has been prepared for your use. The
background information on the tectonic setting of the northern Apennines and all of Italy is the
subject of your readings and independent research through the semester. You will use this DEM
to collect data pertaining to the types and rates of tectonic deformation affecting the northeast vs.
southwest flanks of the northern Apennines.

Availability of DEMs is highly variable depending on the region of interest as well as the desired
resolution. The Italian 250 m DEM will work fine for the kinds of basic, region analyses we will
carry out. But higher resolution models are available and afford much more detailed analyses of
smaller areas. In an effort ot expose you to the acquisition of higher resolution data, the first part
of this skill set will have you find and download a 7.5 minute quadrangle from the United States.
However, you will learn how to make the actual measurements on the Italian DEM.

WHAT IS A DEM?

A DEM is a digital elevation model....with an emphasis on the word model. It is a rasterized
representation of the Earth’s surface. The raster pixels are almost always squares and the size of
the pixel determines the resolution of the final DEM. DEMs are made at a variety of scales by
most developed nations. In all places except the US, DEMS are expensive to purchase. In the
US, there are seeral forms of digital elevation data that are freely available over the web. In
addition, the US is in the process of processing radar data from a recent shuttle mission dedicated
to producing a DEM of the entire Earth between N60 and S60 degrees at the 90 m resolution.
That data set will be freely available next year. In the meantime, here are the common datasets
easily available over the web from the USGS or EROS data center:

NAME               Resolution        Format            Coverage                   Source
GTOPO301           30 seconds        2 byte integers   world    edcdaac.usgs.gov/gtopo30/gtopo30.html
                   or ~ 1 km

1 X 2 degrees2     3arcseconds       USGS DEM          US       edcwww.cr.usgs.gov/doc/edchome/ndcdb/ndcdb.html
(1:250 K)          or ~ 90 m

1 X 1 degrees      3arcseconds       SRTM – signed     World    ftp://e0srp01u.ecs.nasa.gov/srtm/version2/
                   or ~ 90 m         image data                 http://seamless.usgs.gov/ (seamless server)
                                                                http://srtm.csi.cgiar.org/

7.5 minute quad2   30 m              SDTS              US       www.mapmart.com (and others like GIS data depot)
(1:24 K)

7.5 minute quad    10-30 m           ArcGIS            US       http://seamless.usgs.gov/ (seamless server)
(1:24 K)
1 = these data come in a binary file that may need to be processed with a FORTRAN program first. See FJP for more
information.
2 = these data come in a format that Arc/Info can use directly. We will be dealing with SDTS format and will use the
SDTSIMPORT command in Arc. The sister command for USGS DEM format is DEMLATTICE.


WHAT METRICS WILL WE MEASURE?

The specific metrics are highlighted in BOLD throughout the procedures section. In some cases,
you are able to pull numbers off the DEM easily, but it will take a further calculation by hand or
in a spreadsheet to finalize the measurement. All of these metrics are outlined in Keller and
Pinter and I refer you to your text as a reference for the background on how to calculate specific
metrics like hypsometry or residual topography.

PROCEDURES

Anything in italics denotes a non-specific, or non-unique name that you must enter.

1. Log onto your account on a computer. Mount any disks that provide access to data or files.

        Log on and start ArcInfo and/or ArcMap.

2. For DEM data obtained from the seamless server, you can go directly to ArcGIS to display
and work with it. For SRTM or GTOPO30 data obtained from binary data blocks, download the
data and complete the following to unzip it and make it usable for ArcGIS:

        Open up a UNIX Cygwin shell and change directory to your working directory.

         >cd /filesystem/lastname

       Unzip the file.

        >gunzip *.gz

       If necessary, untar the files one by one, building them into an Arc/Info lattice one by one.
        This works best when you have one UNIX window open for working in Arc/Info, and
        another window open for working in UNIX.

        >tar xvf filename2.tar

        notice how several files are extracted from the tar.

        For the PC, unzipping and untaring files may be accomplished by WINZIP or the
         decompression routines in Windows XP.

For using the SRTM binary data specifically on the PC, follow these instructions:
Users of ARC/INFO or ArcView can display the DEM data directly after renaming the file
extension from .HGT to .BIL. However, if a user needs access to the actual elevation values for
analysis in ARC/INFO the DEM must be converted to an ARC/INFO grid with the command
IMAGEGRID. For IMAGEGRID to work there must be a separate header file whose name
(including case) is exactly the same as the image file name. The header file must also carry the
extension .hdr . The contents of a sample file that works with 3 arc-second SRTM file
N37W105.hgt are below. Note that the coordinates of the binary data block, listed in the file
name indicate the latitude and longitude of the lower left hand corner of the data block.
--------------------------------------------------------------------------
BYTEORDER                  M
LAYOUT                 BIL
NROWS                    1201
NCOLS                    1201
NBANDS                   1
NBITS                    16
BANDROWBYTES                          2402
TOTALROWBYTES                         2402
BANDGAPBYTES                          0
NODATA                   -32768
ULXMAP                   -105.0
ULYMAP                   38.0
XDIM                     0.000833333333333
YDIM                     0.000833333333333
--------------------------------------------------------------------------

For a 30 m DEM (1 degree by 1 degree) where the original data are named N40W076.hgt
--------------------------------------------------------------------------
BYTEORDER                  M
LAYOUT                 BIL
NROWS                    3601
NCOLS                    3601
NBANDS                   1
NBITS                    16
BANDROWBYTES                          7202
TOTALROWBYTES                         7202
BANDGAPBYTES                          0
NODATA                   -32768
ULXMAP                   -76.0
ULYMAP                   41.0
XDIM                     0.0002777777777778
YDIM                     0.0002777777777778
------------------------------------------------------------------------

To correct for no data and negative elevations:

outgrid = con (ingrid < 8888, ingrid)
To fill data gaps, you must average good elevation values, and use them to estimate the
elevations in the holes with no data. To do this, use the command line in Arc, or the rater
calculator in ArcGIS to run the following command

      Con (isnull(ingrid), focalmean (ingrid, circle, 2), Ingrid)

3. If you happen to have more than one lattice that needs to be merged into a new lattice use the
following commands in Arc. The equivalent steps in ArcGIS are the Mosaic and Mosaic to new
raster tools:

       Arc: latticemerge merged_grd
       Enter lattices to be merged(Type END or a blank line when
       done):
       =========================================================
       Enter the 1st lattice: name1
       Enter the 2nd lattice: name2
       Enter the 3rd lattice: name3
       Enter the 4th lattice: name4
       Enter the 5th lattice: end

      Now, display the assembled map.

       Arc: grid
       Grid: mape merged_grd
       Grid: image merged_grd

      Fix any problems with how the quadrangles merged.

       Grid: fixed_grd =
       con(isnull(merged_grd),focalmean(merged_grd,
       circle,3),merged_grd)

      Now display the fixed grid

       Grid: mape fixed_grd
       Grid: image fixed_grd

      Make a hillshaded portrayal of the assembled DEM

       Grid: fixed_hs = hillshade (fixed_grd, #, #, #, 2)
       Grid: image fixed_hs

4. Coordinate systems, Projections, and datums (from: http://webhelp.esri.com/arcgisdesktop/
   9.2/index.cfm?TopicName=An_overview_of_map_projections)
Coordinate systems enable geographic datasets to use common locations for integration. A
coordinate system is a reference system used to represent the locations of geographic features,
imagery, and observations such as GPS locations within a common geographic framework. Each
coordinate system is defined by:

       Its measurement framework which is either geographic (in which spherical coordinates
        are measured from the earth's center) or planimetric (in which the earth's coordinates are
        projected onto a two-dimensional planar surface).


       Unit of measurement (typically feet or meters for projected coordinate systems or
        decimal degrees for latitude-longitude).


       The definition of the map projection for projected coordinate systems.


       Other measurement system properties such as a spheroid of reference, a datum, and
        projection parameters like one or more standard parallels, a central meridian, and possible
        shifts in the x- and y-directions.

There are two common types of coordinate systems used in GIS:

       A global or spherical coordinate system such as latitude-longitude. These are often
        referred to as geographic coordinate systems.


       A projected coordinate system based on a map projection such as transverse Mercator,
        Albers equal area, or Robinson, all of which (along with numerous other map projection
        models) provide various mechanisms to project maps of the earth's spherical surface onto
        a two-dimensional Cartesian coordinate plane. Projected coordinate systems are
        sometimes referred to as map projections.

Coordinate systems (either geographic or projected) provide a framework for defining real-
world locations. In ArcGIS, the coordinate system is used as the method to automatically
integrate the geographic locations from different datasets into a common coordinate framework
for display and analysis.

A spatial reference in ArcGIS is a series of parameters that define the coordinate system and
other spatial properties for each dataset in the geodatabase. It is typical that all datasets for the
same area (and in the same geodatabase) use a common spatial reference definition.

An ArcGIS spatial reference includes settings for:

       The coordinate system
      The coordinate precision with which coordinates are stored (often referred to as the
       "coordinate resolution")
      Processing tolerances (such as the cluster tolerance)
      The spatial or map extent covered by the dataset (often referred to as the "spatial
       domain")

A geographic coordinate system (GCS) uses a three-dimensional spherical surface to define
locations on the earth. A GCS is often incorrectly called a datum, but a datum is only one part of
a GCS. A GCS includes an angular unit of measure, a prime meridian, and a datum (based on a
spheroid).

A point is referenced by its longitude and latitude values. Longitude and latitude are angles
measured from the earth's center to a point on the earth's surface. The angles often are measured
in degrees (or in grads). The following illustration shows the world as a globe with longitude and
latitude values.




In the spherical system, horizontal lines, or east–west lines, are lines of equal latitude, or
parallels. Vertical lines, or north–south lines, are lines of equal longitude, or meridians. These
lines encompass the globe and form a gridded network called a graticule.

The line of latitude midway between the poles is called the equator. It defines the line of zero
latitude. The line of zero longitude is called the prime meridian. For most geographic coordinate
systems, the prime meridian is the longitude that passes through Greenwich, England. Other
countries use longitude lines that pass through Bern, Bogota, and Paris as prime meridians. The
origin of the graticule (0,0) is defined by where the equator and prime meridian intersect. The
globe is then divided into four geographical quadrants that are based on compass bearings from
the origin. North and south are above and below the equator, and west and east are to the left and
right of the prime meridian.

Latitude and longitude values are traditionally measured either in decimal degrees or in degrees,
minutes, and seconds (DMS). Latitude values are measured relative to the equator and range
from -90° at the South Pole to +90° at the North Pole. Longitude values are measured relative to
the prime meridian. They range from -180° when traveling west to 180° when traveling east. If
the prime meridian is at Greenwich, then Australia, which is south of the equator and east of
Greenwich, has positive longitude values and negative latitude values.

It may be helpful to equate longitude values with X and latitude values with Y. Data defined on a
geographic coordinate system is displayed as if a degree is a linear unit of measure. This method
is basically the same as the Plate Carrée projection.

Although longitude and latitude can locate exact positions on the surface of the globe, they are
not uniform units of measure. Only along the equator does the distance represented by one
degree of longitude approximate the distance represented by one degree of latitude. This is
because the equator is the only parallel as large as a meridian. (Circles with the same radius as
the spherical earth are called great circles. The equator and all meridians are great circles.)

Above and below the equator, the circles defining the parallels of latitude get gradually smaller
until they become a single point at the North and South Poles where the meridians converge. As
the meridians converge toward the poles, the distance represented by one degree of longitude
decreases to zero. On the Clarke 1866 spheroid, one degree of longitude at the equator
equals 111.321 km, while at 60° latitude it is only 55.802 km. The standard conversion of
degrees longitude to km as a function of latitude is:

                                     km = 111.321 * cos (lat)

Because degrees of latitude and longitude don't have a standard length, you can’t measure
distances or areas accurately or display the data easily on a flat map or computer screen.

A projected coordinate system is defined on a flat, two-dimensional surface. Unlike a
geographic coordinate system, a projected coordinate system has constant lengths, angles, and
areas across the two dimensions. A projected coordinate system is always based on a geographic
coordinate system that is based on a sphere or spheroid.

In a projected coordinate system, locations are identified by x,y coordinates on a grid, with the
origin at the center of the grid. Each position has two values that reference it to that central
location. One specifies its horizontal position and the other its vertical position. The two values
are called the x-coordinate and y-coordinate. Using this notation, the coordinates at the origin are
x = 0 and y = 0.

On a gridded network of equally spaced horizontal and vertical lines, the horizontal line in the
center is called the x-axis and the central vertical line is called the y-axis. Units are consistent
and equally spaced across the full range of x and y. Horizontal lines above the origin and vertical
lines to the right of the origin have positive values; those below or to the left have negative
values. The four quadrants represent the four possible combinations of positive and negative X
and Y coordinates.

When working with data in a geographic coordinate system, it is sometimes useful to equate the
longitude values with the X axis and the latitude values with the Y axis.
Universal Transverse Mercator (UTM) system is a specialized application of the Transverse
Mercator projection. The globe is divided into 60 north and south zones, each spanning 6° of
longitude. Each zone has its own central meridian. Zones 1N and 1S start at -180° W. The limits
of each zone are 84° N and 80° S, with the division between north and south zones occurring at
the equator. The polar regions use the Universal Polar Stereographic coordinate system.

The origin for each zone is its central meridian and the equator. To eliminate negative
coordinates, the coordinate system alters the coordinate values at the origin. The value given to
the central meridian is the false easting, and the value assigned to the equator is the false
northing. A false easting of 500,000 meters is applied. A north zone has a false northing of zero,
while a south zone has a false northing of 10,000,000 meters.




While a spheroid approximates the shape of the earth, a datum defines the position of the
spheroid relative to the center of the earth. A datum provides a frame of reference for
measuring locations on the surface of the earth. It defines the origin and orientation of latitude
and longitude lines.
Whenever you change the datum, or more correctly, the geographic coordinate system, the
coordinate values of your data will change. Here are the coordinates in DMS of a control point in
Redlands, California, on the North American Datum of 1983 (NAD 1983 or NAD83).

-117 12 57.75961
34 01 43.77884

Here's the same point on the North American Datum of 1927 (NAD 1927 or NAD27).

-117 12 54.61539
34 01 43.72995

The longitude value differs by approximately three seconds, while the latitude value differs by
about 0.05 seconds.

NAD 1983 and the World Geodetic System of 1984 (WGS 1984) are identical for most
applications. Here are the coordinates for the same control point based upon WGS 1984.

-117 12 57.75961
34 01 43.778837
In the last 15 years, satellite data has provided geodesists with new measurements to define the
best earth-fitting spheroid, which relates coordinates to the earth's center of mass. An earth-
centered, or geocentric, datum uses the earth's center of mass as the origin. The most recently
developed and widely used datum is WGS 1984. It serves as the framework for locational
measurement worldwide.

A local datum aligns its spheroid to closely fit the earth's surface in a particular area. A point on
the surface of the spheroid is matched to a particular position on the surface of the earth. This
point is known as the origin point of the datum. The coordinates of the origin point are fixed, and
all other points are calculated from it.




The coordinate system origin of a local datum is not at the center of the earth. The center of the
spheroid of a local datum is offset from the earth's center. NAD 1927 and the European Datum
of 1950 (ED 1950) are local datums. NAD 1927 is designed to fit North America reasonably
well, while ED 1950 was created for use in Europe. Because a local datum aligns its spheroid so
closely to a particular area on the earth's surface, it's not suitable for use outside the area for
which it was designed.
5. Project your grid into a useful georeferenced projection such as UTM. These transformations
are built directly into the toolbox of ArcGIS, for ArcInfo, use the following commands in Grid:

       Grid: outgrid_utm = project(outgrid_grd, #, bilinear)

       Please provide the input and output map projections.
       Use INPUT to provide the input projection, OUTPUT
       to provide the OUTPUT projection, and END to finish.

       Project:     input
       Project:     projection geographic
       Project:     units dd
       Project:     parameters
       Project:     output
       Project:     projection utm
       Project:     units meters
       Project:     zone 32
       Project:     parameters
       Project:     end

6. Viewing a raster data set as a hillshade. In ArcGIS, this can be accomplished by accessing
the hillshade tool in the Spatial Analyst Extension. For ArcInfo:

       Grid:    mape outgrid_utm
       Grid:    image outgrid_utm
       Grid:    outgrid_hs = hillshade (outgrid_utm, #, #, #, 2)
       Grid:    image outgrid_hs

7. Measure SINUOSITY of a mountain front. With a hillshaded image showing on the screen.

       Grid: measure length *

     This allows you to interactively measure lengths on the screen. The output is in cell
      dimensions. Remember, your cells are 200 m on a side.

Using ArcGIS, measuring distances and areas is very easy using the measure tool (looks like a
ruler) on the main menu. You can interactively measure distances and areas on the screen with
this tool. Note, you can also define the units of measure.

8. Extract CROSS-SECTIONS from the DEM. Upload those cross-section to a PC or Mac and
   manipulate it in Excel to make the plot.

       Grid:    clear
       Grid:    mape outgrid_utm
       Grid:    image outgrid_utm
       Grid:    surface lattice outgrid_utm
       Grid: surfacedefaults
       Grid: surfacezscale 10 constant
       Grid: surfacedrape mesh fishnet

     This suite of commands will render your grid as a 3-D perspective. The following
      commands will show you how to interactively extract a single profile. Follow the
      instructions on the screen to draw a box where the graph will be drawn as well as the line
      of cross section.

       Grid: clear
       Grid: image outgrid_utm
       Grid: surfaceprofile * * xsection

     Now you have to download the cross section to your PC or Mac

       Grid: q
       Arc: info
       ENTER USER NAME> arc
       ENTER COMMAND > select XSECTION
       ENTER COMMAND > export /filesystem/XSECTION.dat sdf
       ENTER COMMAND > q stop

For two or more exactly parallel cross-section lines that you need to make swath profiles, you
will give the surfaceprofile command a file with the appropriate coordinates. Since you are
operating on a UTM projected grid, you must provide those coordinates as UTM meters. This is
sometime a bit awkward to do because we think easily in long-lat coordiantes. I suggest using
the cellvalue command in Grid to explore your drainage basin and decide which coordinates are
appropriate. You must use the Generate command to create a line coverage that will be used
with surfaceprofile.

       Grid: cellvalue outgrid_utm *
       Grid: arc generate outgrid_lines
       Generate: lines

       Enter lines. Terminate lines by entering END at X,Y prompt
       Terminate input by entering END at ID: prompt

       ID:
       X,Y:
       X,Y:
       X,Y: end
       .
       .
       .
       ID: end
       Generate: q
       Grid: arc build outgrid_lines line
       Grid: surfaceprofile * outgrid_lines xsections

      Export the data from INFO as you did in the previous step.

      Now you are ready to transfer the ascii file xsections.dat to your PC or Mac. Use a FTP
       program to do this. FTP programs reside on virtually every PC and Mac in the
       Department and you should know how to use them. Follow the demonstration in class.

Using ArcGIS and EZProfiler 9.x to extract cross-sections and pseudo-swath profiles v. 1

      Add a raster to your active layer

      Click on the “D” drawing icon. Draw your line on the raster, double-clicking to finish.
       The data get output to a point file. You can adjust the number of points under the “PS”
       settings icon.

      Select the point shapefile in the layers column, then click on the excel icon, your data are
       output to excel.

      Now…if you want multiple profiles, you need to create a new shapefile in arc catalogue,
       make sure it carries the same coordinate system as the raster. Next, go into editor and
       draw the lines, save the lines, stop the editor.

      You can generate perfectly parallel lines with the generate command, and knowing
       precisely where you want the lines to begin and end. Use the identify tool in the main
       toolbar to scan the raster and find the starting and ending points of your lines. Next, get
       those coordinates entered into a simple text file that looks like the one below. Note that
       the numbers “100”, “200”, and “300” are simply the identification number for three
       separate lines.

               100
               228000, 4653000
               249000, 4653000
               end
               200
               228000, 4652000
               249000, 4652000
               end
               300
               228000, 4651000
               249000, 4651000
               end
               end

      Use the generate command, with the lines option active. The output from this command
       is a coverage – similar to a shape file, but not exactly. You need to convert the coverage
       toa shapefile, which can be done several different ways, the easiest of which is to right
       click on the file in ArcCatalog and export it as a shapefile.

      You next need to assign the shapefile the coordinate system of your raster – you should
       know how to do this already – the easiest way is to right click on the file name in
       ArcMap or ArcCatalog and go to properties and import a pre-existing file coordinate
       system.

      Add the new shapefile to the active layers.

      Select the raster in the layers column, then click on the “M” icon, the data get output to a
       new point shapefile. You can adjust the number of points under the “PS” settings icon.

      Select the point shapefile in the layers column, then click on the excel icon, your data are
       output to excel. Note that you can adjust the profiles horizontally and sample the Z-
       values for the min, max, and mean. These are called swath profiles.


Using ArcGIS and EZProfiler 9.x to extract cross-sections and pseudo-swath profiles v. 2

      In ArcCatalog, create a new line shapefile and assign it a coordinate system that matches
       the projected DEM from which you will be extracting the data.

      Add the shapefile to the active layers (it will be blank)

      Begin the Editor and make sure it is editing the new line shapefile you just created

      Click on the pencil tool and then click and draw a line in any orientation. Double click to
       end the line.

      Look under the editor pull-down menu for the option “copy parallel lines”. Using this
       copy utility, you can create a perfectly parallel line on one or both sides of the line you
       drew, and you could repeat the steps as many times as you wish. You also have the
       option of specifying the distance between parallel lines.

      Once the lines are drawn, use the multiple lines extraction button in EZProfiler to extract
       the elevation values and follow steps outlined above to get them into Excel for plotting.

Use of EZProfiler 9.x in ArcMap to extract a long profile.

      Complete the ArcHydro steps of Lab 3 to result in a stream polyline shapefile.

      Begin Edit, select the tributaries and delete them, leaving only the trunk channel.

      Select all segments of the trunk channel and merge them into one polyline using the
       merge option under the Editor pull down menu toolbox.
      Save the edits and exit the editor.

      Make the merged line active in the layers column. Select it using the “select features”
       tool in the toolbar. This is he tool that has a white arrow next to a blue and white box. It
       is typically on the toolbar next to the black arrow. Note….the interesting thing about this
       selection tool is that if your lines are merged, then you can more or less interactively
       select the individual line segments that you might need or want to make the long profiles.

      Follow the steps above for extracting a cross-section, but use the “S” icon for extraction
       of a single line.

      Plot it up in Excel.

Using the ArcInfo macro to extract real swath profiles for a boxed area (Courtesy of Mike
Oskin).

Make sure the file swathprofile.aml is in your working directory.

Arc: &amlpath c:\<enter your path here>

Arc: &run swathprofile.aml <in_grd> <swath_name> <width in map units> <bin
size in map units>

Follow on screen, interactive instructions, noting to not make the swath name too long (not
greater than ~5-6 letters). The <width in map units> should be some number like 2000 for a
two-kilometer wide box in a UTM meter grid. Bin size in map units should be something like 10
to 100. Note that you need to interactively define the box of the swath profile. After the script
runs, instructions on how to get the data out of the ArcInfo table and into ASCII form for
plotting are presented. Note, there is a typo – it is zonalstats, not zonalstate. Basically, an ascii
table with a .txt extension will have been created automatically and you can simply use this to
plot up your swath. Column 3 is the X-distance in x-increments (not summed), column 4 is min
elevation, column 5 is maximum elevation, and column 6 is mean elevation. Get this file into
Excel or SigmaPlot and plot it up.


9. The next step is to EXTRACT YOUR DRAINAGE BASIN. I will assign one to you in
   class. You will have to find that drainage basin by comparing the master topographic map
   (the road map) to the hillshaded portrayal of your assembled dataset.

      The general procedure that you are following here is to fill all artifical holes in the dataset,
       then to define the direction that the water is flowing, which cells are streams and which
       are hillslopes, then define the boundaries of a watershed.

       Arc: Grid
       Grid: mape outgrid_utm
       Grid: image outgrid_utm
    Grid:    fill outgrid_utm outgrid_fill
    Grid:    outgrid_fd = flowdirection (outgrid_fill)
    Grid:    outgrid_acc = flowaccumulation (outgrid_fd)
    Grid:    mape outgrid_acc
    Grid:    image outgrid_acc

   Now it is important to find your stream on this map and zoom in on its mouth. Follow the
    instructions in class illustrating how to do this. You MUST click on a stream precisely,
    not a cell next to the stream.

    Grid: name_ws = watershed (outgrid_fd, selectpoint
    (outgrid_fill, *))

   Now you are ready to clip out your watershed using a cookie cutter shaped just like your
    watershed.

    Grid:    name_clip = gridpoly (name_ws)
    Grid:    arc latticeclip outgrid_fill name_clip name_grd
    Grid:    clear
    Grid:    mape name_grd
    Grid:    image name_grd

   Now you are ready to define the streams just for this basin. So repeat the flowdirection
    and flowaccumulation steps. The critical step here is the con command where you are
    setting the threshold for what will and will not be a stream. In the example below, I’ve
    chosen “100” as that threshold. This means that any cell with 100 units of water will be a
    stream, all other cells are hillslopes. In these instructions name is the name of your grid,
    typically “outgrid”.

    Grid:    name_fd = flowdirection (name_grd)
    Grid:    name_acc = flowaccumulation(name_fd)
    Grid:    mape name_acc
    Grid:    image name_acc
    Grid:    name_strm = con (name_acc > 100, 1)
    Grid:    clear
    Grid:    mape name_strm
    Grid:    image name_strm
    Grid:    name_cov = gridline (name_strm)

   Now you have made a line coverage of the streams. Display this on the screen with the
    drainage basin grid and save it as an eps file.

    Grid:    clear
    Grid:    mape name_grd
    Grid:    image name_grd
    Grid:    linecolor blue
    Grid:    arcs name_cov
       Grid: linecolor red
       Grid: arcs name_clip
       Grid: screensave eps name.ps

10. With the drainage basin and drainage lines clearly displayed, it should be very straight
    forward to investigate DRAINAGE BASIN ASYMMETRY. You will do this interactively
    on the screen. Choose the trunk channel. And then measure the areas on either side of the
    channel.

       Grid: measure area *

       The area is reported in square meters.

      Clean up the display and measure BASIN ELONGATION, defined as the maximum
       lenght, divided by the maximum width.

       Grid:    clear
       Grid:    mape name_grd
       Grid:    image name_grd
       Grid:    linecolor blue
       Grid:    arcs name_cov
       Grid:    linecolor red
       Grid:    arcs name_clip
       Grid:    measure length *
       Grid:    measure length *

11. The next step is to extract data that will be used to determine the drainage basin
    HYPSOMETRY.

       Grid:    clear
       Grid:    mape name_grd
       Grid:    image name_grd
       Grid:    name_slice = slice (name_grd, eqinterval, 10)

      To get the data out of the GIS, once again, you need to go into info and export it. Follow
       the steps in 6. above, but this time, select NAME_SLICE.VAT. And then FTP that
       extracted ascii data file to your PC or Mac for manipulation in Excel. The data will mport
       to Excel as 2 columns. In column A will be the numbers 1 through 10. In column B will
       be the number of cells that correspond to these equal interval slices. Here is what you
       need to do in Excel:

       (a) Fill column C with the values of column A divided by 10. Column C is now h/H and
           will be plotted on the Y-axis.
       (b) Fill Column D with the values of column B *302. This calculates the true area of the
           equal interval slice in square meters.
       (c) Fill column E with the relative percentage area for each equal interval slice. The
           relative percentage is equal to the value in column D divided by the sum of column
           D.
       (d) Fill column F with the cumulative percentage of column E. Column F should stretch
           from nearly 0 in row 1, to very nearly 1 in row 10. Column F is a/A and will be
           plotted on the X-axis.
       (e) Make a plot of column F vs. column C.

12. Lastly, you need to extract the basin’s LONGITUDINAL PROFILE. This is a bit
    involved and you are going to have to be patient with it because it requires a fair bit of
    cleaning up the coverages using ArcEdit. There are two different procedures that work –
    one involves using the SETMASK command in GRID; the other involves using surfaces
    and the SURFACEPROFILE command as instructed above in procedure (8).

Both procedures need a (1) UTM-projected grid of the cut-out filled drainage basin (basin_utm),
    (2) a grid of the flowaccumulation performed on this cut-out drainage basin (basin_acc), and
    a (3) line covereage of the stream network in this drainage basin (strm_cov). Make sure
    your grids are in UTM projection – see procedure 5 to do this. You should already be in
    UTM projection if you’ve been following this skill set from the beginning. Use a large
    threshold to creat (3) with the CON command – a number like 500 is fine. You need to do
    procedure 12b; 12a is provided for your reference.


12a.

       Grid: copy strm_cov strm2_cov
       Grid: q
       Arc: arcedit
       Arcedit: mape strm2_cov
       Arcedit: edit strm2_cov
       Arcedit: de all
       Arcedit: ef arc
       Arcedit: draw
       Arcedit: select many

      Arcedit now allows you to interactively select streams you do not want. As you do, they
       turn yellow. When you are ready to delete them, type the 9 key.

       Arcedit: delete

      Do this until all but the trunk stream remains. At that point, you want to transform it to a
       single line.

       Arcedit:      select all
       Arcedit:      unsplit
       Arcedit:      save
       Arcedit:      q
      Now you take the next step to use this line to clip out the elevations that correspond only
       with the trunk stream.

       Arc: grid
       Grid: mape basin_utm
       Grid: image basin_utm
       Grid: setcell basin_utm
       Grid: stream_grd = linegrid (strm2_cov)
       Grid: setmask stream_grd
       Grid: lp_grd = basin_utm
       Grid: setmask off
       Grid: lp_int = int(lp_grd)

      The long profile information resides in info in the lp_int.vat file. You must go into info,
       as before, to export it to an ascii file.

       Grid: q
       Arc: info
       ENTER USER NAME> arc
       ENTER COMMAND > select LP_INT.VAT
       ENTER COMMAND > export /filesystem/lp.dat sdf
       ENTER COMMAND > q stop

      When the data are exported, there will be two columns. The first are elevations, the
       second are the number of that particular elevation. The assumed distance between any
       adjacent elevation is 200 m (in our case, because the spacing is 200 meters). This of
       course is not precisely true for cells that are in contact diagonally...but it is close enough
       for our purposes. You need to bring the ascii file lp.dat to Excel via FTP. Then you must
       make a plot of distance vs. elevation. The only tricky thing that you must keep track of is
       if a given elevation has more than one value. In other words, if there is the elevation 10
       in column one and the corresponding number 1 in column 2, you can be assured that
       there is only one place on the whole long profile where there is the elevation 10. On the
       other hand, if you find a 10 in column 1 and a corresponding 3 in column 2, that means
       that there are three pixels of elevation 10 adjacent to one another. You need to take that
       into account in plotting up the long profile or it will underestimate the true length of the
       stream.

      Repeat these exact steps (follow the Arcedit commands) on the flowaccumulation grid:
       basin_acc.

12b.

      Follow all of the Arcedit steps of 12a to select a single line for the long profile, then..

       Grid: mape basin_utm
    Grid:    image basin_utm
    Grid:    linecolor red
    Grid:    arcs strm_cov
    Grid:    surface lattice basin_utm
    Grid:    surfacedefaults
    Grid:    surfaceprofile * strm_cov lp

   lp signifies “long profile”, follow the screen instructions to draw the box

    Grid:    clear
    Grid:    mape basin_acc
    Grid:    image basin_acc
    Grid:    arcs strm_cov
    Grid:    surface lattice basin_acc
    Grid:    surfacedefaults
    Grid:    surfaceprofile * strm_cov ap

   ap signifies “area profile”, follow the screen instructions to draw the box

    Grid:    info
    ENTER    USER NAME> arc
    ENTER    COMMAND > select LP
    ENTER    COMMAND > export c:/path/lp.dat sdf
    ENTER    COMMAND > export c:/path/ap.dat sdf
    ENTER    COMMAND > q stop
    Grid:    describe basin_utm
    Grid:    q

   Record the cellsize that is displayed on the screen

    Arc: q

   Now, take the ascii file created by the export command in INFO and process it through
    the FORTRAN program slpar.f This program imports the data, smooths it, and then
    arranges it in columns that can be easily plotted in Excel or SigmaPlot.

   Run the program by compiling the source code, then executing the .exe file.

   Take the output from the program and plot it up making graphs of the long profile, the
    log area-log slope and log area-log distance.

				
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
views:7
posted:12/4/2011
language:English
pages:19