A GIS add-in for automated measurement of sand dune migration using LiDAR-derived multitemporal and high-resolution digital elevation models

Understanding how sand dunes form and change has long been an important research topic in Earth and planetary surface processes (e.g., Wasson and Hyde, 1983; Rubin and Hesp, 2009; Bridges et al., 2012). Since the 1970s, studies on sand dunes have gradually shifted from single dunes to complex dune systems (Livingstone et al., 2007). In addition to many field studies, numerical and cellular automaton models as well as flume and landscape-scale experiments have been tested for better understanding of sand dunes (Alhajraf, 2004; Zhang et al., 2010; Taniguchi et al., 2012; Araújo et al., 2013; Ping et al., 2014).

Although remote sensing has been used for sand dune studies for decades (Hunter et al., 1983; Gay, 1999; Levin et al., 2004; Yao et al., 2007), the majority of these studies are based on visual interpretation and manual measurement methods, which are often labor-intensive, time-consuming, and error-prone, especially in large dune fields. Leprince et al. (2007) proposed an automated method named COSI-Corr (coregistration of optically sensed images and correlation) for ground deformation measurements, which has also been used for sand dune studies (Vermeesch and Drake, 2008; Necsoiu et al., 2009; Hermas et al., 2012). However, the relatively low accuracy associated with medium- and low-resolution images is a major limitation of the COSI-Corr method.

Light detection and ranging (LiDAR) has provided unprecedented data sets for sand dune studies (Mitasova et al., 2005; Saye et al., 2005; Reitz et al., 2010; Baitis et al., 2014; Ewing et al., 2015). However, visual analysis and manual measurement remain a common practice in these studies. Hardin et al. (2014) studied analysis methods for coastal LiDAR time-series data in a geographic information system (GIS). Dong (2015) proposed the pairs of source and target points (PSTP) method and built a model for automated measurement of dune migration using multitemporal LiDAR data. The work of Dong (2015) was implemented in ArcGIS ModelBuilder, which is not a generic GIS software. As a continuing study of Dong (2015), this research note focuses on the development of a generic ArcGIS add-in toolbar using the Python programming language, for automated measurement of dune migration directions and migration rates based on multitemporal and high-resolution digital elevation models (DEM), which can be derived from LiDAR data or other sources. The add-in is freely available (see ^{footnote 1}), and it can be a useful tool for geoscientists studying migration of desert or coastal sand dunes.

Bagnold (1941) observed that sand avalanching and slumping events caused by gravity occur in the inclination direction of the dune slip face. This observation is used as the theoretical foundation of the PSTP method. Although it is difficult or impossible to monitor every step of dune movement, multitemporal data sets such as high-resolution DEMs derived from LiDAR data can be used as snapshots for studying sand dune migration. In the PSTP method, the angle of repose (AOR) is used as a sensitive indicator of dune migration. Dune slip faces are extracted from DEMs based on AOR values, and the centerlines of slip faces are obtained by vectorization. The slip face centerlines extracted from the first DEM and the second DEM are called source lines and target lines, respectively (Fig. 1). For any random target point t_n on a target line, a search radius r is used to search for the nearest source point s_n on the source line. If no source point is found, the target point is deleted; if a source point is found, the direction of the source point relative to the target point is called source direction, which is statistically related to the prevailing wind direction (Dong, 2015). The source direction is further compared with the orientation of the line segment on the source line to validate the target point. If the source direction is not perpendicular to the line segment that contains the source point, the target point is deleted. For example, the source direction from s₁ to t₁ in Figure 1 is not perpendicular to the line segment that contains s₁, and so t₁ is deleted. For remaining target points, the migration distance s_nt_n is converted to migration rate (m/yr) by taking into account the time interval (number of days) between the two DEMs.

An ArcGIS add-in was implemented using the Python programming language. The add-in toolbar has six combo boxes and a command button (Fig. 2). To reduce the size of the toolbar, data acquisition dates are contained in the input DEM names in the form of YYYYMMDD.

The items on the toolbar in Figure 2 are described below:

• (1) DEM1: The first DEM raster, which can be created from LiDAR data or other data sources. The data acquisition date is contained in the DEM layer name in the format of YYYYMMDD, and the YYYYMMDD string can be anywhere in the DEM name as long as it is the first eight numbers, for example, A20090124DEM1. The DEM layer name can be changed by the user in ArcGIS, and it can be different from the actual file name.

• (2) DEM2: The second DEM raster (similar to DEM1).

• (3) AOR: Angle of repose for sand dune slip faces. AOR is usually around 34°, depending on the sand grain size, shape, and moisture content. Users can select/input a range, such as 30–35, for AOR values.

• (4) Min-Dist: The minimum distance between two random points. The unit of distance is the same as the linear unit of the DEM layers.

• (5) Radius: The search radius used to identify the nearest source point around a random target point. The unit of radius is the same as the linear unit of the DEM layers.

• (6) Workspace: The folder for output rasters and shapefile. To ensure the geoprocessing steps are not affected by any existing files, there should be no existing files or folders in the workspace before a user clicks the OK button; otherwise, a warning message will pop up.

• (7) OK: Click OK to run the program. Error messages will appear if there are any errors in the above items.

The pseudo-code for class ButtonOK(object) is listed below:

Obtain the number of days between two data acquisition dates:

• (1) Extract dates from DEM1 and DEM2 layer names, and calculate the number of days between the two dates. The number of days will be used to convert dune migration distance into migration rate (m/yr).

Process input DEMs:

• for each raster layer in the active data frame:

• if the layer name is DEM1 and the layer is not broken:

Perform the following geoprocessing steps:

• (2) Calculate slope raster from DEM1.

• (3) Convert slope raster to binary raster based on AOR values.

• (4) Remove noises in binary raster.

• (5) Thin binary raster.

• (6) Convert thinned raster to polyline shapefile (polyline1.shp).

If the layer name is DEM2 and the layer is not broken:

Perform the following geoprocessing steps:

• (7) Calculate slope raster from DEM2.

• (8) Convert slope raster to binary raster based on AOR values.

• (9) Remove noises in binary raster.

• (10) Thin binary raster.

• (11) Convert thinned raster to polyline shapefile (polyline2.shp).

Perform the following additional geoprocessing steps:

• (12) Create buffer shapefile (buffer1.shp) around polyline1.shp using search radius.

• (13) Create buffer shapefile (buffer2.shp) around polyline2.shp using search radius.

• (14) Clip polyline1.shp using buffer2.shp to create source line shapefile (sourcelines.shp).

• (15) Clip polyline2.shp using buffer1.shp to create target line shapefile (targetlines.shp).

• (16) Get the maximum polyline length in targetlines.shp to estimate the number of random points on each feature.

• (17) Create random points (tmptargets.shp) along polyline features in targetlines.shp.

• (18) Search the nearest point (NEAR_X, NEAR_Y) on sourcelines.shp for each random point in tmptargets.shp. If the nearest source point is found, calculate distance (NEAR_DIST) between the point pair and the angle. Convert the angle to azimuth direction (field “Direction”) as the migration direction, which is reported by the direction from which dune migration originates. If the nearest source point is not found, the distance value is zero for the random point. Figure 3A shows four selected points in tmptargets.shp and their attributes.

• (19) Add x and y coordinate fields to tmptargets.shp and create a new shapefile (tmptargetsXY.shp).

• (20) Make a copy of tmptargetsXY.shp and save it as tmpsourcesXY.shp. Move each point in tmpsourcesXY.shp to (NEAR_X, NEAR_Y) available in the attribute table.

• (21) Split sourcelines.shp at vertices and create a new shapefile source_split.shp. Calculate azimuth direction (Azimuth1) of each line segment in source_split.shp.

• (22) Create source points “sourcepoints.shp” through spatial join using tmpsourcesXY.shp as target feature and source_split.shp as the join feature. This is to put two fields “Direction” and “Azimuth1” in the same attribute table for comparison.

• (23) Make a copy of sourcepoints.shp and save it as targetpoints.shp. Move each point in targetpoints.shp to (POINT_X, POINT_Y) available in the attribute table.

• (24) Compare field values of Direction and Azimuth1 in the attribute table of targetpoints.shp. If the absolute difference (Diff_Azi) between the “Direction” and “Azimuth1” field values is 90 degrees or 270 degrees, the point in targetpoints.shp is valid; otherwise, the point is deleted. Figure 3B shows three valid target points (5055, 5056, and 5057) after deleting an invalid point in the upper-right corner, and the “Direction” field is for the source direction (migration direction).

• (25) Delete unmatched target points (NEAR_DIST = 0) and unnecessary fields.

• (26) Add a float field “m_Rate” for dune migration rates in the attribute table of targetpoints.shp. Calculate m_Rate values using migration distances in the NEAR_DIST field multiplied by 365, and then divided by the number of days between the two DEMs.

Using the above 26 steps, dune migration directions (source directions) and migration rates are automatically calculated for each random point on the dune slip centerlines. Examples are provided in the next section.

Multitemporal airborne LiDAR data acquired in a 2.4 km by 9 km study area in the White Sands dune field in New Mexico (USA) are used for application demonstration. White Sands dune field is a gypsum-sand dune field that is believed to have occurred during the mid-Holocene and progressed through the deflation of pluvial and playa Pleistocene Lake Otero sediments (Langford, 2003; Kocurek et al., 2007). Current sediment input to White Sands dune field is primarily from the deflation of gypsum from modern playas situated to the southwest of the dune field (Ewing and Kocurek, 2010).

Figures 4 and 5 show application samples using two 1-m-resolution DEMs derived from LiDAR data acquired on 24 January 2009 and 6 June 2010, in the 2.4 km by 9 km study area in White Sands dune field. The parameters are shown on the add-in toolbar, and the legends are shown in the table of contents on the left side of the ArcMap window. Using the add-in, the results for 8707 points in the 2.4 km by 9 km study area (Fig. 5) were obtained in less than 2 min on ordinary desktop computers.

The dune migration rates in Figure 5 are in line with some previous results. For example, Ewing and Kocurek (2010) reported that migration rates for the barchanoid and crescentic dunes in the area range from 1 to 7 m/yr and are typically faster in the upwind section. As can be seen in Figure 5, the crescentic dunes near the upwind western margin have the highest migration rates. The slip faces are less developed in the central part of the study area, which has relatively wider interdune areas, while the relatively stable parabolic dunes form in the east part of the study area due to the presence of vegetation. To compare the dune migration distances obtained from the GIS add-in and from manual measurements, distances were manually measured from 30 target points (targetpoints.shp) to source lines (sourcelines.shp). It should be noted that manual measurement of dune migration distance and migration directions could be very difficult without source lines and target lines detected from slip faces of LiDAR-derived DEMs, due to the lack of obvious ground control points on the DEMs. Before taking the manual measurements, the visibility of source points (sourcepoints.shp) on source lines was set to invisible so that human visual interpretation was not affected by computer-generated source points. Table 1 lists the manually and automatically measured dune migration distances at 30 target point locations. The results in Table 1 suggest that the relative error of manually measured migration distances can be less than 2% compared with automatically measured migration distances. However, it should be noted that manual measurements rely heavily on the experience and judgment of the operator, and they are very time-consuming compared with the GIS add-in, which can automatically measure and save dune migration rates and migration directions for over 10,000 points in less than 3 min.

This paper introduced the development of a free GIS add-in for automated measurement of sand dune migration using LiDAR-derived, multitemporal, high-resolution DEMs, including the research background, method, software implementation, and application demonstration. The results are automatically saved in GIS databases, which allow for further analysis of the results and input parameters. For example, dune migration rates at individual target point locations can be used to create continuous raster surfaces through spatial interpolation to better reveal the dune-field scale patterns of dune migration, and source directions can be used to statistically reveal prevailing wind directions and detect abnormal source directions for further investigation. The application demonstration using multitemporal LiDAR data from the White Sands dune field in New Mexico (USA) shows that the results are in line with landform units of the study area. With the help of source lines and target lines detected from LiDAR-derived DEMs, manual measurement of dune migration distances could provide similar results as the GIS add-in. However, the accuracy of manual measurements depends on the experience and judgment of the operator. For large dune fields such as the White Sands dune field, manual measurement of the migration direction and migration rate at numerous individual locations can be very time-consuming. As the first GIS tool for automated measurement of desert dune migration rates and migration directions using LiDAR-derived DEMs, the GIS add-in can greatly improve the accuracy and efficiency of dune migration measurement. Finally, it should be noted that the add-in is for crescent-shaped (barchans), barchanoid ridge, transverse, and parabolic dunes, not for linear dunes that have two slip faces on both sides caused by acute bimodal winds, or star dunes that have multiple slip faces.

Light detection and ranging (LiDAR) data acquisition and processing were completed by the National Center for Airborne Laser Mapping (NCALM, http://www.ncalm.org). NCALM funding was provided by the National Science Foundation Division of Earth Sciences, Instrumentation and Facilities Program, grant EAR-1043051. LiDAR data sets were downloaded from the OpenTopography Facility at the San Diego Supercomputer Center at http://dx.doi.org/10.5069/G9Q23X5P and http://dx.doi.org/10.5069/G97D2S2D. P.D.’s work was supported by the State Key Laboratory of Remote Sensing Science, Chinese Academy of Sciences (Project OFSLRSS201607). Xia’s work was partly supported by the National Natural Science Foundation of China (grant 41461103). Author contributions: Dong designed the research conception; Xia and Dong performed coding and data analysis; Xia and Dong wrote the paper. The authors thank Mauro Spagnuolo and an anonymous reviewer for their comments and suggestions.