Automated Extraction of Coastal Dune High and Dune Low from High Resolution Lidar Digital Elevation Models Using GIS

By Kristy K. Guy 

Open File Report 2005-1344

U.S. Department of the Interior
U.S. Geological Survey



Introduction 

The impacts of extreme storms and hurricanes on the sandy coastlines of the Southeastern and Gulf coasts of the United States are being studied with the aid of high resolution LIDAR (LIght Detection And Ranging) data.  These data are collected along long stretches of coastline before and after major storms.  Included in the data collected is surface elevation.  These high density point elevations are processed into continuous raster surface digital elevation models (DEMs) that reveal geomorphic features such as beaches, dunes, and berms, and structures such as roads and buildings (Figure 1).  With further processing, additional information can be extracted from the DEMs.  Two important variables that are required to investigate, and perhaps one day help predict, storm impacts on coastal areas are 1) the location and elevation of the high point of the coastal dunes, referred to as dune high (Dhi), and 2) the location and elevation of the dune bases, or the "foot" of the dunes, referred to as dune low (Dlo).  

A semi-automated method was developed to extract  Dhi and Dlo from the lidar DEMs.  This method entailed using the DEMs, along with slope and aspect grids derived from the DEMs, to manually digitized lines relatively close to the actual Dhi and Dlo.  Two major drawbacks to this method are that 1) the large volume of data collected makes even this semi-automated process time consuming and tedious, and 2) the reliance on user judgment in difficult-to-interpret areas often produced results that differed from one user to another and from the same user from one time to another. 

An automated method of extracting Dhi and Dlo from lidar DEMs that largely alleviates these problems has been developed for use on the sandy Southeast and Gulf  coasts of the United States.  The method has been written into an Arc AML script that runs from a command line in ArcInfo Workstation, a popular GIS software product.  The output are GIS ready Dhi and Dlo point shapefiles that include several attributes that can assist in post-processing editing as well as elevation.  This report explains the auto_dhidlo.aml script, contains the script itself, and serves as its users' guide.

FIGURE 1 
Figure 1._ A) A color-scaled DEM of a coastal area showing an undeveloped area, roads, houses, high-rise buildings, dunes, and the beach and ocean.  The color scale indicates elevation in meters.  B) An aerial photograph of the same area (NOAA, July, 2005).

The Auto_dhidlo Script

The sandy beaches of ocean coastlines are often backed by one or more dunes or dune ridges.  Dunes defend the land behind them from the high water and waves generated by hurricanes and other extreme storms.  Researchers are interested in how the dunes are impacted by these large storms.  One way to measure the impact is to identify and compare the location and elevation of  Dhi and Dlo pre- and post- storm.  Dhi is the peak of a dune measured perpendicular to the shoreline or shore-normal.  Two Dhis are of interest to researchers, the Dhi closest to the shoreline  and the Dhi of the highest dune.    Dlo is the point where the greatest upward change in slope occurs seaward of the face of the most seaward dune.  It defines the first point at which the dunes would be impacted by high waves and water.(1)

Auto_dhidlo.aml is an ArcInfo AML script that runs from an ArcInfo Workstation command line.  ArcInfo with the GRID extension, or ArcGIS with the Spatial Analyst extension, are required to run the script.  The auto_dhidlo.aml file must be accompanied by eight kernel text files.  These nine files must either reside in the user's current workspace or the ArcInfo session must be set to look for them in another location (&AMLPATH command).  The script will run on either PC or UNIX platforms.  

Auto_dhidlo.aml attempts to find the maximum Dhi (mDhi), the most seaward Dhi (cDhi), and the Dlo, for every meter along the shoreline.  In effect, it uses each row or column, depending on the shoreline orientation, of an input DEM as an elevation transect.   The script is designed to use DEMs that have a resolution of 1 meter by 1 meter and have elevation values in centimeters.  It  will not find Dhi or Dlo points landward of buildings and trees.  It will not find Dhis at elevations higher or lower than a user specified limits.  A Dlo will not be found if a Dhi was not found for a particular row or column.  No Dhi or Dlo points will be located for areas without dunes, such as low overwash areas.

Input to the script includes an ArcInfo GRID (the DEM), an output name, and several user defined parameters described below.  Before running the script, users should take a few minutes to examine the input GRID to determine appropriate values for the minimum allowable Dhi, the maximum allowable Dhi, and the maximum Dhi search distance from the shoreline.   Defaults for these parameters will be set if the user enters zeroes for them, however, the program will likely run faster and produce better results if the user enters values based on the specific location.

The output consists of three Arc shapefiles, a "meta" text file, and a "watch" text file.  There are two shapefiles for Dhi, one for the mDhi and another for cDhi .  There is one shapefile for dune low (Figure 2).  The output shapefile names will end in _mh.shp, _ch.shp, and _l.shp for mDhi, cDhi, and Dlo respectively.  The "meta" file only provides information about the input, output, and parameters for that run.  It is not a standard Federal Geographic Data Committee (FGDC) metadata file (visit www.fgdc.gov for more information).  The "watch" text file can be referred to if problems occur during the script run.  All files are written to the current workspace.  The output files will likely contain some incorrect Dhi and Dlo points.  For example, some buildings may not have been successfully removed allowing Dhis to be located on a roofs.  Or it  the end of a pier or boardwalk may have been identified as a dune and a Dhi is located there. Therefore, the user should expect to edit the output for quality control.  If the output files contain many errors, the user may need to refine the values used for the various parameters and run the script again or consider an alternative method of finding Dhi and Dlo. 

FIGURE 2 
Figure 2.  An example of mDhi (purple dots), cDhi (red dots) and Dlo (blue dots) output points are shown on a color-scaled DEM.  The elevation along the red line is shown in profile in Figure 3.

FIGURE 3 
Figure 3.  The elevation profile of the red line drawn in Figure 2 with the mDhi (purple), cDhi (red), and Dlo (blue) output locations indicated.

Command Line Input

The command line usage and examples follow.  Required parameters are in angle brackets and optional parameters are in curly brackets.  All parameters must be entered in the correct order.  Parameter explanations, including instructions on how to skip parameters, are found below. 
 
Usage:
&run auto_dhidlo.aml <input DEM> <output name> <min Dhi> <max Dhi> <shoreline orientation> <ocean direction> <maximum Dhi search distance> {restraint file} {Dhi training file} {Dhi training item} {rotation angle} {dhi only}

Example with all parameters used:
&run auto_dhidlo.aml d:\fl\dems\stpete_dem stpete2005 250 750 ns west 125 fl_restraint fl_train elev_cm -45 dhi
Note:  Type command on one line. 
	
Example with minimum parameters used (see Figure 4):
&run auto_dhidlo.aml stpete_dem stpete2005 0 0 ns west 0

FIGURE 4
Figure 4.  Example of ArcInfo Workstation command line to run auto_dhidlo.aml.

Required Parameters

<input DEM>
The input DEM is an ArcInfo GRID created from high resolution (approximately 1-meter horizontal) elevation data such as that obtained by airborne lidar.  The GRID cell size should be 1-meter x 1-meter and the vertical units (elevation) in centimeters.  All ocean areas and the background should be zero or "nodata".  Enter the DEM's full pathname if necessary.  The GRID should be less than 49 kilometers long in the shoreline direction.  If not, refer to the sub_size hard-coded variable below.

<output name>
The output name will be used as the prefix for the output shapefiles.  Names longer than 29 characters will be truncated.  Do not enter a path.  The output files will be written to the current directory.  Existing output shapefiles with this prefix will be overwritten automatically.

<min Dhi>
The minimum allowable Dhi elevation should be set according to local conditions.  Enter the value in centimeters.  A default value of 250 will be used if the user enters 0 (zero).  If no Dhi greater than this value is found for a particular row (or column), no mDhi , cDhi or Dlo output points will be generated for that row (or column).

<max Dhi>
The maximum allowable Dhi elevation should be set according to local conditions.  Enter the value in centimeters.  This value limits the Dhi output but also helps in removing buildings and trees.  The user can opt to enter 0 (zero).  One of two things will happen if zero is used; the default value of 1000 cm will be used, or if the user also enters a train file and a train item (see below under Optional Parameters), the program will calculate a maximum Dhi based on the train file.  If the user enters a value here and also enters a train file and train item, the value entered here will be ignored.

<shoreline orientation>
The user must enter ns (for north-south) or ew (for east-west) for the shoreline orientation (Figure 5).  Use whichever orientation more closely approximates the actual shoreline orientation.  If the actual shoreline dramatically changes orientation within the GRID, such as around a cape, the user should split the GRID and run each separately.  Alternatively, the user can run the GRID once with a ns orientation and a second time with an ew orientation and edit the two outputs together.  The user has the option of rotating the GRID so that shoreline is oriented more closely to ns or ew.  However, rotating the GRID is not recommended.   See rot_angle under Optional Parameters below.

<ocean direction>
This is the direction of the ocean in relation to the land area in the GRID (Figure 5).  The ocean direction must be north or south if an ew shoreline orientation is used, or east or west if an ns orientation is used.    

<maximum Dhi search distance>
The maximum Dhi search distance limits the distance inland the script will search for Dhi.  Use the maximum distance, in meters (pixels), measured from the shoreline inland to just beyond the dunes in due north, south, east or west directions, NOT in a shore normal directions (Figure 5).  If the user enters 0 (zero) the default value of 150 will be used.  The maximum search distance allowed is 250.

FIGURE 5 
Figure 5.  The required shoreline orientation, ocean direction, and maximum Dhi search distance parameters are illustrated in their required cardinal directions.   Optional restraint lines, drawn here seaward of a seawall and a road, are also illustrated.

Optional Parameters

{restraint file}
The restraint file is an optional user generated ArcInfo line coverage or shapefile used to limit the search distance for Dhi.  The program will not find Dhi or Dlo landward of this line (Figure 5).  A restraint line can be useful when the maximum search distance parameter cannot be decreased but there are problem areas where the script is finding Dhi too far inland.  The lines do not have to be continuous nor encompass the entire extent of the GRID.  Use the word "none" (lower case, no quotes) if you do not want to enter a restraint file but want to enter any of the following parameters.  Leave blank if neither this nor any of the following parameters will be used.

{Dhi training file}
The Dhi training file is a point coverage or shapefile containing Dhi for the entire input GRID area.  It can cover a larger area but will only use data within the input GRID area.  The script  uses it to calculate appropriate maximum allowable Dhi elevations.  The Dhi file can be output from a previous Dhi analysis, from a different lidar date, and/or from the semi-automated method of extracting Dhi.  It must be in the same projection as the input GRID and the elevations must be in centimeters.  If a training file is entered, any value entered for maximum allowable Dhi elevation parameter will be ignored.  This option can save the user the time it would take to examine the GRID to determine an appropriate value.  Use the word "none" (lower case, no quotes) if you do not want to enter a training file but want to enter any of the following parameters.  Leave blank if neither this nor any of the following parameters will be used.

{Dhi training item} 
A Dhi training item must be specified if a training file is entered.  This is the item or field containing the Dhi elevation data in the training file.  Use the word "none" (lower case, no quotes) if you did not enter a training file but want to enter any of the following parameters.  Leave blank if neither this nor any of the following parameters will be used.
Example:
&run auto_dhidlo d:\fl\dems\stpete_dem stpete2005 250 750 ns west 125 fl_restrain none none -45 dhi 

{rotation angle}
The script looks for Dhi and Dlo from due north, south, east or west directions.  Generally the script will give satisfactory results even when the angle of the shoreline relative to due north-south or east-west is quite large.  However, if the results are unsatisfactory, the user can enter a rotation angle (between -45 (counterclockwise) and +45 degrees (clockwise)) that will more closely orient the shoreline to due north-south or due east-west.  It is strongly recommended that you DO NOT rotate the grid unless necessary.  Due to the nature of raster data, rotating GRIDs changes the values of their cells and can shift their location once re-rotated.  Very often, the resulting point elevations will not match the original input GRID pixel at the same location.  It will, however, match the value of a pixel close by.  At a minimum, the user should double the horizontal resolution of the results (2-meters instead of 1-meter).  

Rotating the GRID will also require more post-processing editing (Figure 6).  If the input GRID is more than 2 kilometers long in the shoreline direction, the script subsets the GRID to improve efficiency.    The subsets occur at 1000 meter intervals (see the hard-coded variable sub_size) beginning from the west if the GRID was rotated to an east-west orientation and beginning from the south if it was rotated to a north-south orientation.  The output file is an assemblage of the output from each subset.  The output data that fell at the ends of each subset, particularly the eastern or northern edges, will need to be edited since the rotated subset ends didn't contain the entire beach/dune width.  

If a rotation angle is entered, the shoreline orientation and ocean direction for the rotated GRID must be entered.  Enter 0 (zero) if you do not want to rotate the GRID but do want to specify the following parameter.  Leave blank if neither this nor the following parameter will be used.

FIGURE 6
Figure 6.  One of the effects of rotating the GRID is illustrated here.  (A) is an unrotated  GRID subset with a northeast-southwest oriented shoreline.  (B) is the same subset rotated during the analysis showing that the ends of the subset do not contain the full width of the data.  (C) shows output for the subset and the areas at the ends of the subset  where the output must be verified and/or edited.  

{dhi only}
The user has the option of finding Dhi only.  By entering "dhi" (lower case, no quotes) here, the script will find both mDhi and cDhi but not Dlo.  This will reduce processing time considerably.  There is no option for finding Dlo only since Dlo is based on Dhi.  Leave blank if you do not want to specify Dhi only.

Hard-Coded Variables

There are several hard-coded variables that the user may want to change under special circumstances.  These variables are set near the beginning of the AML.  They are as follows:

sub_size (default: 1000)
This is the subset length in the shoreline orientation direction in meters (pixels).  The script subsets the input GRID into 1000 meter sections to improve efficiency. The maximum number of subsets allowed is 49.  If the input GRID is more than 49 kilometers long in the shoreline orientation direction, increase the subset length so fewer than 49 subsets will be made.  Note that smaller subsets increase efficiency.  If a training file is used, the maximum allowable Dhi will be calculated for each subset.  Therefore, the output files for GRIDs that were rotated will contain data derived from varying maximum allowable Dhi parameters.  

dlo_lo (default: 100)
The minimum allowable Dlo in centimeters.  This can be changed to better accommodate local conditions.

edge_num (default: 200)
This number is used to help remove houses and trees from the input GRID.  If the script is finding Dhi on too many buildings and trees, decrease this number (try increments of 25).  If Dhi is not being found on steep dunes or scarps, lower this number.

dlo_slope (default: -10)
This number represents a steepness threshold used to help identify dunes.  Lower this number if  steeper dunes only are to be identified (try increments of 3).  Raise this number if flatter  dunes are to be identified.

sl2slsl (default: 7)
This is the maximum distance seaward in meters (pixels) from the steepest part of the dune that Dlo can be found.  Increase this number to allow a longer search distance, decrease to shorten the allowable search distance.  
	
train_max_margin (default: 75)
This is the number of vertical units (centimeters) that is added to the maximum allowable dune high the script calculates from the training file.  It is used only if a training Dhi file is provided.  Increase the number if the input GRID may have dunes higher than 75 cm above what was found in the training Dhi file (for example, where artificial dunes have been built since the training Dhi file was created).  Decrease this number to reduce the allowable elevation of Dhi.

neighbor_num (default: 7)
This is the size of the matrix used to find the number of neighbors an output point has (see "neighbors" in the Attribute sections of this report that follow).  The default is 7 x 7 pixels, 3 pixels in each direction from the point (Figure 8).  Increase the number to include a larger area, decrease it to reduce the search area.

ATTRIBUTES FOR MAXIMUM DHI (*_mh.shp):
x_coord	
UTM easting coordinate.

y_coord	
UTM northing coordinate.

mdhi_elev	
Maximum dune high (mDhi) elevation in centimeters.

FIGURE 7
Figure 7.  The dist2dlo and dist2shr attributes are illustrated on a gray-scale DEM with Dhi and Dlo outputs.  These attributes are measured in a due north, south, east or west direction, not shore normal.  These numbers can assist in post-processing editing.  

dist2shr	
The due north, south, east or west distance from mDhi to the shoreline in meters (pixels).  Do not use for shore normal distances (Figure 7).

Use in editing:  The user may find that the maximum Dhi search distance parameter was set too high and many mDhis were found too far inland.  If so, the user can mass edit out points by selecting them based on this attribute.  This might be a reasonable alternative to rerunning the script.  Similarly, if the minimum allowable Dhi elevation was set too low or many Dhis were found on the beach, mass editing out of points below a certain distance may be useful.

dist2dlo	
The mDhi to Dlo distance measured in a due north, south, east or west direction, in meters (pixels).  If dhi only was specified this value will be zero.

Use in editing:  The user may want to check points with high dist2dlo values.

neighbors	
The total number of mDhis in a 7x7 meter (pixel) matrix surrounding the point (Figure 8).

Use in editing:  The minimum number will be one, which means the point has no neighbors (only itself) within the area.  A value of one or other low numbers may indicate that the point is an outlier and a candidate for deletion.  Seven is the maximum number of neighbors possible (itself and 6 neighbors).  A value of seven or other high numbers can indicate more confidence in the point since other mDhis were found close by.  However, points with no or few neighbors can be correct.

FIGURE 8
Figure 8.  Close-up view of a gray-scale DEM showing individual pixels, output points (red and blue), and the 7 x 7 pixel matrix used to find the number of neighbors for the blue points.

ATTRIBUTES FOR MOST SEAWARD DHI (*_ch.shp):

x_coord	
UTM easting coordinate.

y_coord	
UTM northing coordinate.

cdhi_elev	
Most seaward dune high (cDhi) elevation in centimeters.

dist2shr	
The due north, south, east or west distance from cDhi to the shoreline in meters (pixels).  Do not use for shore normal distances (Figure 7).

Use in editing:  Same as for mDhi.

dist2dlo	
The cDhi to dune low distance in meters (pixels).  If dhi only was specified this value will be zero.

Use in editing:  Same as for mDhi.

neighbors	
The total number of  cDhis in a 7x7 meter (pixel) matrix surrounding the point (Figure 8).

Use in editing:  Same as for mDhi.

ATTRIBUTES FOR DLO (*_l.shp):

x_coord	
UTM easting coordinate.

y_coord	
UTM northing coordinate.

dlo_elev	
Dune low (Dlo) elevation in centimeters.

dist2shr	
The due north, south, east or west Dlo to the shoreline distance in meters (pixels).  Do not use for shore normal distances (Figure 7).  If Dhis were edited out using this attribute, editing out Dlos may be necessary.  This attribute may be useful for this purpose.

neighbors	
The total number of Dlos in a 7x7 meter (pixel) matrix surrounding the point (Figure 8).

Use in editing:  Same as for mDhi.  Caution:  Dlos tend to be more scattered than Dhis.  Therefore, points with no or few neighbors can be correct.

dunesteep	
This is a dune steepness indicator.  The larger the number the steeper the fore dune.

Use in editing:  The user may want to retain Dlos for the steeper dunes only.  The attribute can be used to mass delete Dlos for flatter dunes.

steep2dlo	
The distance between the steepest part of the fore dune and Dlo in meters (pixels).

Use in editing:  The user may want to retain only those Dlos that are close to the steepest part of the dune.  Use this attribute to mass delete Dlos greater than a certain distance.  The default maximum search distance for Dlo from the steepest part of the dune is 7 meters (pixels).  The user may want to eliminate Dlos with a steep2dlo of 7 since this can be considered an arbitrary cutoff.

Contact Information

Kristy K. Guy
Email: kguy@usgs.gov
Telephone: 727-803-8747 x 3063

Learn more on the Web: http://coastal.er.usgs.gov/hurricanes

References

Sallenger, A.H., Jr., 2000, Storm Impact Scale for Barrier Islands, Journal of Coastal Research, v. 16, no. 3, p. 890-895.

Disclaimer

This publication was prepared by an agency of the United States Government.  Neither the United States Government nor any agency thereof nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed in this report or represents that its use would not infringe privately owned rights. Reference therein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. 

No warranty, expressed or implied, is made by the USGS as to the accuracy of the data and related materials and (or) the functioning of the software. The act of distribution shall not constitute any such warranty, and no responsibility is assumed by the USGS in the use of these data, software, or related materials. 

Information contained within this electronic publication is considered public information and may be distributed or copied, and may not be copyrighted. We request that Public re-use of these materials bear a USGS byline/photo/image credit line.