VIEWS: 7 PAGES: 19 POSTED ON: 12/4/2011
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.
Pages to are hidden for
"EES 303"Please download to view full document