# THE NEED FOR A HIGH-ACCURACY, OPEN-ACCESS GLOBAL DIGITAL ELEVATION MODEL

EDITED BY : Guy Jean-Pierre Schumann and Paul Bates PUBLISHED IN : Frontiers in Earth Science, Frontiers in Environmental Science and Frontiers in Ecology and Evolution

#### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers. The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88966-334-7 DOI 10.3389/978-2-88966-334-7

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# THE NEED FOR A HIGH-ACCURACY, OPEN-ACCESS GLOBAL DIGITAL ELEVATION MODEL

Topic Editors:

Guy Jean-Pierre Schumann, University of Bristol, United Kingdom Paul Bates, University of Bristol, United Kingdom

Citation: Schumann, G. J.-P., Bates, P., eds. (2021). The Need for a High-Accuracy, Open-Access Global Digital Elevation Model. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88966-334-7

# Table of Contents

*04 Editorial: The Need for a High-Accuracy, Open-Access Global Digital Elevation Model*

Guy J.-P. Schumann and Paul D. Bates


Dean B. Gesch


Apoorva Shastry and Michael Durand


Ricardo Tavares da Costa, Paolo Mazzoli and Stefano Bagli


Dietmar J. Backes and Felix Norman Teferle

# Editorial: The Need for a High-Accuracy, Open-Access Global Digital Elevation Model

Guy J.-P. Schumann1,2 \* and Paul D. Bates <sup>1</sup>

1 School of Geographical Sciences, University of Bristol, Bristol, United Kingdom, <sup>2</sup> Research and Education Department, RSS-Hydro, Dudelange, Luxembourg

Keywords: digital elevation model, open-access, global, challenge, accuracy, applications

Editorial on the Research Topic

#### The Need for a High-Accuracy, Open-Access Global Digital Elevation Model

Many different scientific and end-user communities across many disciplines and sectors, as well as many government agencies and international aid and development organizations, agree that there is a pressing need for a high-accuracy, open-access global digital elevation model (DEM). Although some notable efforts and negotiations among important players within the private industry, government agencies and international NGOs are already happening, more needs to be done to realize and speed-up progress in creating a high-accuracy, open-access DEM.

The purpose of this collection of articles is to raise awareness further by describing the problems with current global DEMs and illustrating best practices as well as upcoming opportunities.

In an opening opinion article, Schumann and Bates (2018) lay out the limitations and the requirements for a better global, open-access DEM and present a possible way forward. In a reactive commentary, Winsemius et al. (2019) further argue that for the many applications such a DEM would be game-changing, a DEM itself is not enough and, particularly at the local scale and at the impact-level scale, other properties of the environment besides terrain elevation become important as well. These include details of critical infrastructure such as bridges, culverts and flood defenses. The authors therefore argue that a bottom-up data collection approach at the local community level should complement the global consortium approach suggested by Schumann and Bates (2018) for a better global DEM.

In a thorough assessment study, Tavares da Costa et al. (2019) examine the adequacy of ten different free DEMs for watershed studies. They argue that the intrinsic inaccuracies in free DEMs limit progress and knowledge and they demonstrate that while most DEMs generally represent elevation profiles well enough, they do not adequately represent important topographic and geomorphic features and therefore are inadequate for many applications. This of course has practical implications since non-trivial limitations of any particular global DEM currently available may significantly hinder progress in solving a number of environmental and socioeconomic challenges. In a similar study but focusing on coastal areas, Gesch (2018) discusses best practices for elevation-based assessments of sea level rise and coastal exposure risk studies. Using many different low-accuracy global DEMs, but also high-accuracy local DEMs, he demonstrates how accuracy information should be considered in planning and implementation studies. The work shows that current global DEMs are not adequate for high confidence mapping of exposure to fine increments of sea level rise or with shorter planning horizons and thus they should not be used in this way, but they are suitable for general delineation of low elevation coastal zones.

#### Edited and reviewed by:

Nick Van De Giesen, Delft University of Technology, Netherlands

#### \*Correspondence:

Guy J.-P. Schumann gjpschumann@gmail.com

#### Specialty section:

This article was submitted to Hydrosphere, a section of the journal Frontiers in Earth Science

Received: 16 October 2020 Accepted: 22 October 2020 Published: 12 November 2020

#### Citation:

Schumann GJ-P and Bates PD (2020) Editorial: The Need for a High-Accuracy, Open-Access Global Digital Elevation Model. Front. Earth Sci. 8:618194. doi: 10.3389/feart.2020.618194

In view of the fact that the creation of a high-accuracy, openaccess global DEM may still take some considerable time, Hawker et al. (2018) and Shastry and Durand (2019) present ways to improve upon existing global DEMs. Hawker et al. (2018) suggest the creation of new DEMs using a geostatistical approach to stochastically simulate floodplain DEMs from several openaccess global DEMs based on the spatial error structure. This DEM simulation approach enables an ensemble of plausible DEMs to be created, which can be used for many different applications that are based on probabilistic assessments, such as probabilistic flood risk mapping. In the same context, Shastry and Durand (2019) propose to make use of flood inundation extents from remotely sensed observations to obtain better floodplain topography in data-poor areas, given that such observations indirectly provide information about topography. By combining this information with model predictions via a data assimilation approach, better floodplain topography maps can be generated.

As mentioned by the opening article of Schumann and Bates (2018), there are existing technologies to generate DEMs at the required resolutions and accuracies, over various spatial scales. Backes and Teferle (2020) present a multiscale integration of very high-resolution satellite and drone imagery for creating a highaccuracy DEM over Tristan da Cunha, a remote group of small volcanic islands in the South Atlantic Ocean. Their work demonstrates that combining very high-resolution satellite imagery and low-altitude drone-based imagery can produce inexpensive alternatives to high-quality industry-standard DEMs. This might be particularly relevant for regions, such as small island developing states, where existing satellite data might be insufficient and which may not be prioritized in data acquisition campaigns. Following the same idea of using novel technologies to generate topographic datasets, Faherty et al. (2020) show that with intelligent image classification methods, bare Earth DEMs can be generated for large floodplains using the newest developments in Ka-band synthetic aperture radar (SAR) sensor technology. Specifically, they used NASA GLISTIN-A

#### REFERENCES


airborne mission data with single-pass interferometric SAR (InSAR) processing to derive a terrain model over a large portion of the Red River of the North floodplain (ND, United States), with vertical accuracies comparable to that of state-ofthe-art airborne LiDAR but at a much lower cost.

Using the same SAR sensor characteristics onboard the upcoming NASA/CNES surface water ocean topography (SWOT) satellite mission, Langhorst et al. (2019) give valuable insights of how this new sensor technology can improve our global knowledge about important river hydraulic variables, such as water surface elevations, slope breaks and knickpoints in rivers. They present a novel noise reduction method for multitemporal river water surface elevation profiles from simulated surface water ocean topography data on the Po, Sacramento, and Tanana Rivers and obtain average profiles with errors much lower than those of existing DEMs. Findings like this allow of course new advances in riverine research globally that is not possible with the current low accuracy global DEMs.

Despite the exciting approaches described in this Research Topic, DEM data collection remains under-prioritized by satellite agencies and commercial providers given the importance of these data to a range of disciplines and applications. Fundamental progress at the global level for many key environmental problems will only be achieved when a better accuracy and open-access global DEM becomes available.

#### AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

#### ACKNOWLEDGMENTS

We would like to thank all the authors for contributing to this collection.


Conflict of Interest: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Schumann and Bates. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The Need for a High-Accuracy, Open-Access Global DEM

Guy J-P. Schumann1,2,3 \* and Paul D. Bates <sup>1</sup>

<sup>1</sup> School of Geographical Sciences, University of Bristol, Bristol, United Kingdom, <sup>2</sup> Remote Sensing Solutions Inc., Monrovia, CA, United States, <sup>3</sup> RSS-Hydro S.à.r.l.-S, Innovation Hub Dudelange, Dudelange, Luxembourg

Keywords: DEM (digital elevation model), hazard & risk, accuracy, open access, global–local processing

# LIMITATIONS OF CURRENT OPEN-ACCESS DEMS

Digital elevation models (DEM) are recognized as a core spatial dataset required for many environmental applications. However, the availability of comprehensive DEMs for water resources studies is rather limited and limitations of current, free or open-access DEMs are well-known. Freely available and global scale DEMs, such as that from the Shuttle Radar Topography Mission (SRTM) or from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) mission exhibit large vertical errors which are exacerbated over complex topography and they cannot resolve microtopographic variations in relatively flat terrain (Gallien et al., 2011; Chu and Lindenschmidt, 2017). For instance, SRTM mission requirements defined absolute and relative elevation errors of 16 and 6 m, respectively (Rabus et al., 2003). Even though several studies have found actual errors to be considerably smaller than the requirements (see e.g. Schumann et al., 2008; Patel et al., 2016), the inherent errors in the SRTM-DEM are still notably unacceptable. Over the years, several processing algorithms and approaches for merging with other elevation datasets have been proposed to increase accuracy and remove vegetation biases (Baugh et al., 2013; Robinson et al., 2014; O'Loughlin et al., 2016; Yamazaki et al., 2017; Yue et al., 2017). However, whilst such derived versions are widely used, they still typically exhibit errors in the vertical much larger than those acceptable for many applications.

For instance, current global DEMs cannot resolve the detail of terrain features that control flooding (Schumann et al., 2014). Kenward et al. (2000) assessed the effects of different largescale DEMs on hydrologic runoff predictions and highlighted that different DEMs can lead to a difference of almost 10% in runoff predictions. Sanders (2007) examined the suitability of different online DEMs for flood modeling, including the US National Elevation Dataset (NED), Shuttle Radar Topography Mission (SRTM) data, airborne interferometric synthetic aperture radar (InSAR) data, and the Light Detection and Ranging (LIDAR) data. The study concluded that DEMs based on InSAR technology, including SRTM, suffer from radar noise and require prior processing while the NED DEM was remarkably smooth and may overestimate inundation extent. Nevertheless, Sanders (2007) highlighted the utility in SRTM as a global source of terrain data for flood modeling, as successfully demonstrated by e.g., Sampson et al. (2015). In a similar study, Li and Wong (2010) argued further that inundated areas from a flood simulation, albeit a simplistic case in their study, vary significantly across different DEM data sources. In another study dedicated to hydrologic applications, Mukherjee et al. (2013) evaluated the vertical accuracy of open source DEMs (ASTER and SRTM) and concluded that slope and drainage network delineation are largely violated compared to their reference DEM. As illustrated by Walker and Willgoose (1999), even higher accuracy published DEMs from aerial photogrammetry may be significantly different from ground truth, especially when looking at smaller catchments or applications for which localized errors in the elevation can have detrimental effects on local scale applications.

Another serious limitation of most freely available global DEMs, such as the SRTM-DEM acquired in February 2000, is that they are now very dated. Acquisition is not routinely repeated,

#### Edited by:

Nils Moosdorf, Leibniz Centre for Tropical Marine Research (LG), Germany

#### Reviewed by:

Norbert Pfeifer, Vienna University of Technology, Austria

> \*Correspondence: Guy J-P. Schumann gjpschumann@gmail.com

#### Specialty section:

This article was submitted to Hydrosphere, a section of the journal Frontiers in Earth Science

Received: 22 August 2018 Accepted: 23 November 2018 Published: 04 December 2018

#### Citation:

Schumann GJ-P and Bates PD (2018) The Need for a High-Accuracy, Open-Access Global DEM. Front. Earth Sci. 6:225. doi: 10.3389/feart.2018.00225

**6**

whereas most national DEM programs acquiring photogrammetric or LiDAR DEMs, are repeated at regular multi-annual intervals. Consequently, any significant geomorphological change in topography as well as any man-made alteration of the Earth surface since the time of data acquisition will not be captured in global DEMs. Of course, this can have significant consequences when predicting for instance flood events in locations where previous events are known to have altered local topography considerably. More strikingly, there are now nearly 1.5 billion more people on earth compared to when SRTM was acquired, and every single one of them will have contributed to re-shaping the surface topography of our planet in some way. In this context, James et al. (2012) for instance noted that important geomorphological changes can occur over decadal timescales, and even more rapidly during extreme events and significant human alterations of landscapes (e.g., mining etc.).

It is clear that currently available global DEMs cannot be used to accurately simulate or predict local scale processes or impacts thereof. Even after considerable preprocessing to remove significant biases (due to vegetation and other physical structures) and to reduce inherent vertical errors, publicly available global DEMs still suffer from inaccuracies oftentimes orders of magnitude greater than length scales of the processes that are simulated. For instance, the residual vertical error in the SRTM DEM is for most river gradients orders of magnitude higher than the magnitude of the flood waves in those rivers (Bates et al., 2013). Consequently, it becomes very difficult to accurately model floods, and other inherently local phenomena, with current coarse resolution global DEMs, such as the SRTM-DEM. Local-scale skill is key, and this needs higher resolution and accuracy terrain data than currently available at large coverage. Furthermore, large scale hazard models are typically validated for rural areas, but all the major impacts are felt in urban zones where relevant process length scales are much higher; however, it is in urban areas where current global DEMs break down.

# REQUIREMENTS OF A GLOBAL HIGH-ACCURACY OPEN-ACCESS DEM

#### Level of Accuracy

Spatial resolution, absolute error in the vertical, and accuracy in relative gradient (slope) are all important attributes of a DEM, and the requirements of those need to be properly established. **Figure 1** puts available DEM technology in context of the requirements for spatial resolution and vertical accuracy for different scales of application.

Whilst views on the impact of spatial resolution may be more varied, depending on the complexity of the process simulated or indeed the DEM source (Horritt and Bates, 2001; Fewtrell et al., 2008; Li and Wong, 2010), absolute and relative accuracy in the vertical is crucial and even small vertical errors can have significant impacts on the application accuracy, especially at the local scale.

LeFavour and Alsdorf (2005) showed that for calculating a reliable slope from SRTM, reach lengths need to extend sufficient distances to accommodate the height errors. They developed the following simple relationship to determine the appropriate reach length (RL):

$$2\sigma/\text{RL} = \text{S}\_{\text{min}} \tag{1}$$

For example, to be able to reproduce the minimum slope, Smin, along the very shallow gradient mainstem of the Amazon River of 1.5 cm/km found by Birkett et al. (2002) using Topex/POSEIDON radar altimetry, the appropriate RL would need to be of 733 km in order to accommodate the ± 5.51 m height error that LeFavour and Alsdorf (2005) derived. In a similar study of the River Po in Northern Italy, Schumann et al. (2010) showed that when using an adequate RL, SRTM-DEM can be used to estimate river gradients almost as reliably as LiDAR data, albeit with much higher levels of uncertainty in the data.

Supposedly, Airbus Defense and Space's WorldDEMTM based on the TanDEM-X data, provides a global DEM of unprecedented quality, accuracy, and coverage with a vertical accuracy of 2 m (relative) and better than 6 m (absolute) at a 12 m spatial resolution (Riegler et al., 2015). While it is true that the accuracy may surpass that of any global satellite-based elevation model available, WorldDEM is not a free, open-access DEM and probably will therefore not define a new geospatial standard in global elevation models. Also, the errors stated, both in the absolute and relative are hardly sufficient for local-scale applications requiring an accuracy in the vertical of at least 0.5 m (**Figure 1**). Nonetheless, a number of recent studies have demonstrated that with expert processing, TanDEM-X data can be successfully employed to model flood hazard (Mason et al., 2015, 2016; Gutenson et al., 2017). It is worth noting that whilst this DEM may be more accurate than the SRTM-DEM, it is (also) not a bare-earth DEM and so deriving an accurate terrain model from the TanDEM-X source data would require significant time and resources.

Vricon (owned by Saab and DigitalGlobe) produces a very high resolution near-global digital surface and terrain model (DSM, DTM) at 0.5 m spacing using commercial multi-pair satellite imagery from DigitalGlobe, with a vertical absolute accuracy of 3 m and 1 m relative accuracy, according to the data provider (Vricon, 2018). They also generate a 10 m DSM of the same accuracy specifications from the 0.5 m resolution DSM. Again, such a near-global DEM could be game-changing for many applications but is currently not available as an open-access product. As a highly successful open-access flagship example using this type of imagery and technology, the ArcticDEM<sup>1</sup> is an NGA-NSF public-private initiative to automatically produce a high-resolution, high quality digital surface model (DSM) of the Arctic using optical stereo imagery, high-performance computing, and open source photogrammetry software.

# Different Data Acquisition Technology

As illustrated in **Figure 1**, many different proven technologies exist to acquire DEMs, however, more importantly, there is

<sup>1</sup>https://www.pgc.umn.edu/data/arcticdem/

no need for a new technology in order to meet the accuracy requirements shown, even for local-scale applications. This has also been argued by Schumann et al. (2014) who suggested that an advanced global DEM would use existing LiDAR data and stereo satellite images, and new LiDAR elevation data could be acquired on board disaster-relief aircraft or on drones deployed over flood plains. The operation costs would therefore be substantially cheaper than most satellite missions. Schumann et al. (2016) further stated that such a DEM could be achieved from a consortium of industry, governments and humanitarian agencies and would undoubtedly be the environmental equivalent of the Human Genome Project (HGP; NIH, 2010) but at a fraction of the cost (see Sampson et al., 2016).

Hence, the innovation of creating a high-accuracy, open-access global DEM would merely lie in merging and quality controlling different "tiles" from different acquisition technologies and, prior to that, establishing the same metadata for those tiles, on a standardized (vertical) datum and map projection, such as for instance the recently proposed Discrete Global Grid System<sup>2</sup> (DGGS). The goal of DGGS is to enable rapid assembly of spatial data without the difficulties of working with projected coordinate reference systems. DGGS would provide the capability to properly integrate global geospatial, social, and economic information. It also allows communities with data attributed to fundamentally different geographies to integrate this information in a single consistent framework. Also, new satellite sensor technology, such as on ICESat-2 (https://icesat-2.gsfc.nasa.gov), provides unprecedent elevation measurement accuracy which may be used to assist vertical datum matching and rectification (see Hall et al. (2012) for an example of this type of application with ICESat-1).

# Storage Facility

At the White House Climate Data Initiative announcement in March 2014<sup>3</sup> , Google Inc. and many other industry entities, government agencies and NGOs, took significant steps in the right direction. Google Inc. committed one petabyte of cloud storage for climate data, as well as 50 million hours of highperformance computing and challenged the global innovation community to develop a high-accuracy global terrain model to help communities build resilience to anticipated climate impacts.

Although this private industry commitment is absolutely fantastic, it should be noted that privately-owned platforms and servers, even if made completely open with free registration, may not be accessible to all given that there may be access restrictions in place in some counties or at some governmental agencies. Therefore, it may be better to store a high-accuracy, open-access DEM on a fully accessible public server. Having said that, these are legalities that can eventually be sorted out easily, even with some flexibility from private sector industries.

#### Licensing

A high-accuracy, global DEM needs to be open-access and as such needs to be under a fully open license. Broadly speaking, an

<sup>2</sup>http://docs.opengeospatial.org/as/15-104r5/15-104r5.html

<sup>3</sup>https://obamawhitehouse.archives.gov/blog/2014/03/19/climate-data-initiativelaunches-strong-public-and-private-sector-commitments

open license is one which grants permission to access, re-use and redistribute a work with no restrictions. Current low resolution and low accuracy global DEMs, such as the SRTM or ASTER are already distributed under a fully open license. Equally, many governments, in Europe, the U.S., Australia and elsewhere have already adopted an open license for national LiDAR elevation data, so extending this principle to a high-accuracy global DEM is easily achievable.

#### WAY FORWARD

Along the lines of the argument put forward by Schumann et al. (2014) to form a consortium of many different organizations working in a public-private partnership to create an open-access high-accuracy global DEM, McDougall et al. (2008) already suggested several years earlier that a collaborative approach has the potential to deliver improved outcomes with respect to the development of an open-access, high-quality DEM. A collaborative approach should optimize the existing investment across the public and private sectors and reduce potential duplication of effort.

It should be noted that we are not suggesting any highly laborious efforts, such as for instance, a country-by-country collection of already existing national scale high-resolution DEMs, but rather we suggest a consortium effort to look for a feasible solution using DEMs from already available

#### REFERENCES


commercial high-resolution technology, preferably merged with LiDAR (or other technology of similar capability) DEMs where already available as open data. The resulting DEM should be global and be made open-access, so that it can be more effectively utilized to support decision making on important issues such as water management and disaster assessment and mitigation. Negotiations among important players within the private industry, government agencies and international NGOs are already well underway [see for instance Schumann et al. (2014) and associated online comment<sup>4</sup> ] but still there needs to be more substantial efforts to realize and speed-up progress in creating a high-accuracy, open-access DEM.

#### AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

### ACKNOWLEDGMENTS

All opinions expressed in this article are the opinions of the authors and do not represent any government's or other entity's opinion.

<sup>4</sup>https://www.nature.com/articles/507169e


DEM product. Remote Sens. Environ. 182, 49–59. doi: 10.1016/j.rse.2016. 04.018


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Schumann and Bates. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Commentary: The Need for a High-Accuracy, Open-Access Global DEM

Hessel C. Winsemius 1,2 \*, Philip J. Ward<sup>3</sup> , Ivan Gayton<sup>4</sup> , Marie-Claire ten Veldhuis <sup>1</sup> , Didrik H. Meijer 1,2 and Mark Iliffe<sup>5</sup>

*<sup>1</sup> Department of Water Management, Faculty of Civil Engineering and Geosciences, Delft University of Technology, Delft, Netherlands, <sup>2</sup> Deltares, Delft, Netherlands, <sup>3</sup> Institute for Environmental Studies, Vrije Universiteit Amsterdam, Amsterdam, Netherlands, <sup>4</sup> Humanitarian OpenStreetMap Team, Washington, DC, United States, <sup>5</sup> United Nations Statistics Division, Department of Economic and Social Affairs, United Nations, New York, NY, United States*

Keywords: DEM (digital elevation model), open-access, community mapping, urban, flood risk

#### **A Commentary on**

#### **The Need for a High-Accuracy, Open-Access Global DEM**

by Schumann, G. J.-P., and Bates, P. D. (2018). Front. Earth Sci. 6:225. doi: 10.3389/feart.2018.00225

#### Edited by:

*Ke Zhang, Hohai University, China*

#### Reviewed by:

*Simone Tarquini, National Institute of Geophysics and Volcanology (INGV), Italy*

> \*Correspondence: *Hessel C. Winsemius h.c.winsemius@tudelft.nl*

#### Specialty section:

*This article was submitted to Hydrosphere, a section of the journal Frontiers in Earth Science*

Received: *16 January 2019* Accepted: *13 February 2019* Published: *06 March 2019*

#### Citation:

*Winsemius HC, Ward PJ, Gayton I, ten Veldhuis M-C, Meijer DH and Iliffe M (2019) Commentary: The Need for a High-Accuracy, Open-Access Global DEM. Front. Earth Sci. 7:33. doi: 10.3389/feart.2019.00033*

Schumann and Bates (2018), hereafter referred to as SB, argue that there is a need for a global scale, high-accuracy, open-access Digital Elevation Model (DEM). They find current global DEM data unacceptable in accuracy and resolution. They position their argument mainly around hydraulic applications, by referring to required accuracies to estimate hydraulic energy slopes, and the need for highly detailed topography in urban zones, where most flood impacts are felt. A vertical accuracy of 0.5 meter is suggested to be adequate for local scale applications. This DEM, according to SB, can be established through a collaborative effort of industry, governments and humanitarian agencies by merging and quality assuring DEMs from different existing sources such as Light-Detection and Ranging missions, photogrammetry or satellite stereo. They close their plea by stating that the anticipated DEM will "help communities build resilience to anticipated climate impacts."

We support the statement that a DEM with the envisioned accuracy is useful. In fact, it is a general prerequisite to simulating flood behavior at large scales, as indicated by SB. In addition, terrain and elevation data at this level of detail may find its application in other fields such as morphology (e.g., Tarolli and Sofia, 2016), cadastral digitization (e.g., Dorninger and Pfeifer, 2008) and landslide predictions (e.g., Ciampalini et al., 2016). However, at the level of detail suggested (urban zones), an accurate DEM alone is not enough if the objective of flood modeling goes beyond awareness raising or flood zoning. In this commentary, we argue that: (a) the "consortium effort" proposed by SB should focus on regions that will profit most: developing countries; (b) technically, at the local scale where the suggested accuracy (0.5 m) with the suggested spatial resolution (∼5 m) becomes useful, other properties of the environment besides terrain become dominant in flood behavior and, consequently, in flood risk. These include the water infrastructure (bridges, channels, culverts, etc.) and its maintenance state. Such infrastructure information is key in establishing what makes a street, neighborhood or city more flood resilient, hence supporting communities in their decision-making for a sustainable future; and (c) that a bottom-up data collection approach for infrastructure and maintenance states will lead to less inequality in global data coverage

and a reduction in dependency on outside capacities. This then leads to sustainable data collection and ownership in the regions where flood risk information is most needed. We elaborate on these three points below.

First and foremost, the efforts should concentrate on the regions that may profit most (in socio-economic terms) from an open-access DEM, provided it is kept up to date, as SB suggest. These are rapidly urbanizing cities in developing countries, where no local resources are at hand to regularly establish a detailed DEM. Globally, urban population already accounts for over 50% of total population. In developing countries this number is still below 50% but the expectation is that it will rapidly catch up with an expected 90% growth in urban population happening in Africa and Asia until 2050 (United Nations, 2018). Furthermore, in these areas, poor people may be affected to a degree that it traps them into poverty (Winsemius et al., 2018) and therefore there is a strong rationale to link climate and poverty policies (Hallegatte et al., 2016) in order to achieve global sustainability targets.

Within these urban areas in developing countries, flood protection is very limited (Scussolini et al., 2016). Hence pluvial floods can cause frequent small-scale impacts, and river floods will already occur with flow in excess of natural bank-full capacity. According to extreme value theorem, this occurs at the very least once every 1.5–2 years (i.e., the most occurring annual flood, see e.g., Dunne and Leopold, 1978; Savenije, 2003). It is likely, that due to rapid urbanization, overbank flow will occur even more frequently than once every 1.5–2 years, due to decreased soil water holding capacity, rapid erosion and impediments in drainage due to sediment and solid waste deposition around and in front of infrastructure. It is these moderate floods that cause the largest contribution to a community's accumulated risk, because they simply occur most frequently (Veldhuis et al., 2011). In urban environments, these moderate events are much more dictated by details, i.e., local infrastructural impediments (e.g., lack of drainage connectivity, open and closed tertiary drainage, culverts) and their maintenance state (drainage impediments by sedimentation, solid waste accumulation) than solely the surrounding topography (e.g., Veldhuis, 2011; Rodriguez et al., 2013). We see all such issues occurring in our experience in working in for instance Dar es Salaam, Accra, Jakarta and Nairobi. Besides quality of flood simulations, the ability to assess the drainage system and effectiveness of interventions in urban systems is essential. For this, the status-quo, i.e., the present-day drainage infrastructure and its maintenance state are essential to understand, and we even argue that this requires attention first and foremost.

Finally, we strongly advocate that the rather top-down approach to data collection, suggested by SB, is complemented by bottom-up approaches, supported by local initiatives and incentives. The top-down approach can offer a first important short-term alleviation of data poverty. However, it is difficult to sustain this, especially if international commercial considerations are at stake. Why would a company such as Google, map a poor vulnerable neighborhood, including its water infrastructure when there is little commercial gain to be expected? Furthermore, top-down implies that capital is brought in from outside rather than building capacity to generate capital on the inside. Human capital, expressed in raising of knowledge and education levels, is the main corner stone of prosperity and economic growth in the long-term (Solow, 1956; Piketty, 2014) and therefore, a locally sustainable approach to data collection should be sought to accompany global data collection efforts. For instance, data collection through community mapping is now moving from efforts deployed ex-post, i.e., responding to a crisis (Hagen, 2010; Soden and Palen, 2014) to ex-ante, exemplified by the largest community mapping project in the world "Ramani Huria" (Swahili for "Open Map"), taking place in Dar es Salaam—Tanzania. This project aims to build preparedness and resilience by letting citizens map their own environment (Iliffe et al., 2017). These methods are still in their infancy, and sustainable incentivization is a large challenge. But first results show that if communities are organized, detailed and accurate information on infrastructure, such as roads, bridges, channels, culverts, man-holes along with their dimensions and maintenance states can be collected in a very short time and with very few resources<sup>1</sup> , <sup>2</sup> Ramani Huria demonstrates that this can work by empowering academics through curriculum development, gradually leading to expertise within authorities as well. **Figure 1** illustrates how this can result in a plethora of data that is valuable from a sustainability and humanitarian point of view, even if it has relatively little direct commercial interest. The left-hand side shows the status of Google Maps (on 4 January 2019) vs. the right-hand showing OpenStreetMap at the same date. Recent floods led to water authorities now accepting the use of these data, even though they were not collected by their own staff. Ramani Huria can serve as an example for water authorities elsewhere, demonstrating that bottom-up data collection processes lead to good and even more complete data than what can be achieved by authority data collection efforts. The scientific community can facilitate bottom-up data collection approaches such as community mapping through two research directions: first, to find social structures to promote and create incentives to data collection (Buytaert et al., 2014; Paul et al., 2018); and second through new innovative, affordable and sustainable instrumentation (van de Giesen et al., 2014; Tauro et al., 2018), and local quality assurance methods.

Concluding, a global scale higher accuracy DEM with the envisioned accuracy will be very useful both to science and practitioners. But to truly assist the most critically targeted environments in their road to a sustainable future, elevation data have to be accompanied by information on infrastructure, including its dimensions and state of maintenance. We fully agree with SB that targeted regions should be urban areas, as these accommodate most of the global population and economy, and we advocate that focus should lie on developing countries, where the urbanization rates are the highest and data poverty

<sup>1</sup>http://ramanihuria.org/drain-mapping-flood-modelling/

<sup>2</sup>https://medium.com/@h.c.winsemius/community-mapping-for-floodmodelling-2-0-a3b794e02d16

contributed to OpenStreetMap and therefore plotted separately. Top: Map data © 2019 Google; Bottom: Map data © OpenStreetMap contributors, CC BY-SA.

is most extreme. To truly increase resilience here, the topdown data collection process should be considered an important short-term data poverty alleviation step, which needs to be accompanied by a more sustainable, bottom-up procedure to establish the continued gathering of critical information for flood resilience.

#### AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

#### REFERENCES


#### FUNDING

PW is supported by a VIDI grant (016.161.324) from the Netherlands Organisation for Scientific Research (NWO).

#### ACKNOWLEDGMENTS

We thank the World Bank for the continuing support in community mapping efforts in Dar es Salaam, on which parts of this article have been based. The opinions reflected in this article do not represent the opinion of any government or other entity.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Winsemius, Ward, Gayton, ten Veldhuis, Meijer and Iliffe. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Best Practices for Elevation-Based Assessments of Sea-Level Rise and Coastal Flooding Exposure

#### Dean B. Gesch\*

U.S. Geological Survey, Earth Resources Observation and Science Center, Sioux Falls, SD, United States

Elevation data are critical for assessments of sea-level rise (SLR) and coastal flooding exposure. Previous research has demonstrated that the quality of data used in elevationbased assessments must be well understood and applied to properly model potential impacts. The cumulative vertical uncertainty of the input elevation data substantially controls the minimum increments of SLR and the minimum planning horizons that can be effectively used in assessments. For regional, continental, or global assessments, several digital elevation models (DEMs) are available for the required topographic information to project potential impacts of increased coastal water levels, whether a simple inundation model is used or a more complex process-based or probabilistic model is employed. When properly characterized, the vertical accuracy of the DEM can be used to report assessment results with the uncertainty stated in terms of a specific confidence level or likelihood category. An accuracy evaluation has been conducted of global DEMs to quantify their inherent vertical uncertainty to demonstrate how accuracy information should be considered when planning and implementing a SLR or coastal flooding assessment. The evaluation approach includes comparison of the DEMs with high-accuracy geodetic control points as the independent reference data over a variety of coastal relief settings. The global DEMs evaluated include SRTM, ASTER GDEM, ALOS World 3D, TanDEM-X, NASADEM, and MERIT. High-resolution, high-accuracy DEM sources, such as airborne lidar and stereo imagery, are also included to give context to the results from the global DEMs. The accuracy characterization results show that current global DEMs are not adequate for high confidence mapping of exposure to fine increments (<1 m) of SLR or with shorter planning horizons (<100 years) and thus they should not be used for such mapping, but they are suitable for general delineation of low elevation coastal zones. In addition to the best practice of rigorous accounting for vertical uncertainty, other recommended procedures are presented for delineation of different types of impact areas (marine and groundwater inundation) and use of regional relative SLR scenarios. The requirement remains for a freely available, highaccuracy, high-resolution global elevation model that supports quantitative SLR and coastal inundation assessments at high confidence levels.

#### Edited by:

Paul Bates, University of Bristol, United Kingdom

#### Reviewed by:

Juan Pablo Rodríguez Sánchez, University of Los Andes, Colombia Hongkai Gao, East China Normal University, China

> \*Correspondence: Dean B. Gesch gesch@usgs.gov

#### Specialty section:

This article was submitted to Hydrosphere, a section of the journal Frontiers in Earth Science

Received: 31 July 2018 Accepted: 28 November 2018 Published: 12 December 2018

#### Citation:

Gesch DB (2018) Best Practices for Elevation-Based Assessments of Sea-Level Rise and Coastal Flooding Exposure. Front. Earth Sci. 6:230. doi: 10.3389/feart.2018.00230

Keywords: inundation, digital elevation model, uncertainty, accuracy, remote sensing, vulnerability, topography

**15**

# INTRODUCTION

feart-06-00230 December 11, 2018 Time: 16:29 # 2

The effects of sea-level rise (SLR) and other sources of increased water levels along the world's coastlines are pervasive and varied (Williams, 2013). Because of the low-lying nature of many coastal lands, the topography, or elevation in relation to sea level, largely controls their exposure to adverse effects of increased water levels, both chronic conditions (SLR) and episodic events (storm surge inundation or king tide flooding). Elevation data, often in the form of digital elevation models (DEMs), are therefore critical for assessments of exposure, and corresponding vulnerability and risk, to permanent or temporary flooding and other effects of increased water levels along the coast.

For nearly four decades, elevation data, most often in the form of DEMs, have been used to identify low-lying coastal lands over broad areas to conduct assessments of the effects of rising sea levels (Schneider and Chen, 1980; Titus et al., 1991; Titus and Richman, 2001; Small and Nicholls, 2003; Ericson et al., 2006; Rowley et al., 2007; Dasgupta et al., 2008, 2010; Weiss et al., 2011; Haer et al., 2013; Blankespoor et al., 2014; Neumann et al., 2015; Hardy and Nuse, 2016; Kulp and Strauss, 2017; Small et al., 2018; Wolff et al., 2018). Over this period, the quality of DEMs available for use in assessments has improved, especially in terms of spatial resolution and vertical accuracy. For large-area assessments (regional, continental, global), several choices are available for DEMs for the required topographic information to project potential impacts of increased coastal water levels, whether a simple inundation model is used or a more complex process-based or probabilistic model is employed. Previous research has demonstrated that the quality of data, and associated transformations, used for elevation-based assessments must be well understood and applied to properly model potential impacts (Gesch, 2009; Coveney and Fotheringham, 2011; Cooper and Chen, 2013; Cooper et al., 2013, 2015; Gesch, 2013; Schmid et al., 2014; Dahl et al., 2017; Jones et al., 2017; West et al., 2018).

Uncertainty in climate change and coastal hazard assessments has been addressed in numerous studies. Le Cozannet et al. (2015, 2017) and Stephens et al. (2017) discuss the importance of considering uncertainties from all sources in assessments of sea-level change, but they do not specifically mention vertical uncertainty of topographic data used to map potential impact zones on the land surface. A number of studies address the effects of input elevation data uncertainty (Gesch, 2009; Kettle, 2012; Bell et al., 2014; Hinkel et al., 2014; Valentine, 2014) and conclude that higher resolution data with better vertical accuracy significantly improve assessment results. The choice of elevation model for an assessment study can also have a substantial effect on results owing to combined uncertainties of input datasets, especially elevation and population distribution (Lichter et al., 2011; Mondal and Tatem, 2012; Wolff et al., 2016). The result of many SLR assessments is a set of maps that spatially show the areas exposed to inundation or other adverse effects of specific scenarios of sea-level change, and such maps are enhanced by including a description of mapping uncertainty (Kostelnick et al., 2013; Retchless, 2018), often expressed as a confidence level.

Recently, within the coastal hazard modeling community there has been continued recognition that vertical uncertainty should be accounted for quantitatively to improve impact mapping and assessment. National Oceanic and Atmospheric Administration [NOAA] (2010) and Doyle et al. (2015) recognize that DEM vertical accuracy is a critical element in SLR impact assessments and that the uncertainties contributed by elevation data and its associated transformations should be explicitly addressed. Fortunately, there is a rich heritage of technical work on assessing vertical uncertainty (DEM error), its consequences and implications, and its use in improving applications (Hunter and Goodchild, 1995; Fisher and Tate, 2006; Wechsler and Kroll, 2006; Maune et al., 2007; Wechsler, 2007; Höhle and Höhle, 2009), and this body of work can be relied upon as a basis for approaches to rigorous handling of uncertainty in coastal inundation assessments.

The importance of the quality (spatial resolution and vertical accuracy) of the input elevation information in inundation assessments has been well recognized in some studies (Bales and Wagner, 2009; Coveney and Fotheringham, 2011; Zhang, 2011; Fraile-Jurado and Ojeda-Zújar, 2012; Sampson et al., 2016; Wolff et al., 2016; Mogensen and Rogers, 2018; Paprotny et al., 2018; Vousdoukas et al., 2018); however, numerous other studies make no mention of the uncertainty of the critical input elevation layer and the implications for the reliability of the results. There is a variety of DEMs with global or near global coverage mostly derived from remote sensing (Gesch, 2012b; Sampson et al., 2016) available for use in inundation modeling and assessment, yet better data are needed (Simpson et al., 2015). Therefore, it is critical that producers of such assessments understand, characterize, and describe the cumulative uncertainties from the underlying elevation data and associated transformations that propagate to the assessment results (maps and estimates of impacted area and resources). The purpose of this paper is to document the best practices in properly accounting for the vertical uncertainty inherent in all elevation-based SLR and coastal flooding assessments. Additionally, other best practices for such assessments extracted from the scientific record of successful studies are listed and described. Establishment of best practices, or standardized methodology, is recognized as an important advance in improving the usefulness of climate change vulnerability mapping at various scales and for multiple classes of stakeholders and end users (Preston et al., 2011).

# MATERIALS AND METHODS

When an elevation-based SLR or coastal flooding assessment is conducted, especially over broad areas, a simple inundation method known as the "bathtub" model is often used, an approach that has also been referred to as the "single-value surface" (National Oceanic and Atmospheric Administration [NOAA], 2010), "equilibrium" (Gallien et al., 2011), "planar" (Bates and De Roo, 2000), "hydrostatic" (Habel et al., 2017), and "static inundation" (Paprotny et al., 2018) method. To delineate the inundation zone with the bathtub method, the water level is simply raised on a coastal DEM by selecting all areas that are below the specified new water level height. The approach is improved by enforcing hydrologic connectivity

(Poulter and Halpin, 2007; Poulter et al., 2008), ensuring that flooded areas have a direct hydrologic connection to the ocean, which is recommended as a best practice for coastal assessments. Limitations of the bathtub modeling approach have been identified (Passeri et al., 2015; Boyd et al., 2016), including failure of the DEM to represent the detailed hydraulic connections and barriers needed for accurate spatially explicit flood mapping (Gallien et al., 2011, 2013). Bathtub modeling can overpredict flood extent compared to hydraulic and hydrodynamic modeling approaches (Gallien et al., 2011, 2014; Seenath et al., 2016), especially at local scales, yet the simple approach realizes savings in input data and computation requirements (Kovanen et al., 2018) and under certain conditions can perform nearly as well as more complex models (Bates and De Roo, 2000). Another important consideration in using simple inundation models to assess potential impacts of raised coastal water levels is to recognize that not all areas will respond to increased water levels by simply becoming inundated (Gesch et al., 2009; Passeri et al., 2015; Lentz et al., 2016). Instead, some areas will adapt by responding dynamically. In these cases, maps of potential impact zones derived from simple bathtub models will indicate areas that may be affected by sea-level change, but not necessarily transition to permanent open water.

Notwithstanding the known limitations of bathtub modeling, the approach remains widely employed for coastal impact assessments, likely due to its ease of application, especially for initial screening and inventory across large areas. Thus, the analysis and discussion below on properly accounting for vertical uncertainty is directly relevant to simple inundation modeling. However, it is also applicable to more complex physical processbased and probabilistic models because they too invariably need to measure increased water levels, and elevation data are required to do so; therefore, vertical uncertainty must be considered.

# Accounting for Uncertainty in Exposure Assessments

The importance of considering uncertainty in general in climate change assessments, and more specifically vertical uncertainty in SLR and coastal flooding studies of interest here, is well established in the previously cited literature. There is a long record of research on uncertainty in geospatial data, with much of it focused on DEMs and other forms of elevation data (Hunter and Goodchild, 1995; Fisher, 1998; Fisher and Tate, 2006; Wechsler and Kroll, 2006; Wechsler, 2007; Höhle and Höhle, 2009). The approaches to handling vertical uncertainty can be categorized into three methods: (1) ignore it; (2) apply a global error estimate; (3) model the error distribution and then perform spatial error propagation through simulation. Each of these categories is described more fully below. Hunter and Goodchild (1995) use the same construct of three general approaches, and they illustrate by use of an elevation contour example.

Most elevation-based SLR assessments mention the DEM used, but many stop there and ignore the inherent vertical uncertainty or the common description of such, the vertical accuracy (or error). The user of such an assessment is left to guess about the quality of the results, and in the case of a user with little familiarity with elevation data, the implications of vertical error are completely unrepresented and thus cannot be factored into decision making. Some studies at least mention the vertical error in the underlying elevation data, and perhaps generally discuss its implications, but make no attempt to quantify spatially how the uncertainty reflected in that vertical error affects the results (Mcleod et al., 2010; Emrich and Cutter, 2011; Kuhn et al., 2011; Weiss et al., 2011; Haer et al., 2013, 2018; Maloney and Preston, 2014). The assessments that fall into this first category do not explicitly consider uncertainty, and they are labeled as "deterministic" as the delineation of the impact zone has no indication of the quality of that mapping, nor is there any expression of confidence that is associated with the results. The location and extent of the inundation zone is determined simply by where the chosen elevation of the raised water level falls on the landscape.

The second category of approaches to handling vertical uncertainty includes assessments that apply a global error metric, such as the widely used root mean square error (RMSE) or a related measure such as "linear error at 95% confidence" (LE95) (Maune et al., 2007). This method equally applies the full global error estimate everywhere, which assumes that all areas are subject to the full range of vertical error. Thus, this approach can be thought of as a worst-case scenario, and the results reflect a range incorporating the minimum and maximum extremes of error. In practice, the full error is applied both above and below a specified elevation (usually representing a raised water level) by adding and subtracting it to the elevation, respectively, and then using those two new elevations in bathtub modeling to delineate the maximum and minimum impact zones (Gesch, 2013). In essence, each delineation is still a deterministic mapping, thus this approach is called here the "modified deterministic" method. It has the advantage over the straight deterministic method in that it addresses uncertainty by bounding the error range and assigning a label of quantified confidence based on the portion of the full error probability distribution represented by the error metric applied, for instance 68% confidence in the case of RMSE or 95% confidence in the case of LE95 (for an unbiased normal distribution of errors). For users of such assessment results, the stated confidence level indicates how confident the user can be that the true extent of the impact zone is contained within the given range (between the minimum and maximum areas). Examples of the successful use of the modified deterministic approach are found in Gilmer and Ferdaña (2012), Gesch (2013), Nielsen and Dudley (2013), and Enwright et al. (2015).

The third category of approaches to accounting for vertical uncertainty includes methods that model the elevation error distribution and then propagate that error spatially through Monte Carlo simulation (Temme et al., 2009). The result of such an operation is a map containing the spatial distribution of the probable errors, which can be used to indicate the likelihood, or probability, of any location falling above or below a specified elevation, thus this approach is called the "probabilistic" method. There is a long history of treating elevation error probabilistically (Hunter and Goodchild, 1995; Fisher, 1998; Zerger et al., 2002; Wechsler and Kroll, 2006), and the approach has been applied successfully in several recent SLR and flooding assessments

(Cooper and Chen, 2013; Leon et al., 2014; Cooper et al., 2015; Enwright et al., 2017; Fereshtehpour and Karamouz, 2018). In using the probabilistic approach, random error fields that match the error distribution characteristics derived from DEM accuracy assessment are generated and applied spatially. The assumption is that the elevation error is random, but elevation exhibits spatial autocorrelation, thus the error also has spatial autocorrelation (Wechsler and Kroll, 2006). Two techniques have been used to account for spatially autocorrelated errors in propagation through simulation: increasing the spatial autocorrelation in the random error fields by spatial filtering before they are applied to the DEM (Torio and Chmura, 2013; Enwright et al., 2017), and error modeling through sequential Gaussian simulation (Leon et al., 2014; Fereshtehpour and Karamouz, 2018). Other implementations of the probabilistic approach to handling vertical uncertainty in SLR assessments do not explicitly account for spatially autocorrelated elevation error (Cooper and Chen, 2013; Cooper et al., 2015; West et al., 2018), so they could be thought of as a type of worst-case scenario where the error is completely random without any spatial dependence on elevation. For users of assessments following the probabilistic approach, the stated probability indicates the likelihood, or chance, that the delineated area will be impacted or flooded at the specified water level, for example a 95% chance (or at least 95 times out of 100) that the area will be inundated.

The modified deterministic and probabilistic methods of accounting for vertical uncertainty are preferred as best practices over the simple deterministic method that ignores the effects of elevation error. With the preferred modified deterministic and probabilistic methods, maps of impact areas, and the associated inventory of population and resources located within impact zones, have increased value because of attached statements of confidence level or probability. The probabilistic approach is an excellent choice for a study that has ready access to sufficient compute resources required for error propagation through simulation, especially if the study area is large and the input data have high spatial resolution. Alternatively, the modified deterministic approach is suitable and is recommended. The results of the probabilistic method allow selection of different probabilities with which to present maps or statistics of areas exposed to inundation, while the modified deterministic method is less flexible in that it presents results within a bounded range at a specific confidence level. For the probabilistic method, if the full reference dataset used for accuracy assessment of the DEM is available, then a measure of spatial autocorrelation of the errors can be made and subsequent sequential Gaussian simulation can be performed. However, if only a global measure of the DEM vertical accuracy is available, for instance RMSE, then a good choice is to perform spatial filtering of the random error fields to account for autocorrelation as part of the simulation process for error propagation.

In recognition of the importance of considering vertical uncertainty in assessments of SLR and flooding exposure, there is a critical choice of parameters that must be made at the outset of a study: the increment of water level increase to be modeled, and the planning horizon (timeframe of projection into the future). The next three subsections (Minimum Sea-Level Rise Increment, Cumulative Vertical Uncertainty, and Minimum Planning Timeline) describe how the selection of these parameters needs to consider the vertical quality of the input elevation data (and associated datum transformations) and how the choices affect the reliability, or confidence level, of the assessment results.

#### Minimum Sea-Level Rise Increment

Any elevation-based SLR or coastal flooding assessment that uses a DEM, whether a simple bathtub model or a more complex hydraulic model is employed, raises the water level on a geospatial dataset that represents the topography of the study area. Such a process is essentially an elevation contouring process whereby a line of constant elevation (at the selected water level increase) is derived from the spatial arrangement of individual elevation values at discrete locations (usually in a regular grid in the case of a raster DEM). It is easy to define such an elevation contour, especially in a digital geographic information system (GIS), and the vertical increment between adjacent contours, referred to as the contour interval, is a parameter that must be specified in the procedure. A small interval can be applied to any DEM, but doing so does not imply that the derived contours automatically meet published accuracy standards. The interval must not be so small that it falls within the bounds of vertical error of the DEM, as such an operation would place the measurement (elevation increment) "in the noise" of the underlying elevation data. In the case of SLR or flooding assessments, the amount of water level increase from its current elevation to a future projected elevation is analogous to the contour interval, and that increment of increase must be larger than the inherent vertical error of the DEM for the projected future level to have a high level of confidence.

Based on the concept of elevation contour line accuracy, a method has been developed (Gesch, 2012a, 2013) to determine the minimum contour interval, or in the present case, the minimum increment of water level increase that can be used to meet a specified confidence level. Using the minimum increment in an assessment ensures that the chosen study parameter (amount of SLR or flooding level) is truly supported by the DEM and is not too small given the inherent vertical uncertainty. In the U.S., legacy national map accuracy standards applied to topographic contour maps specify that 90% of tested elevations should fall within one-half of the map contour interval (Maune et al., 2007), and this has been called the vertical map accuracy standard (VMAS) with a 90% confidence level, or alternatively "linear error at 90% confidence" (LE90). From the map user perspective, elevations determined from the map on 9 out of 10 points will have true elevations that are within one-half of the contour interval (CI). The contour accuracy standard is expressed in the following equation:

$$\text{VMAS} = \text{LE90} = \frac{\text{CI}}{2}$$

Rearranging this simple equation as

$$\text{CI} = \text{LE90} \times 2$$

allows the contour interval to be expressed as a factor of the elevation data accuracy. Additional accuracy standards developed more recently for digital geospatial data rather than

hardcopy maps, such as the National Standard for Spatial Data Accuracy (Maune et al., 2007) and the American Society for Photogrammetry and Remote Sensing (ASPRS) Positional Accuracy Standards for Digital Geospatial Data (ASPRS, 2015), provide procedures for direct conversion among DEM accuracy metrics and equivalent contour intervals. For instance, much of the coastal and floodplain light detection and ranging (lidar) data in the U.S. was collected to a product specification by the Federal Emergency Management Agency (FEMA) for DEMs with an RMSE of 0.185 m, which was calculated to provide elevation data that would support topographic mapping at a 2-foot (0.61 m) contour interval accuracy (Maune et al., 2007).

In the present case of SLR or coastal flooding, the increment of water level increase is analogous to the contour interval because raising the water level on an elevation dataset is equivalent to a contouring operation, especially when multiple successive water levels are mapped. Thus, the minimum water level increment for modeling can be stated directly as a factor of the elevation data accuracy, expressed as a common vertical error metric (LE90 in the preceding example). Two of the most commonly used DEM error metrics are RMSE and LE95, and direct translations among RMSE, LE90, and LE95 are available (Maune et al., 2007), assuming the errors are from an unbiased normal distribution. Because the error metrics represent a portion of the cumulative probability distribution of errors, a confidence level can be stated for the minimum increment, for example 68% confidence for RMSE (equivalent to the "one sigma" error, or standard deviation of the errors for an unbiased normal distribution), 90% confidence for LE90, and 95% confidence for LE95. Zhang et al. (2011) recognize that the modeled increments of SLR should be tied to the inherent vertical error, and they did so by selecting increments that matched the RMSE of their input DEM, whereas the approach described here has the added advantage of providing a direct method to calculate the proper increments as a function of the DEM accuracy at a specific confidence level. This approach, rooted in contour interval accuracy standards, provides the quantitative basis for the "guideline" (Gesch et al., 2009) and "rule of thumb" (National Oceanic and Atmospheric Administration [NOAA], 2010) that the increment of SLR modeled should be at least twice the vertical accuracy of the elevation data.

As an example to illustrate, consider a DEM with an RMSE of 0.15 m derived from airborne lidar data. Assuming the errors are unbiased and normally distributed, the RMSE can be converted to LE95 by the following formula (Maune et al., 2007):

$$\text{LE95} = \text{RMSE} \times 1.96$$

resulting in an LE95 of 0.294 m. Applying the procedure described above, the minimum SLR increment, abbreviated here and referred to hereafter as SLRImin, is 0.588 m at the 95% confidence level, and 0.30 m at the 68% confidence level. In equation form, SLRImin at 95% confidence is

$$\text{SLR}I\_{\text{min95}} = (\text{RMSE} \times 1.96) \times 2$$

and SLRImin at 68% confidence is

$$SLRI\_{min68} = \text{RMSE} \times 20$$

The SLRImin metric can be interpreted as follows: it expresses the confidence in how well the "contour" line delineating areas with elevations at or below the raised water level is placed vertically. To follow through with the example illustration, there is a 68% chance that a DEM-derived line delineating the inundated area will be placed vertically within ±0.15 m of where the true line is, and likewise, there is a 95% chance that the line will be placed within ±0.294 m vertically of its true location. Because the SLRImin is a direct function of the DEM accuracy, DEMs with lower accuracy would only support much larger water level increments. For instance, a DEM with an RMSE of 5 m would only support a SLR or flooding level increment of 10 m at the 68% confidence level, and using a smaller increment would have a drastically reduced confidence level, progressing to the point where using increments of less than 1 m would result in confidences near 0%.

#### **Cumulative vertical uncertainty**

Digital elevation model error is the main source of vertical uncertainty in elevation-based assessments, but there are other contributors of error, namely the datum transformations required to bring the DEM into a tidal datum reference framework. It is important to include local water level information when mapping potential impacts (Marbaix and Nicholls, 2007) by starting at the high tide line, as the areas below this line are already subject to periodic submersion from the normal range of tides. Mapping impact areas upslope of the normal high water line is recommended here as a best practice. Many DEMs are referenced to an orthometric (mean sea level referenced) datum, and thus require transformation to a tidal datum, often mean higher high water (MHHW), before analysis, and these vertical transformations add more vertical uncertainty. In the U.S., a tool called VDatum (Parker et al., 2003) is widely used for such processing, and it has published uncertainties for the various transformations<sup>1</sup> . Several SLR assessment studies have combined the DEM error and vertical datum transformation error with a root sum of squares (or summing in quadrature) approach to calculate the cumulative vertical uncertainty (Mitsova et al., 2012; Cooper et al., 2013, 2015; Gesch, 2013; Schmid et al., 2014; Enwright et al., 2015), and such a procedure is recognized here as a best practice for elevation-based SLR assessments. The preceding section (Minimum Sea-Level Rise Increment) on minimum SLR increment used the DEM vertical accuracy to calculate the critical assessment parameter, but in practice the cumulative vertical uncertainty, if known, should be used, as demonstrated in Gesch (2013).

#### Minimum Planning Timeline

Another critical assessment parameter is the planning horizon, or the timeframe over which projected increased water levels are mapped to delineate potential impact zones. Like SLRImin, the minimum planning timeline, designated here and referred to hereafter as TLmin, is directly related to the vertical uncertainty of the input DEM, and it incorporates the rate of SLR projected over the time scale of interest. For illustration, assume a linear rate of

<sup>1</sup>https://vdatum.noaa.gov/docs/est\_uncertainties.html

SLR, and then TLmin can be calculated as

feart-06-00230 December 11, 2018 Time: 16:29 # 6

$$TL\_{min} = \text{SLRI}\_{min} \div \text{annual SLR rate}$$

For example, consider the minimum and maximum from the range of likely global SLR scenarios from the Fifth Assessment Report (AR5) by the Intergovernmental Panel on Climate Change (IPCC), 0.28 – 0.98 m by the year 2100 (Church et al., 2013). To simplify the illustration, assume the numbers represent the cumulative SLR over 100 years (2000–2100), so the annual increment for the minimum and maximum are 2.8 mm/year and 9.8 mm/year, respectively. Given a lidar-derived DEM with an RMSE of 0.15 m, TLmin for the minimum scenario is 107 years, and TLmin for the maximum scenario is 31 years. For the minimum scenario, mapping a potential impact zone for any year before 2107 would be unreliable as the cumulative water level increase will not have reached the minimum SLR increment afforded by the elevation data at the specified confidence level. For the maximum scenario, delineating a potential inundation zone would be acceptable for any timeframe after 2031. This simple example uses linear SLR rates, but TLmin can be based on non-linear scenarios as well.

As with SLRImin, a confidence level is associated with TLmin because the input DEM error metric carries a confidence level with it, in this case 68% confidence as RMSE was used, so the result can be noted as TLmin68. Likewise, if SLRImin<sup>95</sup> is used to calculate the minimum planning timeline, the result is noted as TLmin<sup>95</sup> and the confidence level is 95%.

Overall, SLRImin and TLmin are useful to determine what parameters can be effectively used in assessments, especially from a management perspective (Gesch, 2012a). When a specific DEM with a stated accuracy is available, SLRImin and TLmin will be useful for determining what increments and planning horizons (given a SLR rate) will be allowable for high confidence results. Alternatively, if specific targets for modeled increments and planning horizons are known, along with SLR scenarios, then the quality of elevation data (that is, its accuracy) to meet specific confidence levels can be determined.

### Digital Elevation Models

There are numerous global or near-global DEMs available and they have all been used, some extensively, for SLR and coastal flooding assessments. The availability of global DEMs has improved over the last several years and new or refined products continue to appear, so it is important to understand the vertical uncertainty of these DEMs and how that affects their effective use in applications. For this study, DEMs with a medium to high spatial resolution (better than 100 m) and an open data distribution policy are included. Other lower resolution topographic products (250 m to several kilometers) are available but are excluded here as most are derived from the higher resolution global DEMs. The following DEMs are included for analysis:

(1) Shuttle Radar Topography Mission (SRTM) (Farr et al., 2007) data are available for all land areas between 60◦ north and 56◦ south latitude at 1-arc-second (30-m) and 3-arc-second (90-m) grid spacing. The SRTM product specification for vertical accuracy is 16 m LE90, which equates to an RMSE of 9.73 m.


Also included in some of the comparisons are high-resolution, high-accuracy DEMs derived from airborne lidar and stereo imagery to provide context for the global DEM results. These DEMs include the U.S. Geological Survey (USGS) National Elevation Dataset (NED) (Gesch et al., 2002; Gesch, 2007) and 3D Elevation Program (3DEP) (Sugarbaker et al., 2014) lidar DEMs (Heidemann, 2012) for the U.S. An unmanned aerial system (UAS) derived DEM generated with structure from motion (SfM) techniques for Majuro Atoll in the central Pacific island nation of Republic of the Marshall Islands (RMI) (Palaseanu-Lovejoy et al., 2018) is included as an example of newer technologies that are increasingly being used to generate high-resolution, high-accuracy elevation data at local scales.

#### Accuracy Assessment

To obtain the required DEM accuracy metrics that characterize vertical uncertainty, an accuracy assessment was conducted for each of the DEMs. In this case, the value of interest for each DEM is the absolute vertical accuracy, which is calculated from the statistics of a set of reference (or truth) points compared to the DEM. In each case, the elevation value at every reference point is compared to the corresponding DEM elevation (extracted via bilinear interpolation at the exact point location) and the difference in elevations is recorded. The difference represents

the DEM error at that point. The differencing operation is done by subtracting the reference point elevation from the DEM elevation. In this manner, the difference statistics from the full set of point comparisons are easy to interpret; that is, a positive mean error indicates that on average the DEM is too high (the DEM has a positive bias). Conversely, a negative mean error indicates that on average the DEM is too low (a negative bias). This differencing approach is recommended as a best practice for absolute vertical accuracy assessments that use ground truth point data to assess raster DEM datasets.

Prior to comparison of the DEM and the reference data, both datasets must be in the same vertical reference frame so the difference statistics do not contain any artificial biases. The DEMs and reference data used here are a mix of different vertical reference systems: SRTM, ASTER GDEM, AW3D30, NASADEM, and MERIT DEM are referenced to the Earth Gravitational Model 1996 (EGM96) geoid; TanDEM-X is referenced to the World Geodetic System 1984 (WGS84)-G1150 ellipsoid; NED and the U.S. ground truth points are referenced to the North American Vertical Datum of 1988 (NAVD 88) orthometric datum; the Majuro UAS-SfM DEM and ground truth points are both referenced to the International Terrestrial Reference Frame 2008 (ITRF2008) ellipsoid. The ground truth points and each DEM were brought into the same vertical reference frame with a procedure similar to that described in Grohmann (2018), and the VDatum software tool was used in the process.

#### **Elevation reference data**

The reference data for the conterminous United States (CONUS) is an extensive set of high-accuracy geodetic control points produced by the U.S. National Geodetic Survey (NGS) and known as "GPS on Bench Marks"<sup>2</sup> . These points are NGS's best control points, with millimeter- to centimeter-level accuracies, so they are an excellent reference dataset for comparing with DEMs for accuracy assessment purposes. The points have been used extensively for such analyses (Gesch, 2007; Gesch et al., 2014, 2016; Wessel et al., 2018). For the Majuro DEM, the reference data are an extensive set of real-time kinematic (RTK) GPS points collected with survey-grade equipment during the UAS flights (Palaseanu-Lovejoy et al., 2018).

#### **Test areas**

The DEM accuracy assessments and comparisons were conducted across CONUS and in selected coastal locations. Because not all DEMs were available for the full U.S. coastal zone, a set of 17 one-degree by one-degree tiles served as a subset test area where direct comparisons of all DEMs could be made. **Figure 1** shows the GPS on Bench Marks reference data for all of CONUS (23,115 points), the subset of points (3,480) in the low elevation coastal zone (LECZ; defined here as areas less than or equal to 10 m in elevation), and the 17 test tiles. The number of test tiles is due to a limited amount of data available from the TanDEM-X data provider, and these locations are places where USGS scientists have ongoing coastal DEM development and applications activities. Even though the analysis areas are limited

to the U.S. coastal zone, the coverage is extensive enough that the results are deemed applicable for guiding the use of global DEMs in other areas.

#### RESULTS

The initial accuracy assessment was conducted on the DEMs for which coverage was available for all of CONUS: SRTM, ASTER GDEM (version 3), NED, and NED derived from lidar DEM source data (Gesch et al., 2014). **Table 1** shows the accuracy assessment results for the full range of elevations across CONUS included in the reference data, as well as the results for only areas in the LECZ. The LECZ is defined as elevations less than or equal to 10 m above sea level, which is a commonly used elevation threshold to delimit coastal zones (McGranahan et al., 2007; Lichter et al., 2011; Neumann et al., 2015). For SRTM and ASTER GDEM, the elevation accuracy in the LECZ is degraded compared to the accuracy for all

<sup>2</sup>https://www.ngs.noaa.gov/GEOID/GEOID12B/GPSonBM12B.shtml

TABLE 1 | Accuracy assessment results for DEMs over CONUS.


of CONUS, while for NED and NED with lidar source data, the accuracy is slightly improved in the LECZ. Because the primary interest for coastal assessments is in the low-lying areas subject to inundation and other adverse effects of increased water levels, the remainder of the results presented are for the LECZ. This approach follows the recommendation in other studies that emphasize accuracy testing with reference data representative of the area of interest to guard against overly optimistic results (Bolkas et al., 2016). Thus, limiting accuracy assessment to coastal areas (Du et al., 2015) is appropriate for this study.

# Vertical Accuracy and Inundation Assessment Parameters

**Table 2** displays the accuracy information for DEMs over the CONUS LECZ. In addition to the four DEMs in **Table 1**, the specification for 3DEP lidar is added for comparison, as it represents high-resolution, high-accuracy elevation data being collected over broad areas in the U.S. 3DEP is an ongoing program of the U.S. Government to coordinate and collect enhanced elevation data for CONUS in an 8-year cycle (Sugarbaker et al., 2014). Most of the data are being collected to a specification for "quality level 2" (QL2), which requires a vertical RMSE of 0.10 m (Heidemann, 2012). The DEM accuracy metrics included in the table are RMSE, mean error, and mean absolute error (MAE). As noted above, mean error can be indicative of overall positive or negative bias in a DEM and, if a bias is present, can reflect a departure from a normal distribution of the errors (Maune et al., 2007). In these cases, the MAE can be a useful metric to help describe the error characteristics (Chai and Draxler, 2014). The positive mean errors for SRTM and ASTER GDEM do indicate that on average the DEMs are too high relative to the ground (as represented in the reference data point elevations) and the error distribution is not an unbiased normal distribution. The performance of SRTM and ASTER GDEM generally having a positive bias has been noted in previous published accuracy evaluations of these DEMs (Carabajal and Harding, 2006; Gesch et al., 2016).

In these cases of a biased error distribution, an alternative error metric to the RMSE (or its calculated equivalents like LE95) is a sample quantile of the cumulative error distribution (Höhle and Höhle, 2009), which has been widely implemented and used as the "95th percentile" error approach for describing the vertical accuracy (at the 95% confidence level) of DEMs with non-normal error distributions (Maune et al., 2007; ASPRS, 2015). Other sample quantiles can be used, for instance 68


or 90% (Wessel et al., 2018), for the percentile error metric. **Table 2** includes both the 68th percentile and 95th percentile errors for each of the DEMs tested, along with the derived SLRImin measure, calculated from both the RMSE and percentile error metrics. For the percentile error metrics, the corresponding SLRImin measure is simply two times the percentile error. If a DEM accuracy validation is done for a specific area, then there is more flexibility in using the various error metrics to derive SLR assessment parameters (minimum increment and planning timeline), with the ability to deal with factors such as a biased non-normal error distribution with large outliers. In practice, most coastal assessments do not include DEM accuracy testing and characterization, but instead must rely upon DEMs with published accuracy figures, and RMSE is the only metric available with which to determine the proper increment and time horizon parameters. For this reason, comparisons of SLRImin and TLmin presented and discussed below are derived from the RMSE for the DEMs analyzed in this study.

**Table 3** presents the results for the direct comparison of all the tested global DEMs within the 17 coastal test tile locations (**Figure 1C**). All the global DEMs exhibit accuracies that are better than their product specifications, as is often the case when these DEMs have been evaluated over broad areas. Refer to the tables in the **Supplementary Material** for a record of the numerous accuracy tests for SRTM (and its derivatives) (**Supplementary Table S1**), ASTER GDEM (**Supplementary Table S2**), AW3D30 (**Supplementary Table S3**), and TanDEM-X (**Supplementary Table S4**). Also included in the **Supplementary Material** are a summary of evaluations of NED (**Supplementary Table S5**), as an example of higher resolution DEMs with regional to national coverage, and a summary of example accuracies from technologies that produce very high-resolution and very high-accuracy elevation data that have been used in coastal assessments (**Supplementary Table S6**). The results in these tables provide context to the capabilities of global DEMs for coastal assessments and demonstrate the possibilities for very high-accuracy mapping when requirements call for detailed spatially explicit information.

**Table 4** and **Figure 2** show the evaluated DEMs ranked in order of increasing vertical uncertainty. The derived parameters of SLRImin and TLmin needed for coastal assessments have been calculated from the vertical accuracy for each DEM (as stated in the RMSE) and are presented at the 68 and 95% confidence levels. Included in the list of DEMs is terrestrial lidar, which is another example of a high-accuracy source of elevation data that can be used in detailed assessments. In this case, the vertical accuracy (RMSE = 0.05 m) comes from an example coastal DEM used for a flooding assessment on Kwajalein Atoll, RMI (Storlazzi, 2017). The results in **Table 4** and **Figure 2** indicate that only the high-resolution, high-accuracy DEM sources allow a SLRImin of less than 1 m and a TLmin of less than 100 years at high confidence levels, while the national scale or global DEMs do not adequately support such parameters.

**Figure 3** shows how for a given DEM (and its associated vertical accuracy) the confidence level attached to SLRImin



TABLE 4 | Ranking of elevation data sources (in order of increasing RMSE) and corresponding SLRImin and TLmin for use in coastal assessments.

FIGURE 2 | Ranking of elevation data sources (left to right in order of increasing SLRImin). The 1-m sea-level rise (or inundation level) is marked, which indicates that only the high-accuracy DEM sources are suitable for modeling sub-meter increments at high confidence levels.

increases as the increment used for modeling increases. Note that when global DEMs and sub-meter increments of SLR (0.5 m or less) are used, the confidence is very low (in the 0–10% range). A simple formula for calculating the confidence level for SLRImin given the RMSE of the DEM and the modeled increment is described in the **Supplementary Material**.

is included because it is the equivalent of 1 foot, an increment that is commonly used in local and regional inundation assessments in the U.S.

TABLE 5 | The confidence level for delineating the 10-m LECZ and the 20-m coastal zone (CZ) using each global DEM.

feart-06-00230 December 11, 2018 Time: 16:29 # 12


# Delineation of Low Elevation Coastal Zone

In addition to developing spatially explicit flooding or inundation impact zone maps, DEMs are also used to outline general LECZ delineations (McGranahan et al., 2007; Lichter et al., 2011; Neumann et al., 2015). Given the accuracy of global DEMs, and the corresponding SLRImin for each, the confidence levels for delineating the 10 m or less LECZ (that is, the 10-m elevation contour) and the 20 m or less coastal zone (CZ) (the 20-m contour) can be calculated (**Table 5**). The global DEMs suitable for delineating the LECZ (≤10 m elevation) at 68% confidence are TanDEM-X, CoastalDEM (an SRTM derivative) (Kulp and Strauss, 2018), NASADEM, AW3D30, and MERIT. At the 95% confidence level, only TanDEM-X is suitable for delineating the LECZ. All the global DEMs are acceptable for delineating the CZ (≤20 m elevation) at 68% confidence, while at the 95% confidence level SRTM and ASTER GDEM are no longer suitable, but the others are.

# DISCUSSION

When the inherent vertical uncertainty is considered, it is clear that global DEMs are not adequate for modeling fine increments of SLR (<1 m) over short planning horizons (<100 years) at high confidence levels (**Table 4**). SLRImin and TLmin are easily calculated metrics based on the stated DEM accuracy (usually expressed as an RMSE) and are useful to quantify elevation-based coastal assessment parameters and their associated confidence levels. Numerous studies have used global DEMs for SLR and coastal flooding assessments without any regard for vertical uncertainty (elevation error), leading to large uncertainties in the assessment results with low confidence levels. There is ample evidence that SRTM and ASTER GDEM are severely limited for use in coastal assessments (Gesch, 2009; Gesch et al., 2009; van de Sande et al., 2012; Doyle et al., 2015; Griffin et al., 2015; Yan et al., 2015; Kulp and Strauss, 2016; Walczak et al., 2016; Yunus et al., 2016; Santillan and Makinano-Santillan, 2017; Smith et al., 2018). Both SRTM and ASTER GDEM are DSMs that generally overestimate elevations (especially in vegetated and built-up areas), so their use in coastal assessments leads to underestimating areas exposed to a given inundation level (van de Sande et al., 2012; Griffin et al., 2015; Kulp and Strauss, 2016; Smith et al., 2018). The results from the DEM accuracy assessment in this study (**Table 3**) do indicate the overestimation of elevation (positive bias) by SRTM and ASTER GDEM as reflected in the positive mean error for each. A bias (that is, a larger negative or positive mean error) is often an indicator that the error distribution is non-normal; therefore, applying linear scaling factors to the RMSE to derive LE90 or LE95 is not valid, so percentile error methods should be used (Maune et al., 2007; ASPRS, 2015). In the analysis presented here, the RMSE was used to derive SLRImin even when the DEMs had a positive bias, as often only the RMSE is known and the individual error values are not available for percentile error calculations. As shown in **Table 3**, however, SLRImin<sup>68</sup> and SLRImin<sup>95</sup> calculated based on the 68th percentile and 95th percentile errors, respectively, are not much different than those calculated based on the RMSE, so the finding that none of the global DEMs support an increment of 1 m is not changed.

There have been some recent improvements to SRTM (O'Loughlin et al., 2016; Kulp and Strauss, 2018; Moudrý et al., 2018) and ASTER GDEM (Arefi and Reinartz, 2011; Yang et al., 2018) by removing vegetation and other elevated features that caused the positive bias in the original datasets, often with the correction implemented by integrating Ice, Cloud and land Elevation Satellite (ICESat) spaceborne lidar data,. Merges of SRTM and ASTER GDEM (Satgé et al., 2015; Crippen et al., 2016; Yamazaki et al., 2017) have resulted in improved data as well. However, none of these improvements bring the DEMs to the level where they will support high confidence, quantitative coastal assessments with sub-meter water level change increments and planning horizons within the current century. The improvements to SRTM, namely CoastalDEM, NASADEM, and MERIT, all still have an RMSE of about 3 m, which equates to a minimum increment of more than 6 m at 68% confidence (**Table 4** and **Supplementary Table S1**). AW3D30 is in this same class, with an RMSE of slightly more than 3 m (**Table 4** and **Supplementary Table S3**). TanDEM-X does offer an improvement over SRTM, ASTER GDEM, and AW3D30 (Grohmann, 2018) and exhibits very little positive bias (Wessel et al., 2018), although its vertical error does result in a SLRImin of several meters (**Table 4** and **Supplementary Table S4**).

Despite ample evidence of the significant limitations of global DEMs, especially SRTM, for coastal assessments, they have been used extensively for mapping and describing potential impacts of SLR and coastal flooding, often with assessment parameters (small water level increments and short planning horizons) that fall well within the error bounds of the underlying elevation data (Dasgupta et al., 2008, 2010; Hanson et al., 2010; Curtis and Schneider, 2011; Blankespoor et al., 2014; Hardy and Nuse, 2016; Kopp et al., 2017; Runting et al., 2017; Brown et al., 2018a,b; Gebremichael et al., 2018; Haer et al., 2018; Jevrejeva et al., 2018; Lincke and Hinkel, 2018; Nicholls et al., 2018; Prahl et al., 2018; Rasmussen et al., 2018; Schuerch et al., 2018; Wolff et al., 2018). Some of these studies used a model or database in which the global DEM is embedded, such as the Dynamic Interactive Vulnerability Assessment (DIVA) modeling framework (Hinkel, 2005; Vafeidis et al., 2008), so the inherent vertical uncertainty is contained within model or database components. This serves as a

caution that even if a DEM is not a direct input or is not processed or analyzed directly in an inundation modeling exercise, vertical error may be implicit in sub-model or database components, thus modelers should be aware of an obscured source of uncertainty in their assessment. Some assessments have used even coarser global elevation models (1-km spatial resolution) with water level increments in the range of 0.3–2 m (Xingong et al., 2009; Nicholls et al., 2011; Brown et al., 2013; Neumann et al., 2015). All these assessments present statistics about the potential impact zones, including the areas and oftentimes the corresponding population and economic assets that are at risk of adverse effects. However, there is no quantitative statement about the quality of the reported results, for instance confidence level, and no expression of the uncertainty contributed by vertical error. In all cases, the water level change increments used in the assessments are not supported at high confidence levels, as the increments are very small compared to the DEM vertical error (and the derived minimum increment), which calls into question the veracity of the reported results. Reporting of inundated areas (and population and resources contained therein) at intervals of 1 m or less implies a degree of accuracy that is not present in global DEMs, and providing such numbers can be misleading to readers, especially because most will not be familiar with the concepts of vertical uncertainty of elevation models. As much as there is a desire and need for global SLR and coastal flooding analyses with water level increases on the order of a meter or less, current global DEMs do not have the requisite vertical accuracy to derive results with high confidence levels using fine increments, and thus they should not be used for such mapping.

Even though global DEMs are not appropriate for spatially explicit mapping of small increments of inundation with high confidence, they can be used effectively for delineation of general LECZs, and inventorying the population and resources contained within. As the RMSE improves for global DEMs, the confidence level for delineation of 10- and 20-m coastal zones increases (**Table 5**). The LECZ can be the framework for entire studies (McGranahan et al., 2007; Geisler and Currens, 2017), so a quality delineation of such a zone at a known confidence level is critical. For finer vertical slices, and subsequent spatially explicit exposure maps, much higher accuracy elevation data, such as that derived from lidar, high-resolution stereo photogrammetry, and ground survey, are required (Gallien et al., 2013). Airborne lidar, in particular, is an important elevation data source for coastal assessments (Gesch, 2009; Cooper et al., 2013; Runting et al., 2013; Zhu et al., 2015; Enwright et al., 2017), as it can meet the requirements for modeling fine increments of water level changes at high confidence and generally covers larger areas, even regional to national coverage.

# Proper Accounting for Vertical Uncertainty

When vertical uncertainty is properly accounted for, value is added to the results of coastal assessments as additional information is available to inform users. For example, this can take the form of confidence levels attached to the inventory of resources within a potential impact zone or as portrayal of the probability of inundation for a specific location or confidence in the delineation of flooding exposure on a map. Several studies demonstrate best practices for handling cumulative vertical uncertainty in both the selection of assessment parameters (modeled increments and projection timelines) and in presentation of results (graphic and tabular) with expressions of confidence. For implementations of these best practices, see the following examples: Reynolds et al. (2012); Gesch (2013); Nielsen and Dudley (2013); Leon et al. (2014); Enwright et al. (2015, 2017); Dahl et al. (2017); Jones et al. (2017); Santillan and Makinano-Santillan (2017); West et al. (2018).

# Other Inundation Exposure Assessment Best Practices

In addition to rigorous treatment of vertical uncertainty (detailed above in Section "Accounting for Uncertainty in Exposure Assessments" – spatial portrayal of cumulative vertical uncertainty when mapping inundation zones, and in Sections "Minimum Sea-Level Rise Increment" and "Minimum Planning Timeline" – selection of assessment parameters), there are other best practices that will help produce high quality coastal assessments. The following practices have emerged from the scientific record of successful studies, and they are becoming commonplace in the most robust assessments.


population within the potential impact zone. Dasymetric population mapping (Mennis, 2003; Holt et al., 2004) is an effective technique for disaggregating areal population counts to a more realistic distribution of population density across the landscape as a continuous surface, often using land cover or parcel data as ancillary information. The advantages of doing so for coastal assessments have been demonstrated (Mitsova et al., 2012; Merkens and Vafeidis, 2018), and spatially explicit population maps are available over large areas (Mondal and Tatem, 2012; Dmowska and Stepinski, 2017).

# CONCLUSION

Vertical uncertainty is a critical factor to consider and account for in elevation-based assessments of SLR and coastal flooding exposure. Some studies have properly handled the vertical uncertainty, as expressed in the combined elevation data and transformation errors, and in doing so provide valuable additional information to the user about confidence in the mapping and likelihood of projected impacts. However, many other studies ignore the vertical uncertainty stemming from the underlying elevation data and use assessment parameters (water level change increments and planning horizons) that are well within the error bounds and are not appropriate for generating high confidence results, thus leading to questionable delineations of impact zones and inventories of the population and resources contained therein.

The simple methods described herein for selecting coastal assessment parameters (minimum increment of water level change, SLRImin, and planning horizon, TLmin) that are supported at high confidence levels by the vertical qualities of the elevation data are useful for characterizing the capabilities of global DEMs. Application of these methods to current global DEMs (SRTM and its derivatives NASADEM, CoastalDEM, and MERIT; ASTER GDEM; AW3D30; and TanDEM-X) demonstrates that none of these DEMs support coastal inundation or flood assessment at high confidence levels for small water level increments (<1 m) or short planning horizons (<100 years). High confidence assessments of scenarios with cumulative SLR of less than 1 m or planning horizons within the current century require elevation data with much better vertical accuracy than that afforded by global DEMs, which points to high-accuracy sources such as terrestrial and airborne lidar, high-resolution photogrammetry, and ground surveys. These technologies produce high-quality elevation data that facilitate development of detailed spatially explicit inundation maps.

The key finding demonstrated in this study leads the list of best practices to follow in elevation-based coastal inundation assessments.

(1) Account for the inherent cumulative vertical uncertainty in the elevation data by using increments of water level increase and planning horizons that are supported at high confidence levels, and state those confidence levels explicitly in study documentation. The metrics SLRImin and TLmin are direct functions of the vertical accuracy of the DEM used in the study, and they are useful for ensuring that the chosen assessment parameters are appropriate given the error characteristics of the DEM.


As the use of these best practices increases, assessments will improve and become more valuable, especially by having quantified and published uncertainty information (confidence levels and likelihood statements), and results will be directly comparable across different assessments.

In the future, as elevation datasets with large-area coverage improve, analyses utilizing the improved elevation information and the community best practices will result in robust assessments. Ongoing enhancements to widely used methods will also help to improve progress, such as better incorporation of information on physical processes, including tidal regimes (Hanslow et al., 2018) and water level attenuation due to surface roughness (Vafeidis et al., 2017), into bathtub modeling used for broad area screening. However, for elevation-based coastal assessments, the primary factor affecting quality and usefulness of results remains the choice of the elevation model used (National Oceanic and Atmospheric Administration [NOAA], 2010; Doyle et al., 2015; Wolff et al., 2016; Yunus et al., 2016), and how the DEM vertical uncertainty is characterized and accounted for (West et al., 2018). Open-access global DEMs have been a major advance for many Earth science and environmental modeling applications, but the findings from the present evaluation of currently available datasets for detailed assessments of SLR and coastal flooding exposure add to the recent recognition (Schumann et al., 2014; Simpson et al., 2015; Sampson et al., 2016) that the requirement remains for a freely available, high-accuracy, high-resolution global elevation model that supports quantitative coastal inundation hazard assessments at high confidence levels.

# DATA AVAILABILITY STATEMENT

All data generated or analyzed during this study are included in the main text of this publication.

# AUTHOR CONTRIBUTIONS

feart-06-00230 December 11, 2018 Time: 16:29 # 15

DG designed the study, performed the data collection, processing, and analysis, and wrote the manuscript.

# FUNDING

Funding for this research was provided by the Land Change Science Program of the U.S. Geological Survey. TanDEM-X data (©DLR 2017) were provided by the German Aerospace Center (DLR) through the TanDEM-X Digital

# REFERENCES


Elevation Model Announcement of Opportunity (Proposal ID: DEM\_HYDR1176).

## ACKNOWLEDGMENTS

The author thanks Monica Palaseanu-Lovejoy for valuable comments on an earlier version of the manuscript. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart. 2018.00230/full#supplementary-material


Maune (Bethesda, MD: American Society for Photogrammetry and Remote Sensing), 99–118.


rise. Proc. Natl. Acad. Sci. U.S.A. 111, 3292–3297. doi: 10.1073/pnas.122246 9111


"Low-Elevation Coastal Zone" (LECZ). J. Coast. Res. 274, 757–768. doi: 10.2112/ jcoastres-d-10-00072.1




**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Gesch. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Perspectives on Digital Elevation Model (DEM) Simulation for Flood Modeling in the Absence of a High-Accuracy Open Access Global DEM

#### Laurence Hawker<sup>1</sup> \*, Paul Bates<sup>1</sup> , Jeffrey Neal<sup>1</sup> and Jonathan Rougier<sup>2</sup>

<sup>1</sup> School of Geographical Sciences, University of Bristol, Bristol, United Kingdom, <sup>2</sup> School of Mathematics, University of Bristol, Bristol, United Kingdom

#### Edited by:

Matthew McCabe, King Abdullah University of Science and Technology, Saudi Arabia

#### Reviewed by:

Salvatore Manfreda, University of Basilicata, Italy Ben Jarihani, University of the Sunshine Coast, Australia

\*Correspondence: Laurence Hawker Laurence.hawker@bristol.ac.uk

#### Specialty section:

This article was submitted to Hydrosphere, a section of the journal Frontiers in Earth Science

Received: 25 July 2018 Accepted: 30 November 2018 Published: 18 December 2018

#### Citation:

Hawker L, Bates P, Neal J and Rougier J (2018) Perspectives on Digital Elevation Model (DEM) Simulation for Flood Modeling in the Absence of a High-Accuracy Open Access Global DEM. Front. Earth Sci. 6:233. doi: 10.3389/feart.2018.00233 Open-access global Digital Elevation Models (DEM) have been crucial in enabling flood studies in data-sparse areas. Poor resolution (>30 m), significant vertical errors and the fact that these DEMs are over a decade old continue to hamper our ability to accurately estimate flood hazard. The limited availability of high-accuracy DEMs dictate that dated open-access global DEMs are still used extensively in flood models, particularly in datasparse areas. Nevertheless, high-accuracy DEMs have been found to give better flood estimations, and thus can be considered a 'must-have' for any flood model. A highaccuracy open-access global DEM is not imminent, meaning that editing or stochastic simulation of existing DEM data will remain the primary means of improving flood simulation. This article provides an overview of errors in some of the most widely used DEM data sets, along with the current advances in reducing them via the creation of new DEMs, editing DEMs and stochastic simulation of DEMs. We focus on a geostatistical approach to stochastically simulate floodplain DEMs from several open-access global DEMs based on the spatial error structure. This DEM simulation approach enables an ensemble of plausible DEMs to be created, thus avoiding the spurious precision of using a single DEM and enabling the generation of probabilistic flood maps. Despite this encouraging step, an imprecise and outdated global DEM is still being used to simulate elevation. To fundamentally improve flood estimations, particularly in rapidly changing developing regions, a high-accuracy open-access global DEM is urgently needed, which in turn can be used in DEM simulation.

Keywords: digital elevation models, open-access, geostatistics, flood, stochastic simulation, floodplains, hazards

# INTRODUCTION

Digital Elevation Models (DEM) are a gridded digital representation of terrain, with each pixel value corresponding to a height above a datum. Since the pioneering work of Miller and Laflamme (1958), DEMs have grown to become an integral part of a number of scientific applications. DEMs can be created from ground surveys, digitizing existing hardcopy topographic maps or by remote sensing techniques. DEM's are now predominantly created using remote sensing techniques with

Smith and Clark (2005) observing the benefits that a large spatial area can be mapped by fewer people at a lower cost. Remotely sensing techniques include photogrammetry (Uysal et al., 2015; Coveney and Roberts, 2017), airborne and spaceborne Interferometric Synthetic Aperture Radar (InSAR) and Light Detection And Ranging (LiDAR). Spaceborne InSAR is the most common technique to create global DEMs and is the technology behind the most widely used open-access global DEM; the Shuttle Radar Topography Mission (SRTM). An overview of free and commercial global DEMs is given in **Figure 1**. Despite being acquired in 2000, SRTM is still the most popular global DEM because of its accessibility, feature resolution, vertical accuracy and a lower amount of artifacts and noise compared to alternative global DEMs (Rexer and Hirt, 2014; Jarihani et al., 2015; Sampson et al., 2016; Hu et al., 2017). Yet, recently released products such as MERIT and TanDEM-X 90 could change that. Commercial DEMs (Planet Observer, 2017; Takaku and Tadono, 2017; InterMap, 2018) are often prohibitively costly and have a lack of comparison studies. Airborne LiDAR has a higher precision and accuracy owing to its ability to penetrate vegetation and its reduced vulnerability to scatter but is largely limited to a handful of countries and can be expensive to acquire. These characteristics are conducive in creating a high-quality (vertical error <1 m) 'bare-earth DEM,' where objects (e.g., buildings and vegetation) have been removed from the elevation model. Such bare-earth DEMs are essential for applications, such as flood modeling, that rely on the accurate derivation of surface characteristics (e.g., slope).

# Characteristics of DEM Errors

DEM errors occur in both the horizontal and vertical directions. Errors propagate from the input data used in creating a DEM right through to calculating surface derivatives and using DEMs in complex applications (Hutchinson and Gallant, 2000; Fisher and Tate, 2006). Whatever their source, DEMs can appear to provide a definitive and plausible representation of topography which can often lull the user into a false sense of security regarding their accuracy, with many users unaware of DEM errors or how to treat them (Wechsler, 2003, 2007). Whilst accuracy statistics such as RMSE provide an indication of DEM accuracy, they assume error to aspatial (Hunter and Goodchild, 1997; Carlisle, 2005; Fisher and Tate, 2006; Wechsler, 2007). Invoking Tobler's First Law of Geography, whereby he noted that "nearby things are more similar than distant things" (Tobler, 1970), we know error varies spatially so DEM error is spatially autocorrelated. Indeed, Holmes et al. (2000) observe that "although global (average) error is small, local error values can be large, and also spatially correlated." The spatial variation of DEM error is most frequently estimated by calculating accuracy statistics of areas disaggregated by slope and/or landcover class, and more rarely spatial structure of error.

Wise (2000) categorized DEM errors as systematic, blunders or random. These types of errors derive from: (a) deficient spatial sampling and/or age of data; (b) processing errors such as interpolation or numerical errors; (c) measurement errors from poor positional inaccuracy, faulty equipment or observer bias (Wechsler, 2007). Systematic errors occur in the DEM creation procedure by processing techniques that can cause bias or artifacts. Blunders arise from human error (Wise, 2000) or equipment failure (Fisher and Tate, 2006). Random errors occur in any system of measurement due to the wealth of measurement and operational tasks performed to create a DEM (Wise, 2000; Fisher and Tate, 2006), and remain even after known blunders and systematic errors are removed (Wechsler, 2007). An example of random error is speckle noise (multiplicative noise in a granular pattern) (Rodriguez et al., 2006; Farr et al., 2007). Sources of systematic errors and blunders relevant to flood modeling derive from interpolation techniques (Desmet, 1997; Wise, 2007; Bater and Coops, 2009; Guo et al., 2010), erroneous sink filling (Burrough and McDonnell, 1998), hydrological correction (Callow et al., 2007; Woodrow et al., 2016), deficient spatial sampling causing urban features not to be resolved (Gamba et al., 2002; Farr et al., 2007), slope and aspect (foreslope vs backslope) (Toutin, 2002; Falorni et al., 2005; Shortridge and Messina, 2011; Szabó et al., 2015), striping caused by instrument setup (Walker et al., 2007; Tarakegn and Sayama, 2013) and vegetation (Carabajal and Harding, 2006; Hofton et al., 2006; Shortridge, 2006; Weydahl et al., 2007; LaLonde et al., 2010).

The aforementioned errors propagate into errors in surface derivatives including, but not limited to, slope (Holmes et al., 2000), aspect (Januchowski et al., 2010), curvature (Wise, 2011), drainage basin delineation (Oksanen and Sarjakoski, 2005) and upslope contributing area (Wu et al., 2008). As many models rely on these surface derivatives (e.g., change in slope is the dominant control on flow in flood models), error propagation from DEMs can substantially affect results of models that use these surface derivatives. Yet, taking flood models as an example, sensitivity analysis has largely focussed on hydraulic parameters and has under-represented DEM errors (Wechsler, 2007). Davis and Keller (1997) aptly sum up the problem of DEM error with their remark that 'landscapes are not uncertain, but knowledge about them is.'

# Flood Inundation Models and DEM Error

Topography is arguably the key factor for the estimation of flood extent (Horritt and Bates, 2002), but typically flood models use a limited number of DEMs and instead choose to explore the uncertainty associated with other hydraulic parameters (Wechsler, 2007). Studies that do use multiple DEMs either resample DEMs to a coarser resolution to explore the effect of resampling strategies and/or scale (Horritt and Bates, 2001; Neal et al., 2009; Fewtrell et al., 2011; Saksena and Merwade, 2015; Savage et al., 2016b; Komi et al., 2017), or compare flood extents using different DEM products (Li and Wong, 2010; Jarihani et al., 2015; Bhuyian and Kalyanapu, 2018). Generally speaking, the quality of flood predictions increases with higher resolution DEMs. Higher resolution DEMs are more important when modeling urban environments (Fewtrell et al., 2008) so buildings can be captured. Resolution can be less important for rural environments with Savage et al. (2016a) concluding that running simulations finer than 50 m had little performance gain without occurring additional unnecessary computational cost. Too much detail can induce spuriously precise results which


does not represent the uncertainties in making flood predictions (Dottori et al., 2013; Savage et al., 2016a). In data-sparse regions, a limited number of global DEM products dictates that only a single DEM is used, with this most commonly being SRTM (Yan et al., 2015). Whilst understandable, using a single DEM leads to a dangerous situation where spuriously precise estimates of flood extent are presented which do not assess the impact of uncertain topography. DEM simulation overcomes this obstacle

by making available a catalog of statistically plausible DEMs at the native resolution of the global DEM from which it is simulated from.

# CURRENT ADVANCES – CORRECTING DEM ERROR

Here we identify three categories of approaches to correct DEM error: (1) DEM editing; (2) New DEMs created with improved sensing technologies and (3) Stochastic simulation of DEMs. This article focuses on the third approach.

# DEM Editing

DEM error can be reduced by editing – either manually or systematically. Manual editing involves changing pixel values based on additional information or expert judgement. We speculate that this happens frequently but is seldom documented. Systematic editing involves applying algorithms, additional datasets and filters to reduce error. For example, a DEM can be hydrologically corrected through algorithms such as AGREE (Hellweger, 1997), ANUDEM (Hutchinson, 1989), outlet breaching (Martz and Garbrecht, 1999), Priority-Flood (Barnes et al., 2014) and stream burning (Saunders, 1999). Namely, the HydroSHEDs global hydrography dataset makes use of hydrological correction techniques to create invaluable maps such as flow direction, river networks and catchment masks (Lehner et al., 2008). DEMs have also been edited to correct errors from vegetation (Baugh et al., 2013; Pinel et al., 2015; Su et al., 2015; O'Loughlin et al., 2016; Ettritch et al., 2018; Zhao et al., 2018) and to compensate for the positive bias in coastal areas due to vegetation and buildings that lead to an underestimation of coastal flood exposure [e.g., CoastalDEM (Kulp and Strauss, 2018)].

The recent release of the MERIT (Multi-Error-Removed-Improved-Terrain) DEM is the most comprehensive error removal from SRTM to date. Errors are reduced by separating and removing absolute bias, stripe noise, speckle noise and vegetation bias, with the most significant improvements reported in flat regions (Yamazaki et al., 2017). Compared to SRTM, MERIT has fewer artifacts (Hirt, 2018) and a better performance in flood models compared to SRTM (Chen et al., 2018). Whilst a significant improvement on SRTM, MERIT is still fundamentally based on SRTM data and is thus limited by the errors in SRTM. The next new edited DEM will be NASADEM (Crippen et al., 2016) set to be released in 2018, but again this is a reprocessing of SRTM.

# New DEMs

The most widely used open-access global DEM (SRTM) is almost two decades old. Advances in satellite technology, image processing and data storage capabilities, make creating new, more accurate DEMs entirely possible. Two new global DEMs have been recently released, namely optically derived ALOS AW3D30 (Tadono et al., 2014) (∼30 m resolution) and the SAR derived TanDEM-X 90 (Rizzoli et al., 2017) (∼90 m resolution). Both of these products are technically digital surface models, so should only be used with caution in flood models. Looking to the future, new techniques are being explored to create new DEMs. One such example of applying a new technique is the creation of a 2 m pan-Arctic DEM (ArcticDEM)<sup>1</sup> using stereo auto-correlation techniques to overlap pairs of high-resolution optical imagery. Additionally, Ghuffar (2018) demonstrated that a 5 m DEM can be generated from Planet Labs cubesat derived PlanetScope imagery using Semi Global Matching. Alternatively, existing DEMs can be fused together to create new products (de Ferranti, 2014; Yue et al., 2017; Pham et al., 2018). For example ASTER and SRTM have been fused together to create the global EarthEnv DEM (Robinson et al., 2014). High Resolution (<10 m) open-access LiDAR data is becoming increasingly available through initiatives such as OpenTopography<sup>2</sup> , with New Zealand the latest country to release LiDAR data for free. Despite this encouraging step, we optimistically estimate openaccess LiDAR data covers just 0.005% of the earth<sup>0</sup> s land area based on data from OpenTopography and an extensive search of national mapping agencies. Global LiDAR coverage is some way off, with the limited amount of LiDAR data that is currently available almost exclusively found in developed countries.

# DEM Simulation

Stochastic simulation assumes that a DEM is only a single realization amongst a host of potential realizations. DEMs are simulated by altering pixel values in accordance with the spatial error structure. A single true DEM is not created. Instead realizations provide a bound within which the true value is likely to lie. Therefore, DEM error is not reduced as such, but the bounds of error are identified. This idea is relatively well known in the field of geostatistics (Goovaerts, 1997; Hunter and Goodchild, 1997; Deutsch and Journel, 1998; Kydriakidis et al., 1999; Holmes et al., 2000). Using an ensemble of simulated DEMs has been shown to greatly affect the characterization of surface derivatives (Fisher, 1991; Veregin, 1997; Holmes et al., 2000; Endreny and Wood, 2001; Raaflaub and Collins, 2006), landslide hazard (Davis and Keller, 1997; Murillo and Hunter, 1997; Darnell et al., 2008) and flood inundation estimation (Wilson and Atkinson, 2005; Hawker et al., 2018). Research in this area, especially in the flood community, has been largely stagnant for the past decade which seems a shame given the improvement in computational resources and the number of DEMs now available. Here we demonstrate the benefits of DEM simulation in flood inundation modeling based on the recently published work of Hawker et al. (2018).

#### DEM Simulation

In Hawker et al. (2018), DEM simulation is carried out by first quantifying the spatial error structure of a global DEM, and then using the fitted error covariance function to simulate plausible versions of the DEM. The fitted error covariance function was calculated for SRTM and MERIT DEMs by fitting a semivariogram to difference maps (i.e., SRTM/MERIT – reference LiDAR DEM) of 20 floodplain locations. Semi-variograms for

<sup>1</sup>https://www.pgc.umn.edu/data/arcticdem/

<sup>2</sup>http://opentopography.org/

all locations and by landcover class were produced, thus making DEM simulation possible for any floodplain SRTM/MERIT location, with the work available through the R Package DEMsimulation. For more details we refer the reader to Hawker et al. (2018). It must be noted that only a limited number of semivariograms are produced, but is nevertheless useful in creating an ensemble of DEMs.

We test the quality of DEM simulations produced by simulating MERIT DEM by landcover semi-variograms by plotting rank histograms for an ensemble of 2500 DEMs of the Ba catchment in Fiji. Rank histograms (or 'Talagrand' diagrams) (Anderson, 1996; Hamill and Colucci, 1997; Talagrand et al., 1997) are a common tool used to evaluate ensemble forecasts in meteorology, and work by ranking the verification (in our case LiDAR data) relative to the corresponding value in the ensemble in ascending order. An ideal ranked histogram is flat since the observation is indistinguishable from any ensemble member. Typically, a U-shaped rank histogram suggests undervariability in the ensemble, a dome shape over-variability, and excessive population of the extreme ranks bias. Yet, ranked histograms are notoriously difficult to evaluate and can lead to misinterpretations if done uncritically (Hamill, 2001). Nevertheless, we produce rank histograms by taking the mean of LiDAR values which fall within each ensemble pixel for the Ba catchment in Fiji (**Figure 1**). The rank histogram (**Figure 1B**) of all pixels suggests a positive bias in ensemble members as the ranks are clustered to the left. Despite the vegetation correction in MERIT, the rank histogram of mangrove covered pixels shows a large positive bias, whilst cropland has a more uniform shape.

To compensate for errors in observations (LiDAR), we added random observational noise as suggested by Hamill (2001), but this made little difference to the shape and thus is not presented here. Additionally, we compute 3 goodness-of-fit measurements: Pearson X<sup>2</sup> ; Jolliffe-Primo (JP) slope and JP convexity, with the null hypothesis that the rank histogram is flat (Jolliffe and Primo, 2008). These statistics confirm the stronger bias in mangroves (JP Slope) and suggest possible under-sampling with the relatively high JP convexity values. All these results are statistically significant with p-values of virtually 0. Moreover, less than 3% of pixels within the single MERIT DEM were within the error of the LiDAR (≈50 mm), whilst this was 97% for the ensembles. Therefore, the reliability of the DEM simulation is deemed satisfactory but can still suffer from systematic errors from the global DEM product being used. A higher-accuracy global DEM would therefore make this approach even more effective.

#### Simulated DEMs and Flood Inundation Predictions

To demonstrate the usefulness of using simulated DEMs in flood predictions we expand upon work published in Hawker et al. (2018). We simulate a total of 7500 DEMs to use in a LISFLOOD-FP Neal et al. (2012) flood model of Ba, Fiji (**Figure 2A**) for a 50 year return period flood event (Archer et al., 2018). The Ba catchment is predominantly agricultural floodplain, with mangroves present at the coast. DEMs are either simulated using MERIT or SRTM DEMs, and using either an average floodplain semi-variogram or semi-variograms disaggregated by landcover. Flood predictions are compared to four models that use a single DEM – LiDAR at 30 m and 90 m resolution and MERIT and SRTM at 90 m resolution. We assume the LiDAR 30 m model is the benchmark prediction in lieu of a lack of observation data. Flood depth errors are compared against the LiDAR 30 m model are plotted for the deterministic approach using the MERIT DEM and the stochastic approach using an ensemble of simulated DEMs (**Figure 2B**). Whilst the DEM ensemble approach can overpredict flood extent, flood depths are often more accurate as indicated by the more neutral colors of the DEM ensemble flood map given in **Figure 2B**). Further analysis of predicted flood depth (**Figures 2C–F**) indicate the benefit of using ensembles of simulated DEMs in predicting correct water depths. For example, in location 2, the MERIT DEM does not flood, whilst the flood depth in SRTM is large (>4.8 m), but for the ensembles of DEMs the distribution of predicted flood depths are more closely aligned with the flood depths predicted in the LIDAR models. The results also highlight the differences in predictions between DEMs, so we would encourage modelers to use multiple DEMs even if DEM simulation is not used. Yet by using an ensemble of simulated DEMs, we can learn about the distribution of potential flood extent and flood depth, and thus can avoid the spurious precision when using a single DEM.

# FUTURE DIRECTIONS

In this article, we have attempted to reinvigorate the idea of DEM simulation and highlight its value for flood studies. Despite repeated calls to produce a new high-accuracy open access global DEM (Schumann et al., 2014; Simpson et al., 2015), this unfortunately does not seem forthcoming. Everincreasing computing power has made even global flood simulations possible (Sampson et al., 2015), while flood modelers also often run multiple models to explore model parameter sensitivities. However, the impact of DEM error has been largely overlooked in lieu of a lack of suitable stochastic DEM data. DEM simulation overcomes this restriction, making it possible for flood modelers to use a catalog of DEMs. Working in tandem with systematic DEM editing (e.g., MERIT), DEM simulation can fill the gap until a much-needed new highaccuracy open access DEM is produced. Even when this longawaited DEM is eventually produced, DEM simulation will still be an invaluable approach for exploring the effect of DEM error in flood inundation estimates as long as good estimates of the spatial error structure can be made across a sufficient number of locations. We therefore encourage scientists to embrace geostatistics to simulate DEM ensembles and call for increased reporting of spatial dependence by DEM vendors and scientists alike.

# DATA AVAILABILITY STATEMENT

The MERIT dataset can be downloaded after sending a permission request to the developer (Dai Yamazaki, yamadai@rainbow.iis.u-tokyo.ac.jp) from http://hydro.iis.utokyo.ac.jp/~yamadai/MERIT\_DEM/, and is free for research

and education purposes. SRTM data can be freely downloaded from https://earthexplorer.usgs.gov/. The LiDAR dataset is available from The Secretariat of the Pacific Community<sup>0</sup> s Applied Geoscience and Technology Division (SPC SOPAC). LISFLOOD-FP is available for non-commercial purposes from http://www.bristol.ac.uk/geography/research/hydrology/models/ lisflood/downloads/. The code to simulate DEMs can be found at https://github.com/laurencehawker/demgenerator.

### AUTHOR CONTRIBUTIONS

feart-06-00233 December 14, 2018 Time: 14:38 # 7

LH wrote the article and carried out the data analysis. PB, JN, and JR reviewed and provided comments on various drafts as well as advised on the analysis work.

# REFERENCES


# FUNDING

LH was funded as part of the Water Informatics Science and Engineering Center for Doctoral Training (WISE CDT) under a grant from the Engineering and Physical Sciences Research Council (EPSRC), grant number EP/L016214/1. PB was supported by a Leverhulme Research Fellowship and Royal Society Wolfson Research Merit award.

### ACKNOWLEDGMENTS

We would like to thank Leanne Archer at the University of Bristol for her work on the flood model for Ba and Dai Yamazaki for his help with the MERIT dataset. The DEMsimulation R package is available at https://github.com/laurencehawker/DEMsimulation.



Sampson, C. C., Smith, A. M., Bates, P. D., Neal, J. C., and Trigg, M. A. (2016). Perspectives on open access high resolution digital elevation models to produce global flood hazard layers. Front. Earth Sci. 3:85. doi: 10.3389/feart.2015.00085


analysis of a flood inundation model. Water Resour. Res. 52, 9146–9163. doi: 10.1002/2015wr018198


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Hawker, Bates, Neal and Rougier. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Utilizing Flood Inundation Observations to Obtain Floodplain Topography in Data-Scarce Regions

Apoorva Shastry1,2 \* and Michael Durand1,2

<sup>1</sup> School of Earth Sciences, The Ohio State University, Columbus, OH, United States, <sup>2</sup> Byrd Polar and Climate Research Center, The Ohio State University, Columbus, OH, United States

Flood models predict inundation extents, and can be an important source of information for flood risk studies. Accurate flood models require high resolution and high accuracy digital elevation models (DEM); current global DEMs do not capture the topographic details in floodplains, and this often leads to inaccurate prediction of flood extents by flood models. Flood extents obtained from remotely sensed data provide indirect information about topography. Here, we attempt to use this information along with model predictions to produce better floodplain topography. The algorithm we describe is a two-step process: first, we reduce the noise along the observed flood boundaries for all particles. Then, the model predictions from these modified DEMs are assimilated with observations using a particle batch smoother. We implemented the algorithm for a synthetic test case. For the nominal case, we observed a significant improvement in accuracy in terms of RMSE (35% reduction), bias (20%), and standard deviation (40%). We conducted sensitivity analysis by using priors of varying bias (0.5, 1, and 2 m) and standard deviation (1, 2, and 4 m). The bias reduced to ∼0.5 m or below in all the cases: the reduction in bias varied from 11 to 76%. The standard deviation of errors in the final estimate was almost half of the prior: the reduction varied from 40 to 49%. The reduction in RMSE ranged between 35 and 67%. For the case with 2 m bias and 4 m standard deviation (SRTM-like error levels), bias went down to 0.48 m (76% reduction), and standard deviation reduced to 2.24 m (44% reduction). Flood inundation maps produced from the final estimate DEMs also improved on its prior. For the 2 m bias cases, true positive rate (TPR) for peak inundation went from ∼30% to more than 57% in all three cases. The algorithm produces promising results, and this type of analysis can be performed in data-poor floodplains where high resolution DEMs do not exist.

Keywords: digital elevation model, flood modeling, data assimilation, remote sensing, floodplains

# INTRODUCTION

Prediction of inundation extents from flood models is an indispensable source of information for assessing flood risk in the context of hazard studies, but inundation prediction accuracy is often limited by the quality of topographic data available globally. Topographic data, generally in the form of digital elevation models (DEMs), are the primary input data for flood inundation

#### Edited by:

Paul Bates, University of Bristol, United Kingdom

#### Reviewed by:

Andrew Mark Ireson, University of Saskatchewan, Canada Juan Pablo Rodríguez Sánchez, University of Los Andes, Colombia

> \*Correspondence: Apoorva Shastry shastry.7@osu.edu

#### Specialty section:

This article was submitted to Hydrosphere, a section of the journal Frontiers in Earth Science

Received: 06 October 2018 Accepted: 17 December 2018 Published: 11 January 2019

#### Citation:

Shastry A and Durand M (2019) Utilizing Flood Inundation Observations to Obtain Floodplain Topography in Data-Scarce Regions. Front. Earth Sci. 6:243. doi: 10.3389/feart.2018.00243

**43**

FIGURE 1 | (A) True DEM for the Buscot case and (B) maximum flooded extent.

modeling. Airborne light detection and ranging (lidar) DEMs offer the best horizontal resolution and vertical accuracy. However, high resolution lidar DEMs are not available globally and are expensive to obtain. Globally, the best available DEMs are obtained from satellite data: the shuttle radar topography mission (SRTM, spatial resolution 30 m) (Rodriguez et al., 2006) and multi-error-removed improved-terrain (MERIT, spatial resolution 90 m) (Yamazaki et al., 2017) DEM are a few examples of freely available DEMs. The global average for vertical relative error in SRTM is about 6 m (Rodriguez et al., 2006). TanDEM-X DEM produced by DLR (German Aerospace Center) has better spatial resolution (10–12 m) and vertical accuracy (∼2 m) (Krieger et al., 2006; Eineder et al., 2012). However, it is a commercial product and is not freely available. Because the vertical inaccuracies of the SRTM DEM, flood inundation modeling using SRTM produces spatial inconsistencies (Sanders, 2007; Yamazaki et al., 2012; Yan et al., 2015; Fernández et al., 2016; Sampson et al., 2016). Hence, these DEMs are not suitable for flood simulations as they do not represent the topography well (Sanders, 2007). There have been calls to produce global scale high resolution DEMs because of its impact on emergency services and scientific research (Schumann et al., 2014).

Spaceborne remote sensing observations of inundation extent contain indirect information about floodplain topography. Remotely sensed data is widely used to study floods (Schumann and Domeneghetti, 2016) for applications such as flood risk assessment, emergency flood response, and flood mapping, but inference of floodplain topography from inundation has rarely been attempted; Mason et al. (2016) is an example of one such study. Spatiotemporal inundation patterns in floodplains are responses to a combination of factors and input variables like flow rate in the river channels, soil type, vegetation, etc.: floodplain topography determines where floodwaters flow after rivers flood their banks. For this reason, inundation patterns contain indirect information about floodplain topography. As a specific example, the boundaries of flood inundation are essentially the contours of ground elevation, if the water level is assumed to be flat (this idea can easily be extended to account for river slope, as well). Using inundation "contours" in this way is the inverse of the wellknown body of the literature that uses inundation intersected with a high precision DEM to infer water levels (Matgen et al., 2007; Cohen et al., 2018). It is intuitive that these contours could be used to constrain relative variations in floodplain topography: for example, noisy elevations could be smoothed by averaging along inundation edges. Mason et al. (2016) used synthetic aperture radar (SAR) derived flood extents to improve the height accuracy in TanDEM-X DEM. They tested the method in a region where the errors had a mean and standard deviation of ∼ 0.5 and ∼2 m, respectively. The mean difference between TanDEM-X and lidar DEM reduced from 0.5 to 0.3 m, and their standard deviation reduced from 2 to 1.2 m. However, it is not certain the method would work for regions with higher height errors, like SRTM, which has mean and standard deviation of errors in the range of ∼ 1.2 and ∼ 4 m, respectively.

Because inundation images reflect the complex flow paths that water takes during flooding events that can only be captured by flood models, methods such as Mason et al. (2016) (referred to as "smoothing methods," hereafter) are unlikely to extract all possible topographic information from inundation imagery. Because two-dimensional flood models encapsulate floodplain processes, it is natural to attempt to use such models to help extract topographic information from inundation. Indeed, the

objective here is to infer floodplain topography using inundation maps, while flood models do the inverse: predict inundation using floodplain topography. Thus, this is a classic example of an "inverse problem" (Yeh, 1986). Data assimilation is a technique that can be used to solve inverse problems; assimilation combines measurements and models to produce an estimate of the system which is better than just the measurement or model alone.

To our knowledge, assimilation has not been used to estimate floodplain DEMs, though related work has been done. Durand et al. (2008) assimilated water surface elevation observations along with LISFLOOD-FP hydrodynamic model to estimate channel bathymetry. Observations from SAR images of river inundation were assimilated with 2-D shallow water equations to identify optimal Manning's coefficients (Lai and Monnier, 2009; Hostache et al., 2010; Monnier et al., 2016). One obvious impediment to using flood models to infer floodplain topography is the high dimensionality of the problem: in principle, the elevation value for each grid cell in the floodplain model must be estimated. An additional consideration is that flood models are computationally expensive, and many assimilation algorithms require either repeated model runs or ensembles of runs. A final concern is the high degree of non-linearity between the observed inundation and the floodplain topography, as some approaches [e.g., the ensemble-based assimilation of Durand et al. (2008)] are based on linear estimation theory. The non-linearity can be accounted for by using a so-called "Particle Batch Smoother" (PBS), as described in Margulis et al. (2015), which is a non-Gaussian estimator that directly approximates Bayes theorem.

In the present study, we present and test (using synthetic observations) a new algorithm designed to infer floodplain topography using globally available DEMs and inundation imagery. The algorithm consists of two steps: smoothing and data assimilation, that capitalize on the strengths of each method.

#### EXPERIMENT DESIGN

We tested our algorithm for a synthetic case of small domain. This allowed us to explore the sensitivity of the algorithm to errors of various magnitude. We used the Buscot model, a tested example model distributed with LISFLOOD-FP, as our synthetic test case. LISFLOOD-FP (Bates and De Roo, 2000; Bates et al., 2010) is a raster-based inundation model. We considered the DEM used in this model as the truth. The model had a defined river channel, and the flow in it was defined by the inflow boundary condition. The flow into the channel was defined as a triangular hydrograph with a flow of 20 m<sup>3</sup> /s at time 0, 200 m<sup>3</sup> /s at its peak and back to 20 m<sup>3</sup> /s at the end of the 5-day simulation. **Figure 1A** shows the DEM, and **Figure 1B** shows the maximum flooded extent for this test case.

To obtain best results from this method, we require multiple unique flood extent observations. We used 9 flood inundation maps obtained between day 1 and peak inundation on day 3 as observations. Our primary objective was to test the efficacy of the algorithm itself, and we did not focus on studying the impact of less or more flood inundation observations. In order to focus on the effect of prior DEM error on the analysis, we did not add white noise to the classified imagery; we leave for future work how observational uncertainty and temporal revisit would impact the algorithm accuracy. **Figure 2** shows flood inundation area for the simulation period, and the maps that were used as observations.

In the design of the synthetic experiment, we attempt to simulate a realistic situation where we attempt to correct a noisy "prior" estimate of the floodplain DEM. We accomplish this by taking the DEM distributed with the Buscot model to be the "truth." We then create a prior estimate of the DEM by adding errors to the truth. Here we chose to add spatially uncorrelated errors when creating the prior; the level of the errors varies among the various cases. We define the "nominal case" to be addition of 0.5 m bias and 1 m standard deviation of errors, which is referred to as the prior henceforth. We performed sensitivity analysis by considering 8 additional cases, by exploring bias ranging from 0.5, 1, and 2 m, and standard deviation of 1, 2, and 4 m. The case with bias of 0.5 and 2 m is similar to the TanDEM-X errors, which was used by Mason et al. (2016). The cases with bias of 1 and 2 m, and standard deviation of 4 m is similar to SRTM.

We evaluated the performance of the algorithm by calculating bias, standard deviation of DEM errors and root mean squared error (RMSE) for all the pixels modified by the algorithm. We also evaluated the DEM's ability to predict inundation by using true positive rate (TPR), a statistical measure of binary classification. TPR is the proportion of predicted inundation area from the model that is accurate (from observations).

# MATERIALS AND METHODS

Our approach merges smoothing and data assimilation to better extract floodplain topography information from inundation maps. We will use the PBS concept described by Margulis et al. (2015) for a different estimation problem. The PBS is related to the Ensemble Kalman batch smoother used by Durand et al. (2008) to infer channel bathymetry from inundated area, but relies upon a fully Bayesian approach that does not assume linearity between the observations and the quantities to be estimated. The PBS represents the relationships between observations and model parameters, as well as the associated uncertainty, using a set of model simulations, which are referred to as "particles." Each particle represents an independent random realization of the sample. The basic idea is to generate a moderately sized set of particles, where each particle is a random perturbation of the prior DEM. Here we adapt the usual PBS algorithm by first performing smoothing on each of the particles prior to the PBS estimation. Then the LISFLOOD-FP model is run on each of the smoothed particles. Finally, the PBS estimate of the DEM is computed by comparing the LISFLOOD-FP model run output of each particle to the map of inundation observation. The PBS estimate can be thought of as a weighted average of all the DEM particles, where the weights are based upon the accuracy of LISFLOOD-FP in simulating the inundation. **Figure 3** shows an outline of the method we use to estimate an updated DEM using the PBS. One note of clarification: we use weighted averages in two contexts. In the first case, they are used for the smoothing elevations along flood boundaries, and they are referred to as weights (ensemble and regression). In the data assimilation step, we use particle weights to obtain the final estimate, and these weights are referred to as particle weights.

# Generating Particles

It has been established from empirical evidence that DEM errors are not completely random, and often have spatial correlation (Hunter and Goodchild, 1997; Fisher, 1998; Kyriakidis et al., 1999; Carlisle, 2005; Erdogan, 2010 ˇ ). Hence, we added spatially correlated errors to the prior (It should be noted that, in this implementation of a synthetic case, we created the prior by adding uncorrelated random errors to the "true" DEM.), and generated an ensemble of 50 particles. We added errors with zero mean, 1 m standard deviation and 250 m correlation length to the prior. We also added a constant bias to each particle; the constant bias was randomly generated in the same range as the bias of the prior (−0.5 to 0.5 m for the nominal case).

The set of randomly perturbed DEMs thus obtained might be inefficient as the process may generate many unrealistic DEM particles. Data assimilation style approach will require a large set of particles to capture the complex spatial pattern of topography, making the process computationally inefficient. One way to deal with this problem is to make the ensemble of particles more realistic. We make this ensemble of particles more realistic by using the process described in Section "Smoothing Along Flood Boundary".

#### Smoothing Along Flood Boundary

We exploited the indirect information about elevations along the flood boundary to produce a set of realistic particles. For each particle in the ensemble, we extracted the elevations along the flood boundary. We extracted two adjacent pixels along the boundary: one flooded pixel on the edge (the "wet boundary"), and the adjacent pixel on the non-flooded side (the "dry boundary"). Then, we performed linear regression along both boundaries, considering the wet and dry boundaries separately. We consider the length along the boundary going from upstream to downstream as the independent variable and the extracted elevations as the dependent variable. We modified each particle by updating the elevations along the flood boundary as the weighted average of extracted elevations, and its regressions. By doing this, the noise in elevations along the boundary was reduced (**Figure 4**). The weights are expressed as the inverse of RMSE of the prior, and of the regression. We used nine flood inundation maps to perform this operation.

We used the estimate of bias and standard deviation of errors for the prior to calculate RMSE of the prior. We use this RMSE to compute the ensemble weight, wen, as the inverse of RMSE of the prior. For the regressed line along the boundary, the weight is defined as:

$$\omega\_{\text{fit},j} = \frac{1}{rms\left(z\_{en} - z\_{\text{fit},j}\right)}\tag{1}$$

where wfit,<sup>j</sup> is the regressed weight of jth particle, zen is the elevation along flood boundary in the prior, and zfit,<sup>j</sup> is the elevation along flood boundary in the jth particle. **Figure 4** shows RMS of difference between regression (zfit,<sup>j</sup> ) and prior (zen) elevations along the boundary, and its corresponding weight (wfit,<sup>j</sup> ) for four sample particles. The updated elevation along the flood boundary, z + j , is calculated as:

$$z\_{j}^{+} = \frac{\varkappa\_{en}z\_{en,j} + \varkappa\_{\hat{f}t,j}z\_{\hat{f}t,j}}{\varkappa\_{en} + \varkappa\_{\hat{f}t,j}} \tag{2}$$

**Figure 4** shows the elevations along the boundary for the particle, regressed and updated elevations for four sample particles in the nominal case. The flood boundary of each ensemble member was replaced with this value of z + j from that particle to obtain an updated set of particles.

When we modify elevations along the observed flood boundary, we introduce sub-optimality into the analysis, because the errors in the smoothed DEMs are now dependent on errors in the observations. Thus the errors in the LISFLOOD-FP predictions are also correlated with the errors in the observations, whereas the PBS assumes that they are not correlated. However, we assume that the degree of sub-optimality introduced by using the observations is relatively small compared to the large errors in the prior DEM.

#### Data Assimilation

We ran a forward simulation of LISFLOOD-FP using the updated particles (z + j ) to obtain inundation maps for each particle. All the parameters (except the DEM) in the model were the same as in the simulation using "true" DEM. The inundation maps from this set of simulations were used along with observed inundation maps using PBS (Durand et al., 2008; Margulis et al., 2015) to estimate the updated DEM. PBS is a non-Gaussian estimator which directly approximates Bayes theorem. Initially, all ensemble members are assigned equal particle weights. We evaluate the likelihood of each particle by simulating flood inundation maps, and calculating the agreement between the modeled and observed flood inundation. True positive rate (TPR) is used to define the agreement between the modeled and observed inundation. We used exponential probability distribution to update the weights. The value of the probability density function w<sup>j</sup> at any point tj=1−TPR<sup>j</sup> is defined as the particle weight of the jth particle. The probability density function for exponential distribution is defined as:

$$
\lambda w\_j = \lambda e^{-\lambda t\_j} \quad \text{for } t\_j > 0 \tag{3}
$$

where t<sup>j</sup> is the random variable, λ is a rate parameter and λ = 1 µ = 1/σ, µ being the expected value of the distribution and σ the standard deviation. We then update the particle weight according to its likelihood of producing accurate inundation maps. The rate parameter λ is used to represent uncertainty in observations, and we use a value of 0.1 in the current study. Here, we calibrated the value of λ for the data assimilation to produce largest error correction. This means that our confidence in observations is high. The value of λ can be modified to represent the confidence in the observations. Higher value of λ produces less contrast between weights at t equal to 0 and 1. Hence, the weights are close to each other, indicating a lower confidence in observations. The updated particle weight is proportional to the likelihood of the model predicting the observation. Particles that produce inundation maps with high agreement have large particle weights; particles that have poor agreement have small particle weights. The posterior or updated DEM is the weighted average of the prior DEMs, calculated as shown in (4).

$$DEN^{+} = \frac{\sum \text{w}\_{j} DEM\_{j}^{-}}{\sum \text{w}\_{j}} \tag{4}$$

where DEM<sup>+</sup> is the updated estimate, DEM<sup>−</sup> <sup>j</sup> and w<sup>j</sup> are jth particle and its particle weight of jth particle. **Figure 5** shows the agreement in flood extent between the model and observations for a sample of nine particles. TPR, and the corresponding particle weight (wj) assigned to that particle for those nine particles are shown in their respective panel in **Figure 5**. It can be


seen that the particles that produce inundation maps with high agreement are given high particle weights, and vice versa. The final estimate DEM is obtained by using these particle weights to obtain a weighted mean of the updated ensemble of particles.

# Sensitivity Analysis

In most continents, the mean SRTM error is less than 2 m, and standard deviation is ∼4 m (Rodriguez et al., 2006). To ensure that our methodology works in these ranges, we implemented

the algorithm for several cases with varying levels of error mean and standard deviation in the prior. We created eight additional priors to the nominal case. The nine cases had errors of mean 0.5, 1, and 2 m, and standard deviation of 1, 2, and 4 m added to the "true" DEM. We applied the algorithm described in Section "Generating Particles, Smoothing Along Flood Boundary, and Data Assimilation" to the eight additional priors to obtain the final estimate DEM for those respective cases.

# RESULTS AND DISCUSSION

Root mean squared error for the prior for the nominal case (errors of 0.5 m mean and 1 m standard deviation) was 1.10 m. The set of particles generated from this prior had an RMSE of was 1.10 m. The RMSE of height errors along the flood boundaries reduced to 0.78 m in the updated set of particles where we smoothed elevations along the flood boundary. **Table 1** shows height errors in the initial ensemble and updated ensemble of particles. After putting this modified set of particles through the PBS, the RMSE further reduced to 0.71 m. The bias reduced from 0.5 to 0.4 m and standard deviation of errors went down from 0.98 to 0.59 m.

Sensitivity analysis also showed similar trends, improving upon priors in terms of bias, standard deviation of errors and RMSE. We found that the process of smoothing did not have an effect on the bias before and after smoothing. The difference between prior bias and smoothed bias is less than 4 cm for all the cases (**Table 1**). However, the standard deviation of errors went down after smoothing. The standard deviation of error was reduced by 40 to 49% in all the cases (**Table 1**). Generally, as the standard deviation in the prior increased, the amount of reduction increased. The amount of reduction also increased as the bias in the prior increased. **Table 1** shows the values of standard deviation for all the tested cases for the prior and smoothed ensembles.

When this smoothed set of particles was put through a PBS, the bias reduced in the final estimate from the prior and the smoothed ensemble. **Table 1** shows the details about error for the final estimate of elevations along the flood boundary. In general, the relative reduction in bias increased as the prior's bias increased. For the cases with 0.5 m bias, reduction in bias for the final estimate ranged between 11 and 24%. In all other cases (1 and 2 m bias), the bias reduced to 0.5 m or less. The reduction in bias ranged between 52 and 76%. However, there was no effect of the PBS on the standard deviation of error. The reduction is standard deviation of errors was less than 4 cm for all the cases (**Table 1**).

**Figure 6** shows the histograms for errors in all tested cases. It is clear from **Figure 6** that there is a reduction in bias and standard deviation of errors from the prior to the final estimate. A greater number of cells have lower errors, and the number of cells with

TABLE 2 | True positive rate (TPR) at peak inundation.


high errors are reduced. The algorithm improved upon the prior in all the cases, and the reduction in RMSE varied from 35 to 67%.

When we used the final estimate DEMs to predict flood inundation, there was a consistent increase in TPR when compared to the prior. **Figure 7** shows the predicted inundation area for the simulation period. In all the cases, the correctly predicted flood inundation area (i.e., the total area where inundation is in observation and prediction) was greater for the final estimate than the prior. **Figure 8** shows the flood inundation maps obtained from the prior and final estimated DEMs for two cases. The top row shows the peak inundation maps for the case with error of 2 m mean and 1 m standard deviation on 2.3 days. The TPR goes up from 33 to 82%. The second row corresponds to the case with error of mean 2 m and standard deviation 4 m, the case we believe is similar to SRTM DEM. In this case, the TPR increases from 29 to 57%. **Table 2** shows the change in TPR for peak inundation in all the cases. There was greater improvement of TPR with increase in bias of the priors. TPR for the peak inundation obtained from the final estimate DEM for all the cases with 1 m error standard deviation was greater than 82%. For the cases with 2 m error standard deviation, TPR for the peak inundation obtained from the final estimate DEM was between 63 and 75%; for cases with 4 m standard deviation of errors, it ranged between 52 and 62% (**Table 2**). It can be seen than there is a significant increase in TPR from prior to final estimate in all cases.

#### CONCLUSION

We successfully implemented a new algorithm to improve topography information in a floodplain by exploiting indirect information of ground elevations from observed flood extents. In synthetic tests, the algorithm reduced the bias, standard deviation of errors and RMSE. Our primary motivation to produce better topography was to obtain DEMs that are more suitable for flood inundation simulations. The improved DEM we obtained from this algorithm also predicted flood inundation much better than

#### REFERENCES


the prior. We implemented the algorithm for nine different cases with varying mean and standard deviation of errors, and obtained similar trends in the reduction of bias and standard deviation of errors. In fact, the magnitude and percentage reduction in bias increases in cases with higher errors. The results from the synthetic tests show potential, and we believe that the method could be used to improve DEM accuracy. For example, SRTM DEM could be used as prior, along with flood inundation observations obtained from Landsat or radar to obtain a DEM with better elevation accuracy.

Digital elevation models are the primary source of topographic information, and accurate DEMs are hard to obtain in the developing world. Globally available open-source products are easy to obtain, but are not accurate. Hence, they not suitable for flood inundation modeling (Fernández et al., 2016). It is not always physically or financially feasible to obtain lidar DEMs in data-poor regions. An algorithm to improve already available DEMs is one way we can study these regions, and make better predictions. This study has produced promising results, and we believe this algorithm can be applied to real world cases to improve floodplain topography. This will provide us with higher accuracy DEMs in data-poor floodplains which are suitable for flood inundation simulations.

### AUTHOR CONTRIBUTIONS

All authors designed the model and the computational framework and wrote the manuscript. AS performed the analysis and interpretation of data.

#### FUNDING

This work was funded by NASA Headquarters under the NASA Earth and Space Science Fellowship Program (Grant 80NSSC17K0322).

International Geoscience and Remote Sensing Symposium (IGARSS), (Munich), doi: 10.1109/IGARSS.2012.6351130


in Proceedings European Conference on Synthetic Aperture Radar (EUSAR), (Dresden).


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling Editor declared a past co-authorship with one of the authors MD.

Copyright © 2019 Shastry and Durand. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Anticipated Improvements to River Surface Elevation Profiles From the Surface Water and Ocean Topography Mission

Theodore Langhorst<sup>1</sup> \*, Tamlin M. Pavelsky<sup>1</sup> , Renato Prata de Moraes Frasson<sup>2</sup> , Rui Wei<sup>2</sup> , Alessio Domeneghetti<sup>3</sup> , Elizabeth H. Altenau<sup>1</sup> , Michael T. Durand<sup>4</sup> , J. Toby Minear<sup>5</sup> , Karl W. Wegmann6,7 and Matthew R. Fuller<sup>8</sup>†

<sup>1</sup> Department of Geological Sciences, The University of North Carolina at Chapel Hill, Chapel Hill, NC, United States, <sup>2</sup> Byrd Polar and Climate Research Center, The Ohio State University, Columbus, OH, United States, <sup>3</sup> Department of Civil, Chemical, Environmental, and Materials Engineering – DICAM, University of Bologna, Bologna, Italy, <sup>4</sup> School of Earth Sciences, The Ohio State University, Columbus, OH, United States, <sup>5</sup> Cooperative Institute for Research in Environmental Sciences, University of Colorado, Boulder, Boulder, CO, United States, <sup>6</sup> Department of Marine, Earth, and Atmospheric Sciences, North Carolina State University, Raleigh, NC, United States, <sup>7</sup> Center for Geospatial Analytics, North Carolina State University, Raleigh, NC, United States, <sup>8</sup> Nicholas School of the Environment, Duke University, Durham, NC, United States

#### Edited by:

Paul Bates, University of Bristol, United Kingdom

#### Reviewed by:

Kostas Andreadis, University of Massachusetts Amherst, United States Jeffrey Neal, University of Bristol, United Kingdom

> \*Correspondence: Theodore Langhorst tlang@live.unc.edu

#### †Present address:

Matthew R. Fuller, ORISE Postdoc Fellow at U.S. Environmental Protection Agency, Atlantic Ecology Division, Narragansett, RI, United States

#### Specialty section:

This article was submitted to Hydrosphere, a section of the journal Frontiers in Earth Science

Received: 01 October 2018 Accepted: 24 April 2019 Published: 08 May 2019

#### Citation:

Langhorst T, Pavelsky TM, Frasson RPdM, Wei R, Domeneghetti A, Altenau EH, Durand MT, Minear JT, Wegmann KW and Fuller MR (2019) Anticipated Improvements to River Surface Elevation Profiles From the Surface Water and Ocean Topography Mission. Front. Earth Sci. 7:102. doi: 10.3389/feart.2019.00102 Existing publicly available digital elevation models (DEMs) provide global-scale data but are often not precise enough for studying processes that depend on small-scale topographic features in rivers. For example, slope breaks and knickpoints in rivers can be important in understanding tectonic processes, and riffle-pool structures are important drivers of riverine ecology. More precise data (e.g., lidar) are available in some areas, but their spatial extent limits large-scale research. The upcoming Surface Water and Ocean Topography (SWOT) satellite mission is planned to launch in 2021 and will provide measurements of elevation and inundation extent of surface waters between 78◦ north and south latitude on average twice every 21 days. We present a novel noise reduction method for multitemporal river water surface elevation (WSE) profiles from SWOT that combines a truncated singular value decomposition and a slope-constrained least-squares estimator. We use simulated SWOT data of 85–145 km sections of the Po, Sacramento, and Tanana Rivers to show that 3–12 months of simulated SWOT data can produce elevation profiles with mean absolute errors (MAEs) of 5.38–12.55 cm at 100–200 m along-stream resolution. MAEs can be reduced further to 4–11 cm by averaging all observations. The average profiles have errors much lower than existing DEMs, allowing new advances in riverine research globally. We consider two case studies in geomorphology and ecology that highlight the scientific value of the more accurate in-river DEMs expected from SWOT. Simulated SWOT elevation profiles for the Po reveal convexities in the river longitudinal profile that are spatially coincident with the upward projection of blind thrust faults that are buried beneath the Po Plain at the northern termination of the Apennine Mountains. Meanwhile, simulated SWOT data for the Sacramento River reveals locally steep sections of the river profile that represent important habitat for benthic invertebrates at a spatial scale previously unrecognizable in large-scale DEMs presently available for this river.

Keywords: SWOT simulator, DEM, river water surface elevation, elevation profile smoothing, satellite altimetry

# INTRODUCTION

feart-07-00102 May 7, 2019 Time: 16:47 # 2

Accurate measurements of river water surface elevation (WSE) and slope at fine spatial scales are useful for monitoring river discharge (LeFavour and Alsdorf, 2005), calculating stream power in erosional models (Whipple and Tucker, 1999), interpreting underlying geology (Schumm, 1986), identifying knickpoints (Hayakawa and Oguchi, 2006), and characterizing habitat fragmentation for freshwater fish (Dias et al., 2013). In addition, water surface slope is a good predictor of physical habitat classifications (Jowett, 1993), which have different abundances and compositions of benthic taxa (Brown and Brussock, 1991). In situ stream gauges provide high-accuracy measurements of water elevations at discrete points but lack spatial continuity. Even in the most heavily montiored parts of the world, the gauge network is sparse and becoming sparser (Hannah et al., 2011). This problem is worse in less developed areas like the Arctic (Shiklomanov et al., 2002). Meanwhile, existing global DEMs are insufficient over open water due to missing data, large vertical errors, coarse spatial resolution, or limited temporal resolution (Alsdorf et al., 2007).

Remotely sensed measurements of river surface elevation have been made using a variety of methods that each have their limitations (Schumann and Domeneghetti, 2016). GEOSAT, a radar altimeter launched by the U.S. Navy in 1985, was shown to have root mean square error (RMSE) of 0.7 m compared to gauge stations on the Amazon River (Koblinsky et al., 1993). The TOPEX/Poseidon altimeter was designed to measure ocean topography but is also able to observe rivers with widths larger than 1 km with RMSE of 1.1 m (Birkett et al., 2002). The Jason-2 radar altimeter, part of the Jason series succeeding TOPEX/Poseidon, has been shown to have standard errors of 0.28 and 0.19 m over the Ganga and Brahmaputra rivers, respectively, but, like its predecessors, works best in rivers more than a kilometer wide (Papa et al., 2012). Satellite altimeters can provide high-accuracy WSE measurements, but with relatively poor spatial resolution compared to imaging radars. On the other hand, airborne laser altimeters have been used to map WSEs at high spatial resolution. For example, Biron et al. (2013) found elevation errors with standard deviation of 25 cm at 20 m spatial resolution compared to differential GPS measurements. The Shuttle Radar Topography Mission (SRTM) C-band DEM has an estimated standard deviation of error of 5.5 m for open water and has extensive missing data because radar returns depend on the surface roughness (LeFavour and Alsdorf, 2005). TanDEM-X is a relatively new 12 m resolution global DEM from the German space agency (DLR), but it performs poorly over water for similar reasons to SRTM. As a result, the elevations for many northern latitude water bodies in TanDEM-X are collected during winter months when they are ice covered (Wendleder et al., 2013). The Japanese space agency (JAXA) released the ASTER Global DEM Version 2 in 2011, but water elevations are often missing due to failure of the stereo matching technique method over water (Tachikawa et al., 2011). ArcticDEM is a 2 m resolution photogrammetry-derived digital surface model that covers latitudes north of 60◦ plus the Kamchatka peninsula and all of Alaska<sup>1</sup> . Despite the challenges of photogrammetry over water, extracting elevations from the shoreline on the Tanana River produced an elevation profile with height error standard deviation of 0.30 m at 100 m resolution based on ArcticDEM (Dai et al., 2018). More precise data are available from national elevation products like the National Elevation Dataset (NED) in the United States (Gesch et al., 2002), and TINITALY in Italy (Tarquini et al., 2011), but their limited spatial extent prevents global-scale studies. The Multi-Error-Removed Improved-Terrain (MERIT) DEM is created from SRTM, AW3D30, and the Viewfinder Panoramas DEM data, and is processed to decrease random and systematic sources of error (Yamazaki et al., 2017). The accuracy of MERIT over rivers has not yet been evaluated, however, problems of the source DEMs over open water are likely a persistent source of error in this data product.

The Surface Water and Ocean Topography (SWOT) satellite mission, planned to launch in 2021, will provide measurements of elevation and inundation extent of surface waters between 78◦ north and south on average twice every 21 days (Biancamaria et al., 2016). Over its planned 3 years lifetime, SWOT will provide repeat measurements of each river wider than 50–100 m. Any single observation of a river reach will have relatively poor accuracy (∼0.5 m) at the 100–200 m scale required to identify many along-stream topographic features. However, by leveraging multitemporal SWOT data, it is possible to reduce vertical errors and produce global river elevation datasets of unprecedented accuracy in SWOT-observable rivers. Here we use a truncated singular value decomposition to reduce the measurement error in a set of simulated SWOT observations. Subsequently, we use a constrained least-squares estimator to ensure that elevations decrease in the direction of flow, and to reduce individual WSE profiles to an average longitudinal profile. We apply this new method to simulated SWOT WSE profiles of the Po, Sacramento, and Tanana Rivers to evaluate the error reduction. We also compare error statistics for average simulated SWOT elevation profiles and profiles extracted from existing DEMs. Last, we consider case studies in geomorphology and ecology to highlight the scientific value of the more accurate in-river DEMs expected from SWOT.

### MATERIALS AND METHODS

#### Simulator Data

Ahead of the launch of the SWOT satellite, the NASA Jet Propulsion Laboratory (JPL) created a software simulator that approximates the sampling and error characteristics of SWOT (Frasson et al., 2017; Domeneghetti et al., 2018). The SWOT simulator requires a time series of WSEs, usually derived from a hydrodynamic model, for inundated areas. Additionally, a static DEM of the surrounding topography is used to simulate layover errors. The simulator samples the modeled surface elevations temporally and spatially according to the planned SWOT orbit and adds errors from terrain layover, instrument thermal noise,

<sup>1</sup>http://arcticdem.org

TABLE 1 | Physical and SWOT simulation characteristics of the sample rivers.


and satellite positional uncertainties. Errors associated with vegetation and the impact of specular reflections from the water surface are not included in this study. The simulator outputs a pixel cloud of elevation, inundation type, and other properties that are then summarized to regularly spaced nodes along the river centerline (Frasson et al., 2017). We use hydrodynamic models and simulated SWOT overpasses of the Sacramento, Po, and Tanana Rivers described in previous publications (**Figure 1**; Altenau et al., 2017; Frasson et al., 2017; Domeneghetti et al., 2018). The hydrodynamic model of the Sacramento is a onedimensional HEC-RAS model with an average cross sectional spacing 1.9 river widths (**Table 1**; Frasson et al., 2017). The Po River is modeled using a quasi-two-dimensional HEC-RAS model based on bathymetric cross sections an average of 2.5 river widths apart, and combined lidar and SRTM DEMs outside the main channel (Castellarin et al., 2011). The Tanana River simulation uses a two-dimensional LISFLOOD-FP model over a 25 m interpolated bathymetric grid (Altenau et al., 2017).

The three river models vary in spatial and temporal extent, and the frequency of simulated SWOT observations varies according to overlap of the planned orbit and the river location. The Po River is the largest set of data in this study, with a full year of simulation, and the Tanana River is the smallest data set, spatially and temporally (**Figure 2** and **Table 1**). Additionally, the smaller node spacing on the Tanana increases noise at the node level, as fewer simulated point measurements are averaged for each node. Similarly, the noise is relatively high for the Sacramento River simulation as it is the narrowest river and fewer simulated observations are available for each node. We removed one high discharge observation from the Sacramento simulation, so that the remaining profiles better represent the average elevation profile. At high discharge, small-scale details of the elevation profile can be lost, or muted, as larger hydraulic controls extend their influence further upstream (Dingman, 2009). Input rasters to the SWOT simulator are sampled in space and time and summarized to nodes matching the simulator output, which allows direct comparisons of the input and output nodes for error analysis of SWOT.

#### DEM Data

We acquired existing DEMs in order to compare the anticipated errors in SWOT river elevation profiles to current elevation data, including SRTM, MERIT, ASTER, NED, ArcticDEM, TanDEM-X, and TINITALY DEMs, where available. For the

Sacramento River, we acquired lidar data from the California Valley Floodplain Evaluation and Delineation Program and created a 10 m resolution raster from the last returns (California Department of Water Resources, 2013). We have not included any satellite altimeters in our analysis as the spatial scales of these data are too coarse compared to the other DEMs. We sampled elevations from the DEMs using hand-drawn river centerlines densified to the raster resolution. Manually drawing centerlines for each river ensures elevations are sampled from within the channel, which can change location between data sources. Height errors and missing data in the DEMs caused automated methods of centerline delineation to fail in many cases, due to higher elevations within the channel than outside. Centerlines were manually identified based on changes in data quality over open water and from referencing satellite imagery. The sampled DEM elevations are projected onto average centerlines from the SWOT simulations and transformed to flow distance and cross-channel coordinates (Legleiter and Kyriakidis, 2008). This process allows direct comparison of elevations originally sampled from different centerlines using the common flow distance coordinate. Last, we upscale the DEM profiles to the spacing of the SWOT simulated nodes using a windowed mean and interpolate the DEM data at the flow distance of each simulated node (left side of **Figure 3**). The end result of this processing is a separate elevation profile from each DEM projected on to a common coordinate system. To analyze the errors, we compare the DEM profiles to the average hydrodynamic model output for all three rivers, as well as boatmounted GPS profiles of the Sacramento and Tanana Rivers. The GPS data collection is described by Minear and Wright (2016) for the Sacramento River, and by Altenau et al. (2016) for the Tanana River.

#### Profile Smoothing

There are many published methods to make noisy DEM surfaces more suitable for hydrologic research, including filling (Jenson and Domingue, 1988), carving (Soille et al., 2003), spline regression (Harbor et al., 2005), and slope-constrained quantile regression (Schwanghart and Scherler, 2017). One reach definition method being considered for SWOT operations uses a Gaussian smoothing filter to reduce noise on SWOT profiles. The Gaussian weighting function used in the filter has two formulations: one uses a static 2 km standard deviation, and the other uses 1/5 of the reach length and a minimum of 1 km standard deviation. Both of these smoothing operations result in very smooth profiles that accurately capture slopes at 10 km scales (Frasson et al., 2017), but the wide window of the filter means small-scale topographic details such as riffle-pool structures are often lost.

We present new methods to reduce noise from multitemporal river elevation profiles that rely on the commonalities between repeated measurements instead of spatial smoothing. To reduce the noise in the simulated elevation profiles, we first decompose the elevation data using the singular value decomposition and analyze the resulting eigenvectors (full description in section Low Rank Approximation). The data is then recreated using fewer eigenvectors, reducing variability between observed profiles, and highlighting the real variability of the profiles. We hypothesize that repeated measurements of river elevation can be closely approximated at lower rank. In other words, the elevation profiles are linearly dependent on one another. The matrix becomes full rank when the simulated SWOT noise is added. To reduce the noise, we eliminate many of the eigenvectors to create a lowrank approximation (LRA) of the simulated data. This LRA is then further constrained using a least-squares estimator such that node elevations decrease in the downstream direction (full description in section Slope Constraint).

#### Low Rank Approximation

The nodes output by the SWOT simulator are arranged in a matrix such that rows represent nodes along the river centerline, columns represent overpasses, and the values are the simulated elevations. Each simulated SWOT orbit track observes a different set of nodes, which results in an inconsistent number of observations for each node. For example, the middle reaches of the Sacramento river are only observed by half of the simulated overpasses (**Figure 1**). As a result, the matrix has missing values, and the decomposition does not have a unique solution in this case. To overcome this problem, the river is divided into sections such that all nodes in a given section have the same number of observations (similar to the sections defined by shades of blue in **Figure 1**). Before the decomposition, we remove the mean simulated elevation of each node, which allows us to analyze the relationship between observations, instead of the overall slope of the river. We decompose the data matrix A into two sets of eigenvectors and corresponding singular values that represent the weight of each vector pair as follows:

$$A = U \ast \mathcal{S} \ast V^{\mathrm{T}} \tag{1}$$

Where A is the data matrix, U is a set of eigenvectors that describe relationships between river nodes, V is a set of eigenvectors

that describe relationships between overpasses, and S is a diagonal matrix of singular values that relate U and V. The first eigenvector of U represents the most common representation of all the profiles and subsequent eigenvectors describe variations on that profile. The data is recombined from a subset of these eigenvectors:

$$
\tilde{A} = U\_{\mathbf{k}} \ast \mathbf{S}\_{\mathbf{k}} \ast V\_{\mathbf{k}}^{\mathrm{T}} \tag{2}
$$

Where A˜ is the best rank k approximation of A in a leastsquares sense, U<sup>k</sup> and V T k are the first k eigenvectors, and S<sup>k</sup> is a diagonal matrix containing the first k singular values (Eckart and Young, 1936).

The LRA is dependent on the number and combination of eigenvectors picked, and the addition of SWOT-like noise to our data makes this decision more difficult. The magnitude of the SWOT noise at this spatial scale means the singular values will not taper off to zero but instead gradually decrease. To identify a cutoff threshold for the singular values, we use parallel analysis (Horn, 1965) followed by a test for significance between orbits in the V matrix. We calculate an average singular value spectrum from 1000 realizations of random, normally distributed data using the sample standard deviation of the simulated SWOT errors. Singular values from the simulated SWOT data that are greater than the corresponding average singular value from the random data are retained.

As a last step in the factor analysis, we test each eigenvector in the V matrix to check if they distinguish between orbit tracks. Different orbits will observe the same river nodes at different ranges of the radar swath, resulting in different error characteristics for each orbit (Fernandez, 2017). The V matrix contains eigenvectors that scale the effect of the U eigenvectors for each overpass. We test each eigenvector in the V matrix to check for statistically significant distributions for each orbit

track using a Wilcoxon rank sum test at 95% confidence. We interpret any eigenvector in V that shows different distributions when grouped by orbit as errors related to the viewing angle, and not related to real changes in the water surface. The remaining eigenvectors are interpreted to represent the real variability in the elevation profiles as a result of discharge variability.

#### Slope Constraint

We also apply a second method of noise reduction: a slope-constrained least-squares estimator, requiring each node elevation to be less than or equal to the upstream node elevation. The constrained profiles are calculated by solving the following equations:

$$\text{minimize } ||\mathbf{G} \* \hat{z} - z||\_2$$
 
$$\text{such that: } \mathbf{C} \* \hat{z} \le b \tag{3}$$

Where z is the observed heights, zˆ is the constrained heights, G is a matrix that relates zˆ and z, b is a vector of zeros to represent the maximum allowable difference between downstream and upstream nodes, and C is a matrix that calculates the difference between nodes when multiplied by zxˆ . Solving for zxˆ gives the least-squares set of elevations that decrease in the downstream direction.

The formulation of the constrained estimator can produce unwanted results for long sections of observations that display negative slopes. As defined in Equation (3), the solution for these sections has zero slope, and is often followed by a steep slope to compensate (**Figure 4**). For most rivers, this is not a realistic representation of their slopes and would be problematic for some analyses. We present a modification of Equation (3) to counteract this effect:

$$\begin{array}{l}\text{minimize} \left\| \left\| \begin{bmatrix} G\\ \lambda \ast C \end{bmatrix} \ast \hat{z} - \begin{bmatrix} \bar{z} \\ \bar{d} \end{bmatrix} \right\| \right\|\_{2} \\\text{Such that: } \mathbf{C} \ast \hat{z} \leq b \end{array} \tag{4}$$

Where ¯d is the average change in elevation between observed nodes, and λ is a regularization parameter that controls the relative weight of the height and slope terms in the minimization problem. A high λ value penalizes the extreme slopes that can result from Equation (3), whereas a λ value of 0 makes Equations (3) and (4) equivalent. By adding regularization to the estimator, we can come closer to recreating the distribution of slopes we see in the hydrodynamic model (**Figure 4D**), however, the equation now requires careful parameterization. The example in **Figure 4D** was parameterized using our knowledge of the true profiles, and as such is optimistic. Using a λ value that is too low has little effect and using one that is too high results in an overly smooth profile with few details. This is the same problem encountered with other smoothing and local regression techniques: the user is forced to pick the smoothness of the resulting profile. The advantage of Equation (3) over these techniques is there is no decision to be made about smoothness, or assumption of the variability of slopes. The height errors are reduced only by the physical constraint that water flows downhill. As such, we use the output of Equation (3) for the remainder of our analysis.

We use the constrained least-squares estimator from Equation (3) in two ways: to constrain the slope of individual profiles and to estimate an average profile from many observations (**Figure 3**). For our error analysis of the multitemporal simulated SWOT profiles, the constrained least-squares estimator is solved for each overpass separately. These profiles are referred to as the Constrained LRA (CLRA) profiles. We also use the estimator to create an average profile by including all simulated observations in the vector z, and create G to relate all observations at one node to one predicted height in x. This creates one constrained profile from all observations, which we term the average SWOT profile. Combining elevation profiles from different days, which represent a range of discharges, is not a perfect solution to reducing noise, as elevations at different nodes will respond differently to variations in discharge. However, considering the large vertical errors in the raw data, the error introduced at this step are likely comparatively small.

### RESULTS

The application of the CLRA method to our simulated SWOT data sets decreased the mean absolute error (MAE) of every simulated profile when compared to the raw SWOT node elevations (**Figure 5**). Analyzing the errors at an intermediate step shows that both the low rank approximation and the slope constraint result in decreased height errors. Node-scale MAEs were reduced by more than half for all three rivers (**Table 2**). The Sacramento River showed the largest decrease in MAE, at more than 70% reduction. RMSEs of the node heights show similar improvement. The LRA component of the method reduces error variability between days, which is apparent in **Figure 5** as less spread among the days for each river on the vertical axis compared to the horizontal axis.

Decomposing the hydrodynamic model output shows us that the first eigenvector from noise-free data always contains more than 95% of the variance in the data for all three rivers, and that the first four eigenvectors always contain more than 99% of the variance. Our factor analysis method for simulated SWOT data resulted in just one eigenvector being used to recreate the data for all locations. Parallel analysis occasionally indicated two eigenvectors should be used to approximate the data, but the test for significant difference between orbits eliminated the second eigenvector in every case.

We applied the constrained least-squares estimator to all the DEM profiles in **Figure 6** and for evaluation in **Table 3**, which reduced errors in all cases. We calculate MAE and RMSE statistics for the DEM-derived and average SWOT profiles for all three rivers compared to the average hydrodynamic model elevation profiles, and for Sacramento and Tanana Rivers compared to the GPS profiles. Comparisons with the hydrodynamic model favor the SWOT simulation, as any inaccuracies in the hydrodynamic model will not affect the SWOT error estimates but will increase errors for the DEMs. Comparisons to the GPS profiles increase errors of the SWOT simulation profiles due to inaccuracies in

the hydrodynamic models. Even though there is no fully fair direct comparison, we can see the difference between SWOT and existing DEMs is quite large (**Figure 6**). The ASTER-derived profiles are relatively coarse, despite the improvements from the constrained least squares estimator. Estimating river surface slopes from these ASTER profiles would require very large reaches. The difference between SRTM and MERIT is small for the Sacramento and Po Rivers, which is to be expected as MERIT is derived from SRTM data in these areas (Yamazaki et al., 2017). The error removal methods used to create MERIT offer some improvement over SRTM for the Sacramento River, but errors are still quite large compared to the SWOT and lidar profile.

# APPLICATIONS

# Physical Habitats

Riffle-pool sequences in rivers represent distinct physical habitats that also vary ecologically with different invertebrate community composition, density, and biomass (Brown and Brussock, 1991). Yang (1971) defined riffles and pools according to their energy gradient, which can be approximated by the water surface slope, but others have based classifications on channel slope, as it is independent of discharge (Richards, 1976). In practice, classifications of riffle and pool habitats are often subjective, but quantitative analysis of these classifications shows a threshold on

FIGURE 5 | MAE before and after the CLRA method for each observed profile. Points below the 1:1 line indicate the CLRA method decreased the height errors.

TABLE 2 | Error statistics before and after the CLRA method.


surface slope can identify riffles and pools with accuracies of 70 and 79%, respectively (Jowett, 1993). We evaluate the ability of SWOT to find locally steep sections of rivers in comparison to existing elevation data.

We use relative steepness, originally used to identify knickpoints, to identify sections of the simulated Sacramento profiles that are potential riffles or runs. We selected the Sacramento for this application because the upstream reaches show cyclical low-high slope sections (e.g., **Figure 4**). We use the relative steepness metric in place of a threshold on slope, which would vary between rivers. Relative steepness is defined as the trend in surface slope with increasing reach lengths (Hayakawa and Oguchi, 2006). We use reach lengths from 400 to 20,000 m. High values of relative steepness indicate places where the local slope is very different from the slope over greater lengths; values near zero indicate little change in slope over the range of reach lengths. Similar to Hayakawa and Oguchi (2006), we use a threshold of one standard deviation greater than the mean relative steepness value to classify steep sections. We compare locations where relative steepness exceeds our threshold for the average SWOT and DEM profiles compared to the hydrodynamic model output.

Applying the relative steepness metric to both the hydrodynamic model and the average SWOT profile, we see that we can correctly identify many sub-kilometer high-slope sections of the Sacramento River (**Figure 7**). The hydrodynamic model for the Sacramento has an average cross sectional spacing of 258 m, so we are confident the 200 m spacing of the elevation nodes will accurately represent the river morphology. A comparison of the GPS and hydrodynamic model profiles confirms the model's accuracy (mean absolute difference: 13 cm). As such, we evaluate the steepness of the DEMs by comparing them to the hydrodynamic model. Of the 99 nodes we classified as steep in the modeled profiles, we correctly identified 74 nodes in the average SWOT profile. Twenty-five nodes were missed in the classification of the average SWOT profile, and 18 were incorrectly classified as steep. Generally, we identify all of the steep sections of the river, despite having false negatives and positives on the edges of the sections. Near 65 and 115 km flow distance we have false positives far from a true steep point, but the relative steepness of the true profile is near the threshold (**Figure 7**). We also perform this test on the water surface DEM profiles of the Sacramento River and calculate standard positive predictive values, false negative rates, and false positive rates (**Table 4**). The results of the steepness classification for the DEMs mirror the height error results in **Table 3**. SWOT outperforms the existing DEMs in all three metrics. SRTM, MERIT, and NED perform similarly, correctly identifying 20, 16, and 17 nodes. lidar captures the steep sections best of the existing DEMs, correctly identifying 42 nodes.

The relative steepness metric we use on the Sacramento River shows us that we can find the local steep sections of rivers despite the level of noise anticipated from SWOT. While slope is not a perfect predictor of physical habitat, it is a good indicator as shown by Jowett (1993). More sophisticated habitat models could be used as well, such as the in-stream habitat classification models from Demarchi et al. (2016) that use lidar and multispectral imagery. The centimeter-scale errors we report for average SWOT profiles suggest SWOT could be used as a coarse-resolution alternative to lidar for similar models, with the advantage that it will be available over nearly the entire globe.

#### Tectonics

Fault movements and the growth of associated folds can disrupt a river network in equilibrium, as rivers are sensitive to topographic changes. Large displacements can cause rerouting of river networks, whereas small displacements may result in changes to hydraulic geometry variables such as channel concavity, meander wavelength, and floodplain width as the river adjusts to external topographic gradients imposed by a regional tectonic deformation field in an effort to maintain, or reestablish an equilibrium longitudinal profile (Schumm, 1986). In large

alluvial rivers that flow over areas of known tectonic activity, we hypothesize that SWOT will be able to observe anomalies in the elevation profile of rivers potentially related to buried fault and fold structures that are unobservable in existing DEMs.

We test our hypothesis by examining the concavity of the average SWOT elevation profile for the Po River, located in northern Italy between the Apennines and southern Alps, as an example of a low-relief fluvial plain (Po Plain) with no obvious topographic expression of active tectonics, where such a signal may be recorded in spatial variations in channel elevations. The Po River flows eastward along the axis of the foreland basin formed in front of the northward propagating Apennine accretionary wedge (Castellarin, 2001). At the latitude of the river, the frontal portion of the Apennine fold and thrust belt, referred to as the Ferrara-Romagna arc in the eastern Po Plain, is buried by several kilometers of Quaternary clastic sediment and there are no obvious surface topographic expressions of the underlying faults or folds; although the structures are wellimaged by seismic refraction studies and have been intersected in hydrocarbon exploration wells (e.g., Bigi et al., 1992). The buried frontal portion of the Apennine accretionary wedge is seismically active, as exemplified by two damaging earthquakes in May of 2012, a Mw 6.1 followed 9 days later by a Mw 5.8, that were sourced on blind, shallow (<10 km) thrust faults of the Ferrara arc that project north and up-dip toward the Po River (Burrato et al., 2012).

Previous investigators have identified and interpreted planview anomalies in drainage patterns of the Po and its tributaries as evidence for the subtle perturbation of surface topography above blind thrust faults (Burrato et al., 2003), like the ones responsible for the deadly 2012 earthquakes (Burrato et al., 2012). In our analysis of the elevation of the Po River channel we first smooth the average SWOT profile using a 10 km wide windowed mean filter, and then calculate the second derivative of the elevation profile over 20 km. We find there is a spatial coincidence between the concavity of the average SWOT elevation profile of the Po River and the underlying structural geology. The highest calculated longitudinal profile concavity and convexity (negative concavity) values along the study reach are centered around 80 and 65 km flow distance, respectively. Both of these profile

TABLE 3 | Error statistic comparison of simulated SWOT profiles and profiles extracted from existing DEMs.


Profiles are compared to the hydrodynamic models used in the SWOT simulation, and GPS data where available.

inflections coincide with the intersection of the river and the outermost buried thrust faults of the Ferrara arc portion of the Apennine accretionary wedge (**Figure 8**). Away from the Ferrara arc faults, the river profile shows less extreme concavity values.

Mapping the concavity of the Po River profile, calculated from SWOT simulations, on top of the structural geology of the Po Plain suggests that SWOT data can be a useful tool in tectonic geomorphology. The high concavity section of the Po in **Figure 8** is spatially coincident with an anomalous reach (anomaly #19) TABLE 4 | Classification statistics for the relative steepness threshold.


identified by Burrato et al. (2003) as evidence for active blind thrust faulting beneath this portion of the Po Plain. Burrato et al. (2003) did not find any other anomalies in our simulated section of the Po, consistent with our interpretation of the SWOTderived concavity values. **Figure 9** shows the concavity of the available DEMs in addition to the hydrodynamic model and the average SWOT profile for the Po River. Existing DEMs show much higher concavity values than the hydrodynamic model, and with no signature of the anomalous reach identified in the average SWOT profile at a flow distance of 80 km (**Figure 8**) and in Burrato et al. (2003). Looking at **Figure 6**, the absence of a convincing concavity signature in the DEM profiles is not surprising. The vertical error in the DEMs, even at scales of tens of kilometers, obscures the effects of tectonics on the elevation profile of the Po River. Ultimately, our conclusions are based on the hydrodynamic model of the Po, but through the lens of the SWOT simulator. The features we are evaluating in this application are on the scale of tens of kilometers, so we expect the bathymetric cross section spacing of 1.2 km to accurately capture these features. The hydrodynamic model was corrupted and subsampled by the SWOT simulator in order to evaluate the potential for SWOT to see the same curvature features. Our simulated Po River profile suggests that SWOT will be able to resolve anomalous curvatures in similar rivers around the world where existing DEMs cannot.

# CONCLUSION

feart-07-00102 May 7, 2019 Time: 16:47 # 12

The CLRA method presented here greatly improves the node-level errors of river elevation profiles from SWOT by taking advantage of repeated measurements. The singular value decomposition allows us to reduce variation between profiles, largely caused by noise, while retaining some of the real variation in slope due to partial controls. Partial controls are channel features that have less effect on the surface slope with increasing discharge (Dingman, 2009). The results presented in **Figure 4** demonstrate that the CLRA profiles are able to capture changes in slope with changing stage. The elevation drop around 20 km flow distance in **Figure 4** is less pronounced when the river is at high stage, which is observable in both the simulated truth and the final CLRA profiles.

The Po River has the lowest error of all three rivers analyzed, both before and after the CLRA method. We believe the relatively low error in the simulated SWOT profiles is due to the Po being the widest river in our study, which means more simulated measurements are averaged for each node along the river centerline. While the final MAE is the lowest for the Po, the reduction in MAE by the CLRA method is greatest for the Sacramento at 70.45%. We would expect the 52 simulated overpasses of the Po to make the LRA component of the methods more effective than on the Tanana, where only 12 overpasses were simulated, but the change in MAE is approximately the same. It is reasonable to assume that some rivers can be represented with fewer eigenvectors than other rivers, but it is currently unclear what characteristics of the river, or of the SWOT noise, controls performance of the CLRA method. Future application to additional rivers will improve understanding of the algorithm performance.

SWOT is anticipated to provide WSE measurements more accurate than existing publicly available data. lidar data can provide high-resolution, high-accuracy elevation measurements over water surfaces, but is not available in most areas. Where lidar is available, it is often only at one point in time, which limits the study of many dynamic riverine processes. In comparison, SWOT will observe rivers wider than 50–100 m between 78◦

#### REFERENCES


north and south an average of twice per 21 days, or about 35 times per year (Biancamaria et al., 2016). We show that an average of SWOT observations will have accuracies better than available DEMs scaled to 100–200 m along-stream resolution. The development of a new CLRA method reduces errors in the multitemporal river profiles, without sacrificing spatial resolution like many smoothing algorithms. Creating an average SWOT profile reduces the errors even further and provides a static profile that can be used in a similar manner to existing DEMs. We show the ability of average SWOT profiles to capture changes in slope, for example, that are ecologically important for the identification of physical habitat variability for aquatic organisms. We also demonstrate that changes in river curvature identifiable in averaged SWOT profiles, may provide evidence for the subtle deformation of the Earth's surface by buried thrust faults beneath an alluvial plain, promising to provide geomorphologists a new dataset from which to decipher active tectonics.

#### AUTHOR CONTRIBUTIONS

TL and TP analyzed the data and wrote the manuscript. EA created the Tanana River hydrodynamic model and provided SWOT simulation data. AD created the Po River hydrodynamic model. RW and RF provided the SWOT simulation data of the Sacramento and Po Rivers. MD helped develop the to methods. MF advised on the physical habitats section. KW contributed to the tectonics section. All authors reviewed the final manuscript.

# FUNDING

This research was funded by Subcontract #1553389 from the SWOT Project Office at the NASA/Caltech Jet Propulsion Lab.

### ACKNOWLEDGMENTS

We thank JN and KA for their input and help improving the manuscript.

radar altimetry. J. Geophys. Res. Atmos. 107, LBA 26-1–LBA 26-21. doi: 10.1029/2001JD000609



**Disclaimer:** This work is not a product of the U.S. Government or the U.S. Environmental Protection Agency, and the author/editor/speaker is not doing this work in any governmental capacity. The views expressed are those of the authors only and do not necessarily represent those of the U.S. Government or the EPA. The findings and conclusions in this report have not been formally disseminated by the U.S. EPA and should not be construed to represent any agency determination or policy.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Langhorst, Pavelsky, Frasson, Wei, Domeneghetti, Altenau, Durand, Minear, Wegmann and Fuller. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Limitations Posed by Free DEMs in Watershed Studies: The Case of River Tanaro in Italy

#### Ricardo Tavares da Costa1,2 \*, Paolo Mazzoli<sup>1</sup> and Stefano Bagli<sup>1</sup>

<sup>1</sup> GECOsistema Srl, Cesena, Italy, <sup>2</sup> DICAM, School of Engineering, University of Bologna, Bologna, Italy

Topography is a critical element in the hydrological response of a drainage basin and its availability in the form of digital elevation models (DEMs) has advanced the modeling of hydrological and hydraulic processes. However, progress experienced in these fields may stall, as intrinsic characteristics of free DEMs may limit new findings, while at the same time new releases of free, high-accuracy, global digital terrain models are still uncertain. In this paper, the limiting nature of free DEMs is dissected in the context of hydrogeomorphology. Ten sets of terrain data are analyzed: the SRTM GL1 and GL3, HydroSHEDS, TINITALY, ASTER GDEM, EU DEM, VFP, ALOS AW3D30, MERIT and the TDX. In specific, the influence of three parameters are investigated, i.e., spatial resolution, hydrological reconditioning and vertical accuracy, on four relevant geomorphic terrain descriptors, namely the upslope contributing area, the local slope, the elevation difference and the flow path distance to the nearest stream, H and D, respectively. The Tanaro river basin in Italy is chosen as the study region and the newly released LiDAR for the Italian territory is used as benchmark to reassess vertical accuracies. In addition, the EU-Hydro photo-interpreted river network is used to compare DEM-based river networks. Most DEMs approximate well the frequency curve of elevations of the LiDAR, but this is not necessarily reflected in the representation of geomorphic features. For example, DEMs with finer spatial resolution present larger contributing areas; differences in the slope can reach 10%; between 5 m and 12 m H, none of the considered DEMs can faithfully represent the LiDAR; D presents significant variability between DEMs; and river network extraction can be problematic in flatter terrain. It is also found that the lowest mean absolute error (MAE) is given by the MERIT, 2.85 m, while the lowest root mean square error (RMSE) is given by the SRTM GL3, 4.83 m. Practical implications of choosing a DEM over another may be expected, as the limitations of any particular DEM in faithfully reproducing critical geomorphic terrain features may hinder our ability to find satisfactory answers to some pressing problems.

Keywords: digital elevation models, hydrogeomorphology, landforms, terrain descriptors, topography

# INTRODUCTION

One of the most critical elements in the hydrological response of a river basin is its topography. Among other implications, topography can significantly control the distribution of environmental variables (Sørensen and Seibert, 2007) and play a crucial role in the modeling of runoff generation and routing (e.g., Zhang and Montgomery, 1994). Its complexity can greatly influence predicted

#### Edited by:

Guy Jean-Pierre Schumann, University of Bristol, United Kingdom

#### Reviewed by:

Ahmed M. ElKenawy, Mansoura University, Egypt Zaidoon Abdulrazzaq, Independent Researcher, Baghdad, Iraq

#### \*Correspondence:

Ricardo Tavares da Costa ricardotavarescosta@gmail.com

#### Specialty section:

This article was submitted to Hydrosphere, a section of the journal Frontiers in Earth Science

Received: 31 August 2018 Accepted: 16 May 2019 Published: 04 June 2019

#### Citation:

Tavares da Costa R, Mazzoli P and Bagli S (2019) Limitations Posed by Free DEMs in Watershed Studies: The Case of River Tanaro in Italy. Front. Earth Sci. 7:141. doi: 10.3389/feart.2019.00141

**66**

discharges and how water flows over the floodplain, making it largely responsible for the accuracy of flood maps (e.g., Horritt and Bates, 2001). In fact, the great progress experienced in the modeling of both hydrological and hydraulic processes in the last decades cannot be dissociated from the advances in terrain information in the form of DEM datasets.

Digital elevation models (DEMs) are digital elevation datasets representing the Earth's surface. They may be termed more precisely as digital terrain models (DTM) or bare earth DEMs when not accounting for vegetation and buildings; otherwise, they are called digital surface models (DSM). DEMs are distributed as gridded data, where each cell, or pixel, contains a value representing the local terrain elevation. DEMs can be produced from a variety of data sources, most commonly from survey data, digitized maps, and remote sensing. Every source has its trade-offs, but remote sensing is the most versatile one, as it is able to deliver products with different areal coverages, resolutions and accuracies in an operational way – prime reason for the traction it gained in the last decades. Obviously, this popularity might not have been attained had advances in remote sensing not been on a par with those in computational power and software, as well as with the release of other important datasets, land use and land cover, for example.

Within the large variety of remote sensing techniques, we focus on the two most disruptive technologies for generating DEMs: synthetic aperture radar (SAR) and light detection and ranging (LiDAR). The greatest assets of SAR appear to be the use of long wavelengths by the active sensor, which allows imaging under all-weather conditions, and the variation of wavelengths, allowing for different degrees of reflection, e.g., by canopy or ground surface. Its major drawbacks are the presence of geometric distortions and shadows due to the side-looking antennas. On the other hand, LiDAR can be ineffective during heavy rain or in the presence of clouds, although it does not share SAR's shortcomings of geometric distortions and shadows. Both systems are capable of fast data collection, can be mounted on airborne aircrafts or satellites, can operate during day and night and can resolve the presence of canopy. They are also relatively expensive to operate and maintain (up to millions of euros). The acquisition of a LiDAR DEM can typically cost on average 140 €/km<sup>2</sup> , while a SAR acquisition can cost approximately 60 €/km<sup>2</sup> (Croneborg et al., 2015).

The open access licensing of DEMs is capable of breaking financial barriers that are frequently experienced by users, contributing to a faster advancement of science and innovation across fields. However, publicly released, freely available datasets differ in characteristics such as spatial resolution, digital terrain processing decisions and vertical accuracy, which may introduce a range of errors in the modeling of hydrological and hydraulic processes. A number of authors have assessed these errors and conjectured new ways to move forward. For example, Sanders (2007) analyzed the sensitivity of flood modeling to DEM characteristics. The author studied the Santa Clara and Buffalo Bayou rivers, US, using elevation data from LiDAR, interferometric SAR, the USGS national elevation dataset and the SRTM DEM. It was concluded that LiDAR was the best source of terrain data for that particular application and that, although very useful, the free SRTM DEM had serious limitations related to noise and data gaps. Jarihani et al. (2015) evaluated the SRTM and ASTER GDEM datasets in terms of vertical accuracy against survey marks and altimeter data, spatial resolution and digital terrain processing decisions. They demonstrated the significant impact that an underlay DEM has on flood modeling and found that the ASTER GDEM presented higher vertical accuracies in the Diamantina/Cooper river basins in Australia, while hydrologically reconditioned DEMs performed better when compared against vegetation-smoothed or unprocessed counterparts. More recently, Archer et al. (2018) compared flood modeling outcomes in a river basin in Fiji, using a commercial version of the TanDEM-X dataset (12 m spatial resolution), its vegetation-smoothed derivatives, the SRTM and the MERIT datasets against LiDAR data. The authors found that the TanDEM-X with vegetation smoothed by image classification of the amplitude map and progressive morphological filtering outperformed other datasets.

In this paper, the limiting nature of publicly released, freely available DEMs is evaluated, using LiDAR data as benchmark. However, it is done in the context of hydrogeomorphology, in other words of the study of landforms caused by the action of water, rather than focusing explicitly on flood modeling. In specific, for each DEM dataset, the upslope contributing area, the local slope, and the H and D geomorphic terrain descriptors are computed and the differences produced in terms of their cumulative frequency curves within the Tanaro river basin, in Italy, are evaluated.

The terrain descriptors analyzed are frequently used to characterize hydrological or hydraulic processes. For instance, the upslope contributing area can be associated with runoff volume, while the local slope reflects surface flow velocities (Chow, 1959), infiltration rates (Fox et al., 1997), erosional power (Knighton, 1999), drainage density (Tarboton et al., 1992), and response times (Maidment, 1993). In addition, the combination of the upslope contributing area and local slope values can be used to predict soil water content and runoff producing areas (see the topographic wetness index by Beven and Kirkby, 1979), as well as the location of channel initiation points (Montgomery and Dietrich, 1989).

The H and D terrain descriptors have also found numerous applications; for instance, Westerhoff et al. (2013) used H as topographic correction of water mapping based on SAR imagery, Nobre et al. (2016) matched a stage height to an H contour to obtain a proxy of flood extents, Elshorbagy et al. (2017) reclassified both H and D and used the product of their classes to define levels of flood hazard, Rebolho et al. (2018) and Zheng et al. (2018) used H to estimate reach-average hydraulic geometries and derive synthetic rating curves, and, finally, Clubb et al. (2017) and Nardi et al. (2019) used similar approaches to Manfreda et al. (2015) to delineate floodplains and terraces. Moreover, the terrain descriptor D can also be associated with the width function (defined as the flow path distance of any given point in a catchment to the outlet; Kirkby, 1976; Lee and Delleur, 1976), which represents a fully distributed residency time (Rodríguez-Iturbe and Rinaldo, 1997) used in the modeling of the hydrological response of a catchment. In particular, under the

assumption of constant velocity, the width function can be used to estimate a geomorphological instantaneous unit hydrograph (e.g., Mesa and Mifflin, 1986; Moussa, 2008).

In addition to investigating the intrinsic characteristics of freely accessible DEMs and given the importance of river networks in the modeling of hydrological and hydraulic processes, DEM-based river networks are delineated and compared to a photo-interpreted river network. Last, the vertical accuracy of each DEM is quantified in relation to LiDAR data.

### Free Digital Elevation Models

Open access licensing has allowed DEMs to be distributed online, free of charge to the public. This access policy not only allowed the lowering of research costs and the lifting of financial barriers to their use, but also promoted equality between individuals and institutions. As a possible indication of this, **Figure 1** plots the number of scientific publications mentioning the term "digital elevation model" per year, starting before the very first public release of a global DEM (named the Global 30 arc second Elevation Data; GTOPO30, EROS/USGS/USDOI, 1997) in 1996. It can be seen that references to "digital elevation model" have increased from less than 60 publications per year to about 800 in 2018. We note that access to global DEMs before 1996 was either restricted or inexistent and that this trend in publication records probably reflects the use of DEMs in a range of fields, amongst which the modeling of hydrological and hydraulic processes (e.g., Kumar et al., 2000). However, benefits to science and technology of publicly releasing global DEMs may be constrained, as the most popular free DEMs are now seriously dated or lack the desired spatial resolution, accuracy, correction and conditioning to keep enhancing research.

# MATERIALS AND METHODS

### Tanaro River Basin

With inception in the Ligurian Alps close to France and located in north-western Italy, the Tanaro river is the most significant right-side tributary to the Po River in terms of length (c.a. 276 km) and drainage area (c.a. 8000 km<sup>2</sup> ), presenting a highly variable discharge (Degiorgis et al., 2012). The Tanaro river basin is characterized by steep mountainous terrain and a nearly

flat alluvial region (**Figure 2**). The river itself is highly prone to flooding; indicatively, during the 1994 historical Piedmont flood and landslide, 44 persons lost their lives, 2000 were displaced and a whopping 8.8 billion € in damages were estimated (Luino, 1999).

The Tanaro river basin was chosen as case study due to its peculiar characteristics and history of disastrous events. In this work, the DEMs listed in **Table 1** were clipped with the Tanaro river catchment polygon obtained from the Italian Environmental Agency (ISPRA – Istituto Superiore per la

feart-07-00141 June 1, 2019 Time: 10:30 # 5


TABLE 1|Digital elevation models (DEM) currently available free of charge with spatial resolutions below 3 arc seconds.

Datasets are usedas provided, except for the coordinate systems that were transformed to WGS84/EDM96 when not already referred to this reference system and a transformation grid was easily available.

Protezione e la Ricerca Ambientale) in shapefile format. The clipped DEMs are used to extract the terrain descriptors within the study area, namely (1) the upslope contributing area, (2) the local slope, (3) the flow path elevation difference to the

nearest stream, H, and (4) the flow path distance to the nearest stream, D, used for comparison in terms of cumulative frequency curves. In addition, river networks are delineated from the clipped DEMs for visual inspection and comparison with the EU-Hydro photo-interpreted river network (EEA, 2015b). As a final step, vertical accuracies of the free DEMs are reassessed using LiDAR data.

# Description of DEM Datasets

In **Table 1**, an overview of some of the most common DEMs with spatial resolutions of the order of 3 arc seconds (c.a. 90 m) or less is provided. These datasets are currently in the public domain or available upon request mostly for research or other non-commercial purposes. **Table 1** is organized by ascending order of year of public release. In this work, all datasets in **Table 1** are taken into consideration:


# LiDAR Dataset for Italy (Benchmark)

As benchmark for the assessment of vertical accuracies, a LiDAR dataset that partially covers the Tanaro river basin (footprint in **Figure 2**) is used and was obtained from the Italian Ministry of Environment, Land and Sea (Ministero dell'Ambiente e della Tutela del Territorio e del Mare). The LiDAR data was resampled to the corresponding spatial resolution of each DEM in **Table 1**. The LiDAR dataset (with spatial reference WGS84/ITALGEO95) is available to the general public upon formal request and upon payment of a processing fee (to visualize its areal coverage, please visit the Italian National Geoportal – Geoportale Nazionale).

#### Terrain Descriptors

The extraction of the selected terrain descriptors from the free DEMs follows a simple workflow using the TauDEM toolbox (see **Figure 3**; Tarboton, 2015). A clipped DEM is first corrected by identifying sinks and by raising cell elevation values to the level of the lowest pour point in the eight surrounding cells of the structured grid. This is deemed necessary in order to avoid interference with flow routing. From this corrected layer, flow directions from each cell to one of its eight neighbors are determined by following the steepest descent (also known as convergent eight direction flow model, abbreviated as D8 flow directions) and a counter-clockwise coding from 1 (flow to the East) to 8. Using the D8 flow model, the local slope or tangent of the angle of incline, θ, is calculated as the drop, 1y, over distance, 1x, between a cell and its neighbors in the flow path:

$$\tan\left(\theta\right) = \frac{\Delta y}{\Delta x} \tag{1}$$

In turn, the upslope contributing area is obtained by simply accumulating cells following the D8 flow directions.

To delineate river networks from the free DEMs, channel heads are first identified by imposing a threshold of 10<sup>5</sup> m<sup>2</sup> (Giannoni et al., 2005) on an area-slope criterion that characterizes the transition between transport mechanisms (Montgomery and Dietrich, 1988, 1989). This criterion is defined as the product of upslope contributing area and local slope

raised to an exponent k responsible for drainage density changes, assuming k = 1.75 throughout this study (Giannoni et al., 2005). Starting at the channel heads, river networks are delineated following the D8 flow directions to the outlet.

With the river network delineated from the free DEMs, the computation of H and D is programmed in Python following the description in Manfreda et al. (2015) and making use of the Geospatial Data Abstraction Library (GDAL) for raster I/O. H is obtained by calculating the elevation difference between each cell in the DEM raster and the connected river network cell, following the D8 flow directions. D is obtained by counting the number of cells from each position in the DEM raster to the connected river network cell, still following the D8 flow directions. Furthermore, for each unique flow path, adjacent cell counts need to be distinguished from diagonal ones, so that lengths can be obtained by multiplication with the spatial resolution or with the product of spatial resolution and <sup>√</sup> 2, respectively. The total flow path distance from each location in the raster to the stream is simply the sum of corresponding adjacent and diagonal lengths.

#### Accuracy Assessment

In order to assess the vertical accuracy of the free DEMs, three common error measures were selected to be used with continuous variables, in this case the elevation data. Different error measures are report in this study as they may complement each other (Chai and Draxler, 2014). The systematic error or statistical bias is defined as the simple difference between DEM, byi , and LiDAR, y<sup>i</sup> , elevations (here assumed as the truth):

$$\text{BIAS} = \widehat{\wp\_i} - \wp\_i \tag{2}$$

where i the index of an individual cell in a flattened raster. The MAE is defined as:

$$\text{MAE} = \frac{1}{n} \sum\_{i=1}^{n} \left| \widehat{\boldsymbol{\wp}}\_{i} - \boldsymbol{\wp}\_{i} \right| \tag{3}$$

where n is the total number of cells. The MAE represents the average absolute difference between DEM and LiDAR elevations and gives an indication of the magnitude of error. Finally, the RMSE is defined as:

$$\text{RMSE} = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left(\widehat{\nu\_i} - \wp\_i\right)^2} \tag{4}$$

where the mean square error is the second moment of the bias. The RMSE also measures the magnitude of error, but with a higher sensitivity to outliers, thus putting stronger emphasis to unfavorable conditions (Chai and Draxler, 2014). Its normalized version that is less sensitive to outliers is given by:

$$\text{NRMSE} = 100^{\*} \frac{\text{RMSE}}{\nu\_{\text{max}} - \nu\_{\text{min}}} \tag{5}$$

Finally, the linear correlation between <sup>b</sup>y<sup>i</sup> and <sup>y</sup><sup>i</sup> is also reported and measured using the Pearson's correlation coefficient (PC) defined as:

$$\text{PC} = \frac{\sum\_{i=1}^{n} \left(\boldsymbol{\wp}\_{i} - \bar{\boldsymbol{\wp}}\right) \left(\widehat{\boldsymbol{\wp}}\_{i} - \overline{\widehat{\boldsymbol{\wp}}}\right)}{\sqrt{\sum\_{i=1}^{n} \left(\boldsymbol{\wp}\_{i} - \bar{\boldsymbol{\wp}}\right)^{2} \sum\_{i=1}^{n} \left(\widehat{\boldsymbol{\wp}}\_{i} - \overline{\widehat{\boldsymbol{\wp}}}\right)^{2}}} \tag{6}$$

PC takes values between −1 and 1, with PC = −1 corresponding to a perfect inverse correlation, PC = 1 corresponding to a perfect direct correlation and PC = 0 corresponding to no linear correlation.

#### RESULTS

In this section, the outcomes of the method proposed above are examined. The cumulative frequency curves of terrain descriptors presented in **Figure 4** show important differences within the Tanaro river basin and the reader can also refer to **Figure 5** for the average values of the terrain descriptors. Note that the upslope contributing area and D have not been computed for the LiDAR data, as both indicators proved to be meaningless within the limited extent of the LiDAR footprint.

The cumulative frequency curves of elevation within the LiDAR footprint (**Figure 4A**) show that a large number of free DEMs approximate well the curve obtained from LiDAR data at 10 m spatial resolution. Namely, the EU DEM, the VFP DEM, the SRTM GL1 and GL3 and the MERIT DEM provide the best approximation and are very closely followed by the ASTER GDEM and the AW3D30 DEM. Minor differences may be observed between these datasets at lower elevations. In the lower part of the LiDAR footprint the HydroSHEDS hydrologically conditioned DEM has a different cumulative frequency curve than the previous mentioned datasets; in particular, below a certain elevation value, a slightly higher frequency may be expected for the HydroSHEDS DEM in comparison to the LiDAR (or lower elevation values for a certain frequency). Nevertheless, for the remaining 80% of the LiDAR footprint, the HydroSHEDS DEM provides a reasonable approximation of the LiDAR frequency curve. The TDX and the TINITALY DEMs present surprisingly similar cumulative frequency curves between themselves but, at the same time, significantly different from the LiDAR data and the remaining free DEMs; in specific, below a certain elevation value, a lower frequency may be expected for the TDX and the TINITALY DEMs (or higher elevation values for a certain frequency).

In terms of the upslope contributing area, the frequency distribution within the entire Tanaro river basin (and therefore not merely within the LiDAR footprint as before), **Figure 4B**, shows that all curves start to converge at around 150 accumulated cells (i.e., areas up to 1.2 km<sup>2</sup> ) and that 90% of the cells have contributing areas below such value. Within the remaining upper 10%, DEMs with finer spatial resolution present higher upslope contributing areas than DEMs with coarser spatial resolution, as more cells are accumulated downstream in the former.

For the cumulative frequency curves of the local slope within the LiDAR footprint (**Figure 4C**), it is possible to observe that for slopes steeper than 22% (c.a. 40◦ ) the TDX DEM gives the best overall approximation of the LiDAR data, followed by the SRTM GL1, the EU DEM and the ASTER GDEM. In the same range, the MERIT and the VFP DEMs overlap and follow closely the cumulative frequency curve of the HydroSHEDS and the SRTM GL3 DEMs, with higher frequencies relative to the LiDAR; while the remaining DEMs present lower frequencies relative to the LiDAR. For slopes values below 22%, every DEM curve presents a lower frequency relative to the LiDAR. Of all datasets, the AW3D30 and the TINITALY DEMs present the least representative approximation of the LiDAR curve, with lower frequencies. In general, all datasets have more than 90% of the cells with slope falling below 45◦ , approximately. Instead, if one looks at the average values of local slopes (**Figure 5**) within the LiDAR footprint the best approximation is given by the HydroSHEDS, followed by the SRTM GL3, the VFP and the MERIT DEM.

**Figure 4D** shows the cumulative frequency curve of H. In the case of the DEMs, the curve is S-shaped with a sharp increase from 50 (c.a. 10% frequency) to 150 m (c.a. 90% frequency), approximately. Contrarily, the LiDAR curve does not present a perfect S-shape; instead, there is a bump around 12 m. The best approximation of the curve from LiDAR data is given by the ASTER GDEM with a perfect match at low H values, followed by the SRTM GL1 and the TDX DEM with very close representation at high H values. Nevertheless, between 5 m and 12 m not a single DEM can faithfully represent the curve obtained from LiDAR data. In turn, the AW3D30 and the EU DEMs have similar curves that, for the same H value, consistently present a higher frequency relative to the MERIT, the SRTM GL3, the VFP (all three of them overlapping in the graph) and the HydroSHEDS DEM. The ASTER GDEM, the SRTM GL1 and the TDX DEM fall more or less between the previous cases, while the TINITALY, for a same H value, presents the lowest frequency of all DEMs and has the curve that is furthest apart from that of the LiDAR. In general, all datasets have more than 90% of the cells with H falling below 150 m, approximately. By looking at the average values of H (**Figure 5**), no single DEM within the LiDAR footprint can approximate the LiDAR value of 59 m, the closest one being the AW3D30 DEM, but still with a difference of 28 m.

The cumulative frequency curve of D is also S-shaped, with a sharp increase after 400 m (c.a. 10% frequency) and up to 5 km (c.a. 90% frequency), approximately. By looking at how the DEMs compare between themselves, it was found that the AW3D30 DEM and the ASTER GDEM have frequency curves that consistently present a higher frequency for the same D value relative to the MERIT, the SRTM GL3, the VFP (all three overlapping) and the HydroSHEDS DEM. The SRTM GL1 and the TDX DEM fall more or less in between the curves of the remaining free DEMs.

In **Figure 5**, we look into sample regions of the Tanaro river basin in Italy, where the DEM-derived river networks are overlaid. To the side of each sample region, the corresponding portion of the EU-Hydro photo-interpreted river network is presented for visual comparison. In the first sample region

FIGURE 5 | Overlay of river networks for three distinct regions in the Tanaro river basin, Italy, derived from free digital elevation models (DEM). For each region, the corresponding EU-Hydro photo-interpreted river network is shown on the right. (A,B) are sample regions representative of flatter terrain; and (C) is a sample region representative of mountainous terrain. On the bottom left, the locations of the three regions within the Tanaro river basin are marked with red boxes. On the bottom right, the average values of the terrain descriptors are presented.

(**Figure 5A**), located in a flat area, it is clear that no single dataset can faithfully reproduce the photo-interpreted river network. Although all DEM-derived river networks are able to represent the river confluence in the region, they generally fail to match the location: they appear further upstream compared to the photo-interpreted river network. Furthermore, there is a tendency of a few datasets, namely the ASTER GDEM, the AW3D30 and the SRTM GL1, to create what seem to be spurious tributaries in the region.

In the second sample region (**Figure 5B**), another flat area, the only two datasets capable of representing the EU-Hydro are the HydroSHEDS and the EU DEM. The VFP, which once again matches perfectly the SRTM GL3, as well as the SRTM GL1 and the MERIT DEM, fails to properly locate the river confluence (further upstream compared to the EU-Hydro). The same happens with the TDX DEM, but with a more unrealistic meandering. Last, the AW3D30, the TINITALY DEM and the ASTER GDEM completely fail to represent the river network in the sample region, displaying a significant offset from the photo-interpreted river network, unrealistic meandering, unrealistic placement of the river confluence and even no confluence at all in the case of the AW3D30 DEM.

In the third sample region (**Figure 5C**), located in a mountainous area, all datasets are capable of a faithful representation of the photo-interpreted river network and specifically of the river confluence. However, a number of datasets, namely the TINITALY, the AW3D30, the TDX,

FIGURE 6 | (A) Map of statistical bias between the SRTM GL3 digital elevation model (DEM) and the light detection and ranging (LiDAR) data in the Tanaro river basin, Italy; (B) map of statistical bias between the HydroSHEDS DEM and the LiDAR data; (C) DEM elevations plotted against LiDAR elevations; (D) cumulative frequency of absolute BIAS plotted against the elevation difference to the nearest stream, H.

the EU DEM and the ASTER GDEM, present possibly spurious tributaries in the sample region.

In **Figure 6A**, the spatial distribution of bias for the SRTM GL3 is presented. In this figure there appears to be a trend toward overestimating elevation in mountainous areas and a trend toward underestimating elevation in flatter terrain. All other DEM datasets follow more or less this same pattern, with the exception of HydroSHEDs DEM (**Figure 6B**), where there is a pronounced overall trend for underestimation, with a more limited overestimation in mountainous regions. **Figure 6C** shows that all free DEMs are linearly correlated with the LiDAR data, even though different degrees of dispersion can be noted, particularly in the TINITALY, HydroSHEDS, and the AW3D30 DEM. **Figure 6D** shows the cumulative bias as a function of H, which should help to differentiate between flatter and mountainous terrain, or floodplains and hillslopes, at any given elevation. It becomes clear that there is a lower bias in flatter terrain, where small H values are found, and higher bias in hillslopes. For all DEM datasets, the cumulative bias appears to be similar in the floodplains (c.a. 30 m H), with the VFP and SRTM GL3 presenting the smallest value and the AW3D30 DEM the highest one, while the remaining DEMs present values that fall in between. In hillslopes (above c.a. 30 m H), a similar behavior can be found, with the difference that the SRTM GL3 DEM accuracy tends to degrade faster than that of the VFP DEM. The SRTM GL3 DEM cumulative bias is matched by the MERIT DEM at c.a. 100 m H and by the HydroSHEDS DEM at c.a. 300 m H. The ASTER GDEM, the EU DEM and the TDX DEM share similar curves. Similarly, the TINITALY DEM, the SRTM GL1 and the HydroSHEDS DEM have curves that nearly overlap.

In **Table 2**, results in terms of MAE and RMSE are presented. The lowest MAE is presented by the MERIT DEM, followed closely by TINITALY, the SRTM GL3, the SRTM GL1, the AW3D30, TDX, VFP and the EU-DEM. The ASTER GDEM and the HydroSHEDS DEM present the highest MAE among all datasets, more than double the MAE of the MERIT DEM. In terms of RMSE (and NRMSE), the lowest value is presented by the SRTM GL3, followed closely by the MERIT, SRTM GL1, TINITALY and TDX. The EU DEM, the VFP and the AW3D30 DEMs are found to have RMSE values with c.a. 2.5 m more than

TABLE 2 | Vertical accuracy assessment of the free digital elevation models (DEM) in the Tanaro river basin, Italy, expressed as mean absolute error (MAE), root mean squared error (RMSE), normalized RMSE (NRMSE), and Pearson correlation (PC), with the light detection and ranging data (LiDAR) data used as benchmark.


that of the SRTM GL3. The ASTER GDEM and the HydroSHEDS DEM, as with the MAE, present the highest value among all, more than double the RMSE of the SRTM GL3. Discrepancies found among datasets are explained to a limited degree by the lack of transformation of the LiDAR data and of the EU DEM to the EGM96 model, from the national and regional geoids, respectively. Unfortunately, an easily accessible transformation grid was not available to perform this step.

# DISCUSSION AND CONCLUSION

In this analysis, the use of different DEMs is shown to result in some noteworthy discrepancies in the extracted terrain descriptors. In fact, the ability of each free DEM to capture a particular erosional or depositional landform caused by the action of water may vary with spatial resolution, hydrological reconditioning and vertical accuracy. This can impact any subsequent hydrogeomorphic analysis that requires a faithful representation of geomorphic features or the modeling of hydrological and hydraulic processes.

By looking at the cumulative frequency curves of elevation, it is found that the EU DEM, the VFP, the SRTM GL1 and GL3, the MERIT, the ASTER GDEM and the AW3D30 DEMs approximate the LIDAR data curve well. However, this is not necessary or sufficient to assert a faithful representation of geomorphic features. On the other hand, it is also observed that the HydroSHEDS DEM is significantly different from the LiDAR data at lower elevations and that the TDX and TINITALY DEMs is consistently different from LiDAR throughout the entire range of elevations.

In the case of the upslope contributing area within the whole Tanaro river basin, it is found that for areas below 1 km<sup>2</sup> (c.a. 90% of the basin) DEMs with a finer spatial resolution seem to present lower area values than DEMs with a coarser spatial resolution, while for areas above 1.2 km<sup>2</sup> (c.a. 10% of the basin) DEMs with finer spatial resolutions seem to present larger values than DEMs with coarser ones. The same conclusion can also be drawn by looking at the average values, where the TINITALY DEM, with a spatial resolution of 10 m, presents the highest upslope contributing area. The smaller fraction of cells with higher upslope contributing areas can generally be associated with cells belonging to the river network, and channel initiation may occur further upstream in DEMs of finer resolution.

The local slope differs from the LiDAR data to some degree, except in the case of HydroSHEDS, followed closely by the SRTM GL3, the VFP and the MERIT DEM. The best approximations are given by the TDX, the SRTM GL1, the EU DEM and the ASTER GDEM. The average values of local slope show that differences can reach about 10%.

The cumulative frequency curve of H is characterized by a sharp increase around 5 m, followed by a slower increase around 12 m and finally by another sharp increase from 30 m upward. DEMs with spatial resolution of 30 m, in addition to the TDX DEM, give the best approximation of the LiDAR H curve. In particular, the ASTER GDEM perfectly matches the LiDAR cumulative frequency curve at low H values. TINITALY,

followed closely by HydroSHEDS at lower H values, presents a significantly different curve from LiDAR, while average values of H are significantly different between every DEM and the LiDAR.

The cumulative frequency curve of D is characterized by a sharp increase between 400 m and 5 km. The AW3D30 and the ASTER GDEM present consistently higher cell counts per D contour with respect to the remaining DEMs, while HydroSHEDS and the MERIT DEM present consistently lower cell counts; the SRTM GL1 and the TDX DEM fall more or less between all the other curves.

Regarding the DEM-based river network, it is confirmed that in the Tanaro river basin the river network extraction can be more problematic in flatter terrain. In particular, the location of river confluences and meanders can be significantly misrepresented. Over flat areas, coarser spatial resolution DEMs tend to better approximate the photo-interpreted river network, while hydrologically reconditioning a DEM (e.g., the HydroSHEDS) also seems to help.

In terms of BIAS, a tendency to overestimate elevation values in hillslopes and a tendency to underestimate elevations in floodplains appear to exist, except for the HydroSHEDS DEM that tends to underestimate elevations in a more generalized way. In spite of the BIAS found – less in the case of the SRTM GL3 and VFP DEMs – it is shown that DEMs are highly correlated to the LiDAR data, with the TINITALY, HydroSHEDS and AW3D30 showing some noticeable dispersion.

Finally, vertical accuracy measures were computed from each DEM. The lowest MAE has been obtained by the MERIT DEM, 2.85 m, while the ASTER GDEM and the HydroSHEDS DEM have presented the highest MAE of all the datasets, 11.16 and 6.61 m, respectively. The lowest RMSE, has been presented by the SRTM GL3, 4.83 m, or 0.21% NRMSE, while the ASTER GDEM and the HydroSHEDS DEM have presented, once again, the highest RMSE values among all datasets, 15.78 and 9.39 m, or 0.40 and 0.68% NRMSE, respectively.

In practice, differences found may affect several aspects of the modeling of hydrological and hydraulic processes, resulting in diverse outcomes. For instance, differences in upslope contributing area may influence scaling regimes that depend directly on it (Dodov and Foufoula-Georgiou, 2005) or may have an impact on the faithful representation of a river network and of channel hydraulic geometry. The slope may affect channel initiation and the flow of surface water. H may significantly influence the outcomes in low-complexity flood modeling (e.g., Rebolho et al., 2018; Zheng et al., 2018), flood mapping (Degiorgis et al., 2012; Manfreda et al., 2014) and flood detection (e.g., Westerhoff et al., 2013), while D may impact hydrological modeling (e.g., Mesa and Mifflin, 1986; Moussa, 2008).

Concerning the Tanaro river basin, finer spatial resolutions have not always improved the representation of the morphology. With respect to the hydrological reconditioning of DEMs, and even though only the HydroSHEDS DEM was given as an example, the procedure has appeared helpful in achieving a better definition of the river network, but it might also hinder the representation of fluvial landforms, something evidenced by Callow et al. (2007). In the case of vertical accuracy, elevation-based terrain features were affected by the magnitude of the vertical error. Indeed, error metrics can be misleading, since low values may hide significant problems in the representation of the terrain in different parts of the basin, as shown by the BIAS.

Ideally, a good DEM should aim for an optimal balance between spatial resolution, vertical accuracy and DEM reconditioning. Considering the numerous applications that free DEMs have and the fast evolution of instruments and techniques experienced over the past years, one might expect that new and improved datasets, as well as their derivatives, would be released more regularly. The lack of systematic updating constitutes a major hurdle. Operational flood forecasting systems that authorities rely on for early warning and which are as good as their margin of error, as well as flood maps, are only a few examples of applications that may suffer from the lack of timely DEM updates. It is shown that there is still substantial work to be done an there is some indication that a new generation of free DEMs can form the foundations on which to sustain and even accelerate progress in several fields. Therefore, we would like to use this opportunity not only to call for high-accuracy, open-access global DEMs, but also for the systematic processing of new and existing datasets, based on state-of-the-art techniques, with an exhaustive error quantification, the inclusion of comprehensive documentation and the production of terrain derivatives that can prove indispensable for saving valuable resources. Last, it is also recommended the fusion of existing free DEMs and LiDAR data, where the latter is available (e.g., Italy and United Kingdom), as this process may enhance the quality of the digital terrain representation.

# AUTHOR CONTRIBUTIONS

RT conceptualized the work, carried out the analysis, and wrote the manuscript. PM and SB reviewed and provided the comments.

# FUNDING

This work was supported by the European Union's EU Framework Programme for Research and Innovation Horizon 2020 under Grant Agreement No. 676027.

# ACKNOWLEDGMENTS

We would like to thank the European Union's Horizon 2020 Research and Innovation Program under the Marie Skłodowska-Curie System-Risk ETN (https://www.system-risk.eu/) grant, the European Union's Earth Observation Programme Copernicus, NASA for the release of the SRTM DEMs, INGV for TINITALY, the USGS for the HydroSHEDS dataset, NASA/METI for the release of ASTER GDEM, the European Union's Copernicus Programme for the EU-DEM and EU-Hydro, Mr. de Ferranti for the VFP DEM, JAXA for the AW3D30 DEM, Dr. Dai Yamazaki and the University of Tokyo for the MERIT DEM, DLR for the

TDX DEM, and the Italian Ministry of Environment, Land and Sea for the LiDAR dataset for Italy. We would also like to thank Dr. David Tarboton (Utah State University) for TauDEM

#### REFERENCES


utilities, the Python Software Foundation for the Python programming language, the Open Source Geospatial Foundation for the GDAL.


National Aeronautics and Space Administration (NASA) Jet Propulsion Laboratory and Ministry of Economy, Trade, and Industry (METI) of Japan [Dataset]. Available at: https://search.earthdata.nasa.gov/search (accessed July, 2018).


**Disclaimer:** This publication reflects only the authors' views and the European Union is not liable for any use that may be made of the information contained therein.

**Conflict of Interest Statement:** RT is employed by company GECOsistema Srl. PM and SB are owners of the company GECOsistema Srl and declare no competing interests.

Copyright © 2019 Tavares da Costa, Mazzoli and Bagli. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Bare Earth DEM Generation for Large Floodplains Using Image Classification in High-Resolution Single-Pass InSAR

#### Drew Faherty1,2 \*, Guy J.-P. Schumann2,3 and Delwyn K. Moller<sup>2</sup>

<sup>1</sup> Department of Geological Sciences, California State Polytechnic University, Pomona, CA, United States, <sup>2</sup> Remote Sensing Solutions, Barnstable, MA, United States, <sup>3</sup> School of Geographical Sciences, University of Bristol, Bristol, United Kingdom

The increasing frequency and severity of flood events demand improved accuracy of hydrodynamic models to better mitigate the societal and economic consequences of these disasters. Ideally, these models derive elevation data from high-accuracy LiDAR, which often captures the subtle variations in elevation and microtopography – elements that are critical to high-resolution hydrodynamic modeling. Due largely to the cost of acquisitions, these data, however, are not widely available at high spatial resolutions for large-scale areas, i.e., floodplains and coastal regions, especially outside of the United States and Western Europe, where flood-prone communities already lack the infrastructure and resources to manage these disasters. High-resolution interferometric synthetic aperture radar (InSAR) systems may offer a viable and cost-effective alternative to LiDAR, with comparable spatial resolutions and vertical accuracies. Like any swath mapping sensor, systematic errors in calibration knowledge increase as the leverarm increases when viewing further off vertical. InSAR-derived digital elevation models (DEMs) are often corrected using LiDAR or ground control points, although this limits the application of InSAR to those areas where these data exist. In this paper, we present a novel approach for using near-range InSAR data to correct inherent systematic bias that propagates across adjacent, overlapping swaths for the same airborne acquisition. This eliminates the need for LiDAR in generating InSAR DEMs, while maintaining a vertical error of approximately 0.26 m relative to LiDAR, at a spatial resolution of 30 m. Data was acquired using the NASA/JPL airborne, single-pass, Ka-band GLISTIN-A (Glacier and Ice Surface Topography Interferometer) InSAR over the Red River Basin (RRB), ND. The accuracy of the final DEM is evaluated using National Geodetic Survey (NGS) GPS and a high-resolution LiDAR DEM.

Keywords: digital elevation model, floodplain mapping, InSAR, remote sensing, flood hazard

# INTRODUCTION

GLISTIN-A (Airborne Glacier and Land Ice Surface Topography Interferometer) is a Ka-band (35.66 GHz), single-pass interferometric synthetic aperture radar (SPInSAR), inspired by the Shuttle Radar Topography Mission (SRTM) heritage (Rabus et al., 2003) and developed to provide all-weather, high-resolution swath ice surface topography (Moller et al., 2011) not available through

#### Edited by:

Nick Van De Giesen, Delft University of Technology, Netherlands

#### Reviewed by:

Ioannis N. Daliakopoulos, Technological Educational Institute of Crete, Greece Lorenzo Iannini, Delft University of Technology, Netherlands

> \*Correspondence: Drew Faherty dfaherty88@gmail.com

#### Specialty section:

This article was submitted to Hydrosphere, a section of the journal Frontiers in Earth Science

Received: 08 January 2019 Accepted: 27 January 2020 Published: 20 February 2020

#### Citation:

Faherty D, Schumann GJ-P and Moller DK (2020) Bare Earth DEM Generation for Large Floodplains Using Image Classification in High-Resolution Single-Pass InSAR. Front. Earth Sci. 8:27. doi: 10.3389/feart.2020.00027

**81**

existing spaceborne LiDAR or radar sensors. During engineering flights, GLISTIN-A also mapped parts of the Central Valley region of California and it became clear that the sensor's capability to map land heights and water areas should be exploited. To this end, Schumann et al. (2016) undertook a feasibility study to examine this in detail. Indeed, their DEM case study over the California Central Valley confirmed for the first time that GLISTIN-A's Ka-band technology can be used for floodplain mapping and flood studies on small to regional scales, as it generates floodplain DEMs with similar accuracies to that of state-of-the-art airborne LiDAR (mean bias of 5.6 cm, with a standard deviation of ±30 cm), therefore, furthering the international agenda of a high-accuracy, open-access global DEM (Schumann et al., 2014).

Digital elevation models (DEMs) are essential data sets for disaster risk management and humanitarian relief services, as well as many environmental process models. While in some countries, such as the United States, United Kingdom, Australia, and several others, there are available LiDAR DEMs, free global DEMs do not currently meet basic hydrodynamic requirements (Schumann and Bates, 2018). Thus, in many regions of the world, where flooding and other hazards are a major concern, we are limited in modeling and mapping primarily due to the lack of LiDAR-quality DEMs, which in turn limits our ability to unlock the full potential of these models in support of flood applications. As demonstrated by Schumann et al. (2016), however, the GLISTIN-A NASA airborne mission offers a novel single-pass InSAR technology that is cost-effective and can generate floodplain DEMs with a RMSE of approximately 0.3 m for an average spatial resolution of 30 m. While Schumann et al. (2016) study demonstrated the capability of GLISTIN-A to generate a LiDAR-quality DEM, a more targeted flood use case is required to evaluate how this DEM performs when used to condition a 2-D hydrodynamic model and simulate an actual out-of-bank flood event. In their study, Schumann et al. (2016) stated in conclusion that "in order to more fully assess the suitability of the Ka-band InSAR instrument to acquire high-accuracy land elevation data, we suggest to repeat this type of analysis over a more natural floodplain and a larger river, and for an area more prone to regular flooding." This study addresses this type of analysis by acquiring land elevation data over the Red River of the North wide-area floodplain using the GLISTIN-A airborne SPInSAR.

# MATERIALS AND METHODS

#### Test Site

As a test location, we chose the Red River of the North (RRN) basin (**Figure 1**), located in North Dakota and Minnesota in the United States, and in southern Manitoba, Canada. The Red River of the North flows north through Winnipeg, Manitoba, and is a tributary of the Nelson River basin, which carries runoff into Hudson Bay from much of southern Canada, from the Great Lakes to the Continental Divide. Due to its extremely low relief, the region of the RRN depicted in **Figure 1** is prone to flooding.

According to NOAA's Spring Outlook released mid-March 2017, Northern North Dakota (the Souris River, Devils Lake, and the northernmost reaches of the Red River) was predicted to have the greatest risk of major flooding last spring. Snowpack was heavy in the northern plains (containing up to four inches of liquid water that could have increased with additional storms through April), and if long term warm-up had coincided with spring rains, already-saturated soils would not have been able to absorb the increased water, which may have led to increased runoff and potential flooding. Apart from the flood-prone nature of this region, the choice of this test site was also governed by the fact that (i) GLISTIN-A on its March/April 2016 Greenland mission presented a unique overpass opportunity that we leveraged in full by acquiring InSAR data over the area depicted in **Figure 1** and (ii) flood hazard prediction during the snow-melt season in this region is a known challenge for the NOAA National Weather Service (NWS) North Central River Forecast Center (see e.g., Tuttle et al., 2017).

#### GLISTIN-A InSAR Data

InSAR is able to retrieve surface topography by displacing two antennas in the cross-track direction to view the same surface. The interferometric combination of data, received on the two antennas, allows one to resolve the path-length difference from the illuminated area to a fraction of a wavelength. From the interferometric phase, the height of the target (or phase center) can be estimated; therefore, an InSAR system such as GLISTIN-A is capable of providing not only the position of each image point in along-track and slant-range, as with a traditional SAR, but also the height of that point through the use of the interferometric phase (Rodriguez and Martin, 1992; Rosen et al., 2000).

GLISTIN-A (**Table 1**) is a Ka-Band (8 mm wavelength), cross-track, single-pass InSAR, developed for high-precision, high-resolution ice-surface topography mapping, with a swathwidth greater than 10 km (Moller et al., 2017). Using standard processing, the imagery maintains a spatial posting of 3 m. The short-wavelength minimizes interferometric penetration of the electromagnetic wave into surface media. While airborne laser altimetry has negligible penetration into surface media, it is limited in swath width (up to 500 m) and therefore, in spatial coverage.

# SAR Data Calibration

Instrument calibration and the InSAR processing was performed at the Jet Propulsion Laboratory (JPL). Range pixel location can be determined to better than a tenth of a pixel by oversampling the slant-plane imagery. Because range measurements to the two interferometric channels may differ, a differential range correction is computed, by measuring range offsets between the two channels. Differential range measurements are accurate to better than a hundredth of a range pixel and insure proper range registration of the channels during interferogram formation. After determining the common and differential range corrections, the data are reprocessed, and strip map DEMs and orthorectified imagery are generated.

However, once the initial instrument calibration is complete, residual imperfect knowledge or instrumental drift can affect the

interferometric height measurement and introduce a systematic trend in the data. Rodriguez and Martin (1992) provide a generalized summary of such error sources for a single-pass interferometer and their impact of the topographic measurement. To first order the impact of calibration uncertainty can be lumped into three categories:


Moller et al. (2019) predicts that quadratic height error sources (number 3 above) are negligible based on system design and measurement and this is borne out by a comprehensive validation of 2 years of GLISTIN-A data (2016 and 2017) over coastal Greenland. Validation with overlapping lidar data verify that residual systematic trends manifest as a small linear height characteristic propagating across the swath. Specifically, Moller et al. (2019) observed up to meter-scale nadir biases and evenly distributed residual slopes with a standard deviation of ∼8 millidegrees or 0.18 mm/m).

**Figure 2** shows the interferometric geometry where for simplicity the baseline, B is horizontal, the platform is at altitude h<sup>p</sup> above a reference surface, and it is imaging a point of height h<sup>0</sup> above that reference, located at incidence angle θ0, and at slant range ρ. The cross-track distance of that point is therefore x<sup>0</sup> = ρsin(θ0). It is well known that for (non-interferometric) SAR imagery terrain features can be displaced across-track due to layover or foreshortening, but the additional information provided using interferometry allows features to be correctly geolocated. However, any errors in the calibration knowledge will introduce an error in reconstructing the surface location of the interfered return. **Figure 2** illustrates this whereby the rotated dashed triangle shows how a baseline angle error, 1θ, will introduce a positional error in both the cross-track, 1x, and height. 1h, of the reconstructed scene. With perfect system knowledge (the solid triangle):

$$h\_t = \frac{\rho}{\cos(\theta\_0)} = \frac{x\_0}{\tan(\theta\_0)}.$$

TABLE 1 | Summary of key GLISTIN-A operating parameters.


But with residual errors manifesting as a baseline angle error estimate, 1θ, the resulting height estimate is:

$$h\_t' = \frac{\varkappa\_0 + \Delta x}{\tan\left(\theta + \Delta \theta\right)}$$

.

The vertical and horizontal geolocation error is a result of imperfect geometrical knowledge. This illustrates the mechanism by which a slope is introduced: that is a systematic height change that is a function of cross-track distance, x, and where the slope 1h/1x = ht-h 0 t /1x. For generating a DEM of large floodplain areas, the impact of 1x can be considered negligible at high incidence angles and tolerable for small incidence angles.

For the RRN acquisitions, the flight plan was configured to have near/far overlaps in successive passes. **Figure 3** shows the iterative methodology applied to generate the area DEM. A linear trend, or an effective tilt of the swath can be characterized with two parameters: a constant offset and a slope. To enforce consistency without external control we iteratively correct for slope only by forcing the overlapping extent of the far-range of one swath to equal the near range of another. This is done by correcting the slope difference for the overlapping extent (∼2 km across track overlap) calculated as the mean for the entire line extent (>50 km). The resultant domain may have a mean vertical offset with respect to the vertical reference datum but represents a self-consistent digital surface model (DSM) with which to conduct image classification and subsequently, move toward a SAR-derived, self-consistent bare-earth DEM.

# InSAR Image Classification and DEM Generation

The accuracy of any hydrodynamic model is largely governed by the accuracy of the underlying DEM (Bates, 2004; Schumann et al., 2016). The presence of vegetation, especially in regions along river channels, may obscure or distort bank heights and microtopography, which diminishes height accuracies along the channel and in turn, the overall accuracy of the model (Baugh et al., 2013). These elements play an integral role in hydrodynamic modeling and estimates of bank overtopping. It is therefore imperative that tall features, such as trees and buildings, are successfully classified and removed from the DSM. Information available from the processed InSAR imagery include the measured surface height, the relative power (a function of the backscatter "brightness" of the scene and the radar's illumination geometry), and the interferometric correlation, γ. The interferometric correlation, γ, is a measure of the similarity of the signal received at the two antennas. For a single-pass system, such as GLISTIN, it can be written as the product of a number of factors:

$$\chi = \chi\_{\mathbb{B}} \chi\_{\mathbb{B}} \chi\_{\mathbb{V}}$$

where γ<sup>β</sup> is the geometric decorrelation due to the cross-track separation, B, of the antennas; γ<sup>n</sup> is the decorrelation due to thermal noise and γ<sup>v</sup> is the decorrelation due to the interaction of the electromagnetic wave within the scattering volume. Since the wavelength of GLISTIN is just 8.4 mm, it is dominated by surface scatter for bare terrain; however, γ<sup>v</sup> and thus γ decreases significantly with the presence of trees or vegetation. Buildings, on the other hand, are expected to exhibit a discernible height relative to the ground, while maintaining a high correlation when compared with trees. Indeed, using the same SPInSAR sensor in a snow depth study, Moller et al. (2017) show that the correlation image data alone is robust enough to successfully classify trees.

Following the reasoning above, we used the SPInSAR height image as well as the SPInSAR correlation image for developing the SPInSAR classification algorithm (see flowchart in **Figure 3**). Prior to conducting image classification, noise reduction was performed using an average filter, with a 7 × 7 moving window on the 3 m posted imagery. Note that analysis was conducted for imagery that is ground projected relative to the aircraft flight track (the "sch" coordinate system (Zebker et al., 2010), as opposed a geographic map projection) so that range and incidence-angle dependent behavior is reflected in the classification thresholds for robustness. Vegetated areas and buildings were classified using both height and correlation image data. Mean heights were first calculated for 300 × 300 tiles within the larger swath to both capture local variation and to minimize errors introduced by large groves of trees or buildings. A classification algorithm was initiated for any pixels that exhibited a height greater than 2σ of the mean height for a tile. While the success of this initial criteria is likely owed to the exceedingly low relief of the Red River Basin (RRB), this attribute is not unique to this area, and is often a defining characteristic of flood-prone regions.

Feature classification was performed in a 3 × 3 moving window using the coefficient of variation (CV), or relative standard deviation (RSD), to distinguish between vegetation and built structures:

$$C\_v = \frac{\sigma}{\mu}$$

where σ is the standard deviation of correlation and µ is equal to the mean correlation value for the 3 × 3 window. If C<sup>v</sup> > 0.01 for a 3 × 3 window, pixels are classified as

trees; otherwise, they are broadly classified as buildings. Highresolution, Google Earth imagery was used to identify areas of buildings, forests and patches of trees to evaluate the relative accuracy of the proposed classifier and its parameterization in identifying these features. The classification scheme performed particularly well at identifying tall structures for near and midrange look angles, however, delineation of trees from buildings presented a challenge (**Figure 4**). Here, it is important to note that the main aim of the classification procedure is not to classify tall objects in order to retrieve or count the number of trees or built structures, but rather to create a bare-earth surface, hence classification error may actually be less important. It should be noted though that the threshold value attributed to C<sup>v</sup> for identifying and classifying features should be adjusted for areas of different topographic characteristics and for other feature types for which SAR correlation will be different. For instance, in the case presented here, uniform roofs of large farming structures (abundant in the study region) exhibit high SAR correlation values and thus have a very low coefficient of variation, the value of which was determined using a trial-and-error iterative process. In the present effort, the goal is to derive a bare-earth DEM and so the need to accurately distinguish between different land cover types is certainly less important; however, for other applications this may be the actual objective and thus a sensitivity analysis of SAR correlation statistics for various land cover types may be necessary to conduct in a region of more diverse land cover.

Classified pixels were finally lowered to the minimum height of a 15 × 15 neighborhood. Spatial averaging may be used to further improve the height precision in the InSAR DEM, however, at the cost of resolution (Rodriguez and Martin, 1992; Schumann et al., 2016). As discussed, inundation models are highly sensitive to noise variation in topography, therefore, to improve DEM height precision for applications in flood modeling, the resulting gridded DEM was aggregated to a 30 × 30 m resolution; however, in order to get a better understanding of the influence of topographic complexity on the performance of the proposed classification scheme, a sensitivity analysis ought to be conducted, in which the size of the spatial kernel used to compute local height minima needs to be tied to the degree of topographic complexity (variation), particularly when moving into higher relief terrain. In fact, in Moller et al. (2017) a similar tree classification methodology was found to be very effective in extremely high relief terrain. The sensitivity of the classifier to the terrain relief will be a function of the window size over which height statistics are assessed and this evaluation/sensitivity study has not been done (yet). While there is a great deal of high-relief data already acquired and available, little is over populated areas.

# RESULTS

# Accuracy Assessment

In order to derive a bare-earth DEM, we first removed any residual systematic cross-track tilts using the methodology discussed in see Section "GLISTIN-A InSAR Data" and illustrated by **Figure 3**. At a spatially aggregated DEM resolution of 30 m, this resulted in an overall negligible bias of – 2.5 cm. Convergence of this methodology validates our assumptions with respect to system stability over relevant acquisition time-frames (hours).

Next, an intelligent classifier based on informed relationships between InSAR height and correlation data was used to

spatial aggregation to a bare-earth DEM, with tall vegetation (trees, etc.) and built structures successfully removed. The Google satellite image shown for visual comparison highlights classified groves of trees along the meandering Red River of the North. Note that shorter vegetation and the state of agricultural fields on the satellite image are not directly comparable to the SAR data that were acquired at a different time and during late winter.

distinguish between bare-earth, buildings, and tall vegetation (**Figure 4**). Based on our evaluation of the classified SAR image using high-resolution satellite imagery, the classifier appeared to perform well at identifying tall features, however, it was not successful at distinguishing vegetation from built structures.

We derived a mosaiced, consistent, large-scale bare-earth DEM, which was collected as a series of five flight lines. When assessed at the more typical large-scale DEM resolution of 30 m, we concluded an overall negligible bias of −2.5 cm. Noted that this bias can be readily improved if required by increasing the signal-to-noise ratio of the data. Instrument customization, operational customization, and trade-offs enable this analogous-to-acquisition of LiDAR design, where point density and resolution are traded against acquisition cost. Using LiDAR terrain data and GPS ground control points from the National Geodetic Survey (NGS, **Figure 5**), we demonstrate suitability to a wide range of applications, particularly flood mapping, by achieving a RMSE of approximately 0.5 m (**Table 2**).

The height errors between SPInSAR and LiDAR reported in **Table 2** are all computed at the location of the GPS network points, so as to avoid using LiDAR as a "truth" reference but instead as a direct comparison to the SPInSAR technology. It is, however, worth noting that in terms of overall height performance across the study region, LiDAR may be expected to perform better on average given the inherent noise level in the raw SAR data; although performance discrepancies between the two technologies may be acceptable depending, of course, on the application.

# DISCUSSION

# LiDAR vs. GLISTIN-A

For obvious reasons, obtaining a high-resolution, high-accuracy LiDAR terrain (i.e., bare earth) model over floodplains and lowlying coastal areas is preferable – in particular, where such areas have very shallow gradients, such as the RRN floodplain. With a high signal pulse density per unit area and almost no random noise component, airborne LiDAR technology is capable of resolving minimal topographic variations and microtopographic features that define and control non-trivial floodplain flow pathways. In fact, LiDAR, with a typical vertical RMSE of 15– 20 cm, helped shift the flood modeling community from a traditionally data-poor environment to a data-rich environment, which led to a large number of model advances and a rethinking of model design (Bates, 2004). However, as already argued by Schumann et al. (2016), over large, wide floodplains, LiDAR can quickly become prohibitively expensive due to its typically small swath width and thus, requires many overlapping flight lines to cover the total area of interest at high resolution. For large area coverage, wide-swath, single-pass airborne InSAR, therefore, presents notable advantages and has the potential to become much more cost-effective.

In this context, an airborne campaign with overpasses at different times (e.g., at different seasons) using a GLISTIN-A type sensor can be easily envisaged and becomes indeed very competitive in resource allocation and cost. A GLISTIN-A flight pass was repeated over the same study site in snow

and ice conditions. This data set, coupled with the dry baseline DTM, allows identification of river ice jams as well as volumetric estimation of the floodplain snowpack. While examining the latter capability for this site is outside the scope of this study, snowpack volume estimation with GLISTIN-A has been successfully demonstrated over the Sierra Nevada (CA, United States) by Moller et al. (2017).

#### Flood Mapping Potential

The above-mentioned attributes make GLISTIN-A type technology highly relevant and applicable to flood mapping and modeling. It is well known that while variations of in-channel water levels (determined by local flow conditions) drive the

TABLE 2 | Calculated errors for the SAR-derived DEM and LiDAR relative to GPS ground-control points, as well as SAR vs. LiDAR.


timing and amount of river bank overtopping and subsequent flooding of adjacent low-lying land, it is variations in floodplain topography that control floodplain flow paths and the inundated area during a flood event; thus, microtopography and floodplain features, such as buildings, walls, trees, etc., become significant, particularly when interested in localized flow conditions and associated floodplain inundation at the small scale (Mason et al., 2003). Microtopographic features and variations in microtopography are only included in flood inundation (i.e., 2-D hydrodynamic) models when high-resolution, highprecision data on floodplain heights are available. In most cases, however, their effects are parameterized in models of grid resolutions that are typically orders of magnitude larger than the microtopographic controls (Dottori et al., 2013).

Within this context, Mason et al. (2003) state that the typical vertical accuracy of floodplain heights from airborne LiDAR (10 to 20 cm RMSE) provides a realistic lower limit for DEM quality as beyond this, the sensor signal becomes indistinguishable from "noise" generated by background microtopographic features. In line with this reasoning, it is fair to assume that a vertical error of 2σ of the lower limit (∼40 cm) sets an adequate upper limit target value.

In the case study presented here, the DTM generated from the GLISTIN-A InSAR data had a mean absolute error of 0.421 m and a RMSE of 0.537 compared to geodetic ground control points. This makes it suitable for flood mapping and modeling, with much better accuracy than can be obtained with conventional wide-swath, imagery-derived DTMs, such as from satellite InSAR technology or commercial multi-pair optical imagery, which typically exhibit much larger vertical errors.

#### CONCLUSION

In this study, we demonstrate the success of a new InSAR calibration technique that removes the need for ground control points and LiDAR by leveraging areas of swath overlap in the same airborne acquisition. Furthermore, we showcase the potential of InSAR image classifiers to generate floodplain DEMs with vertical accuracies suitable for high-resolution hydrodynamic modeling. Averaging to a spatial resolution of 30 m after classification, a mean error in the vertical as low as 12 cm was obtained. More research on classification algorithms is indeed required to automate the DEM-generation process and to make such classifiers readily available in commercial software. The results of this case study, however, demonstrate the potential of airborne InSAR as a cost-effective alternative to LiDAR for large-scale flood mapping, thus diminishing barriers

#### REFERENCES


to acquisition in the developing world and furthering the prospect of a free, global, high-resolution DEM. The InSAR DEM presented here will be examined in a follow-on study in the RRN floodplain, where its performance will be evaluated using a real-case, 2D flood model application.

#### AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

### FUNDING

DF's time on this work was supported by the COAST internship program. The GLISTIN-A flights were funded under NASA's Terrestrial Hydrology Program. The GLISTIN-A data are available from the JPL UAVSAR data website (https:// uavsar.jpl.nasa.gov/cgi-bin/data.pl) and the DEM can be shared upon request.

#### ACKNOWLEDGMENTS

The authors thank the NASA-JPL UAVSAR team for the GLISTIN-A data acquisition and calibration.


**Conflict of Interest:** GS and DM are employed by Remote Sensing Solutions.

The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Faherty, Schumann and Moller. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Multiscale Integration of High-Resolution Spaceborne and Drone-Based Imagery for a High-Accuracy Digital Elevation Model Over Tristan da Cunha

#### Dietmar J. Backes\* and Felix Norman Teferle

Geodesy and Geospatial Engineering, Institute of Civil and Environmental Engineering, University of Luxembourg, Luxembourg City, Luxembourg

#### Edited by:

Paul Bates, University of Bristol, United Kingdom

#### Reviewed by:

Ahmed Kenawy, Mansoura University, Egypt Brian C. Gunter, Georgia Institute of Technology, United States

> \*Correspondence: Dietmar J. Backes DJBackes@outlook.com

#### Specialty section:

This article was submitted to Hydrosphere, a section of the journal Frontiers in Earth Science

Received: 10 October 2018 Accepted: 07 July 2020 Published: 08 September 2020

#### Citation:

Backes DJ and Teferle FN (2020) Multiscale Integration of High-Resolution Spaceborne and Drone-Based Imagery for a High-Accuracy Digital Elevation Model Over Tristan da Cunha. Front. Earth Sci. 8:319. doi: 10.3389/feart.2020.00319 Very high-resolution (VHR) optical Earth observation (EO) satellites as well as lowaltitude and easy-to-use unmanned aerial systems (UAS/drones) provide ever-improving data sources for the generation of detailed 3-dimensional (3D) data using digital photogrammetric methods with dense image matching. Today both data sources represent cost-effective alternatives to dedicated airborne sensors, especially for remote regions. The latest generation of EO satellites can collect VHR imagery up to 0.30 m ground sample distance (GSD) of even the most remote location from different viewing angles many times per year. Consequently, well-chosen scenes from growing image archives enable the generation of high-resolution digital elevation models (DEMs). Furthermore, low-cost and easy to use drones can be quickly deployed in remote regions to capture blocks of images of local areas. Dense point clouds derived from these methods provide an invaluable data source to fill the gap between globally available low-resolution DEMs and highly accurate terrestrial surveys. Here we investigate the use of archived VHR satellite imagery with approx. 0.5 m GSD as well as low-altitude drone-based imagery with average GSD of better than 0.03 m to generate high-quality DEMs using photogrammetric tools over Tristan da Cunha, a remote island in the South Atlantic Ocean which lies beyond the reach of current commercial manned airborne mapping platforms. This study explores the potentials and limitations to combine this heterogeneous data sources to generate improved DEMs in terms of accuracy and resolution. A cross-validation between low-altitude airborne and spaceborne data sets describes the fit between both optical data sets. No co-registration error, scale difference or distortions were detected, and a quantitative cloud-to-cloud comparison showed an average distance of 0.26 m between both point clouds. Both point clouds were merged applying a conventional georeferenced approach. The merged DEM preserves the rich detail from the drone-based survey and provides an accurate 3D representation of the entire study area. It provides

the most detailed model of the island to date, suitable to support practical and scientific applications. This study demonstrates that combination archived VHR satellite and low-altitude drone-based imagery provide inexpensive alternatives to generate high-quality DEMs.

Keywords: digital elevation model, drone photogrammetry, VHR satellite imagery, Tristan da Cunha, multiscale data fusion

#### INTRODUCTION

Today digital elevation data at medium spatial resolution, e.g., 30 m GSD, are freely available at almost global scale from the Shuttle Radar Topography Mission (SRTM) and the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Map (GDEM). Conventionally, high-resolution topographic data (bellow 0.5 m GSD) are captured by dedicated mapping cameras and LiDAR sensors deployed on manned aircraft. Although these technologies became very effective during the last decades, they are still costly and can be difficult to deploy in remote and hard-to-reach regions. High-resolution DEMs derived from VHR satellite images and drone-based data can provide a valuable alternative where airborne data are not available.

The definition of the terms low-, medium- or high-resolution varies between scientific disciplines with respect to the deployed platform and sensors (e.g., spaceborne, high-, medium- or lowaltitude airborne) as well as derived data products (images, DEMs and point clouds). Furthermore, technical progress in sensor and computer technologies seem to shift the definitions toward smaller GSDs over time. In the scope of this study, we use the classification selected by Dowman et al. (2012) for satellite imagery and expand this definition to DEMs. The geometric resolution of low-resolution data is larger than 30 m GSD, medium resolution data ranges from 30 – 5 m GSD, high-resolution data ranges from 5 – 1 m GSD and very highresolution data shows a spatial resolution of less than 1 m. Since the spatial resolution of low-altitude drone imagery, e.g., 0.03 m, is at least one magnitude higher than VHR satellite images (0.30 m) these data are referred to as drone-based data.

Earth Observation missions such as Landsat and Sentinel provide a free and frequently available source of remote sensing imagery at medium resolutions (Dowman et al., 2012), which enable numerous applications such as land cover mapping, environmental monitoring and change detection (Aschbacher and Milagro-Pérez, 2012; Ban et al., 2015). Most applications also require digital elevation data, which are usually provided by SRTM and ASTER GDEM. However, many studies have reported shortcomings regarding completeness, artifacts, elevation errors and co-registration problems in these data sets (Shortridge, 2006; Tighe and Chamberlain, 2009; Guth, 2010; Shortridge and Messina, 2011; Rexer and Hirt, 2014). To overcome this problem Robinson et al. (2014) proposed a quasi-global, voidfree, multiscale smooth DEM fused from ASTER and SRTM data at 90 m GSD. Higher quality and resolution DEMs at 5 m GSD or better are only commercially available at relatively high costs, e.g., Airbus WorldDEM (based on Tandem-X) or Vricon DSM (Schumann et al., 2016).

The extraction of DEMs from satellite imagery has been studied and demonstrated for decades. With the launch of the Satellite pour l'Observation de la Terre (SPOT), it was possible to generate 3D digital elevation models from satellite imagery for the first time. Using a cross-track stereo configuration, satellite images were acquired from neighboring orbits form stereo-pairs that are covering a defined region of interest. Early research showed encouraging results and verified the suitability for topographic mapping applications (Gugan and Dowman, 1988; Dowman, 1989; Muller and Dowman, 1989; Theodossiou and Dowman, 1990). However, the cross-track stereo configuration bears many challenges regarding image geometry and temporal decorrelation that influence the quality of the extracted DEM. Despite algorithmic improvements of image matching and terrain extractions, as well as improved satellite coverage, the challenges identified at the time still apply. More recent missions deployed dedicated optical sensors which allowed the acquisition of in-track stereo images. SPOT5 Highresolution Stereo (HRS) and the Advanced Land Observing Satellite (ALOS) Panchromatic Remote-sensing Instrument for Stereo Mapping (PRISM) showed improved performance (Fraser, 1997; Michalis and Dowman, 2004; Rudowski, 2004). Today, fleets of highly agile optical remote sensing satellites can collect multispectral VHR imagery with up to 0.3 m GSD. Furthermore, rapid technological advances in orbit modeling and satellite altitude control enable a very accurate determination of the direct georeference. For example, WorldView scenes have an absolute positioning accuracy of less than 5 m circular error (Digital Globe, 2016). Many of these satellites are also able to collect intrack stereo-pairs that allow high-accuracy DEM extraction. The generation of DEMs from VHR spaceborne stereo images over urban areas has been investigated by Gong and Fritsch (2016).

With the launch of the next generation of VHR optical satellites, it has also been recognized that the differences in characteristics between airborne and spaceborne images are decreasing (Dowman et al., 2012 p. 107; Chen et al., 2016). Thus, the quality of DEMs derived from satellite images ought to increase toward the quality of DEMs derived from airborne mapping sensors. It is expected that fleets of increasingly efficient satellites will collect vast amounts of image data, which will enable the extraction of detailed DEMs in almost any part of the world with increasing temporal frequency.

In this study, we report on DEMs derived by digital photogrammetric methods using three comprehensive geospatial software suites capable of processing VHR imagery. They all support common spaceborne sensor models and data formats, providing preprocessing capabilities, triangulation modules as well as surface reconstruction using dense image matching (DIM) algorithms to extract high-resolution digital surface models

(DSMs). In "dense" reconstruction, one elevation is generated per pixel by matching epipolar constraint images. Energy based methods, which aim to find the best global solution (Pierrot-Deseilligny and Paparoditis, 2006; Hirschmuller, 2008; Strecha et al., 2010; Bulatov et al., 2011) produce the best results to date whereby semi-global matching (SGM) by Hirschmuller (2008) has been widely adopted in recent years. Digital surface models extracted by these methods usually include fine details, but they are also more prone to noise. Multiple overlapping satellite images can be treated as so-called multi-view stereo (MVS) configuration, which can increase the reliability of the extracted elevation models.

At the same time, low-cost, consumer-grade drones or unmanned aerial systems (UAS) have become a powerful and flexible tool to collect large blocks of images suitable to generate highly detailed point clouds, mesh models, DEMs and orthophotos with a geometric resolution of a few centimeters. Such systems are only able to cover small areas due to their limited flying time and usually deploy low-precision navigational sensors and consumer-grade cameras (Colomina and Molina, 2014). The limited endurance and data quality requirements need to be considered carefully during data capture and analysis. Many studies to investigate best practices and the accuracy potential of drone surveys have been conducted (Haala et al., 2012; Day et al., 2016; Gerke and Przybilla, 2016; Cramer et al., 2017; Gindraux et al., 2017). However, as drone technology continues to improve, e.g., the introduction of real-time kinematic (RTK) systems, the capabilities and accuracy potential of drone surveys is still evolving. Post-processing is often done by so-called structure from motion (SfM) engines, which are usually employed as black-box systems. They allow minimal user input and are only providing final results as a point cloud or mesh model without quality control. Drone photogrammetry packages, such as Pix4D (Mapper or Engine), do provide more flexibility and detailed reports with photogrammetric quality metrics. Similar to highly efficient SfM engines and simultaneous location and mapping (SLAM) algorithms, highly automated digital photogrammetric methods have advanced rapidly during recent years. However, "the names of SLAM, SfM and photogrammetry do mean different things, notwithstanding the presence of conceptual and algorithmic overlaps. It should be acknowledged unequivocally that neither SLAM nor SfM encompasses the full process of metric-quality image-based 3D measurement, object reconstruction and mapping that is today's highly automated photogrammetry" (Fraser, 2018).

Over the last decade, numerous scientific studies have exploited free medium data sources. However, many applications, e.g., landslide and flood modeling, would benefit from even higher quality DEMs that provide better resolution with a higher level of detail and improved positional accuracy. This study explores the potential to combine existing VHR satellite and low low-altitude drone imagery to generate high-resolution DEMs on a practical example, the remote island of Tristan da Cunha. After the introduction of the study site and the available data sets, an applied workflow to generate and merge DEMs is described in the method section. The results are introduced and analyzed in detail. A discussion revisits best practice procedures, tools and methodologies, reflects on challanges and its applicability to similar study sites. The conclusion summarizes the results and shows future perspectives.

# STUDY SITE AND DATA SET

Amongst a number of important research installations, Tristan da Cunha hosts a recently established continuous Global Navigation Satellite System (GNSS) station and tide gauges in support of global reference frame and sea-level studies. The installation of these sensors provided an opportunity to conduct a local survey to map the local area around the harbor of Edinburgh of the Seven Seas and enable this study.

The acquired data set consists of freely available topographic mapping data (approx. 1/10,000), a large stack of archived VHR satellite imagery from Maxar/DigitalGlobe, a low-altitude drone survey as well as geodetic ground truth data (**Table 1**). While the topographic mapping data and the satellite images provide full coverage of the island, the drone-based images and geodetic ground truth only cover a limited local area around the harbor (**Figure 1**).

# Tristan da Cunha

The small volcanic island of Tristan da Cunha is known as one of the most remote locations on the Earth. It is located at 37◦ 040 S and 12◦ 19<sup>0</sup> W in the South Atlantic Ocean approximately 1950 km west of Cape Town, South Africa, and 2900 km east of Buenos Aires, Argentina. Tristan da Cunha has a circular shape with a diameter of 10-12 km, covers a surface of 96 km<sup>2</sup> and its highest point is Queen Mary's Peak with 2062 m (**Figure 1A**). The last eruption of the volcano dates to 1961. The island is covered by low vegetation including grass and shrubs. Fertile but narrow coastal plains on the north-western and southern regions are followed by steep slopes that rise to approximately 500 – 800 m above the sea level. At this height the topographic gradient flattens to form a high plateau before rising to the highest peak of the island, which has the classic conical shape of an active volcano. Tristan da Cunha is also home to a population of approximately 280 people, living in the settlement of Edinburgh of the Seven Seas at the north-western coastal plain **(Figure 1B**). The settlement is served by a small harbor which lies exposed to strong winds and rough seas. Woodworth (2020) provides a detailed overview on the development of the harbor and the local climatic conditions.

Due to its exposed location the island has a cold temperature oceanic climate with gusty winds and gale forces storms. A persistent cloud coverage provides regular precipitation in excess to 1500 mm at sea level and over 3000 m on the mountain. An increase in coastal erosion and the frequency of landslides have been reported by the local population, which raises concerns for sustaining its inhabitants in the long-term future. Since only medium-resolution elevation data exist, there is a need for accurate 3D mapping. However, its remote location and rough weather conditions provide exceptional challenges for terrestrial, aerial as well as spaceborne data acquisitions for mapping and environmental monitoring.

FIGURE 1 | (A) Satellite image of Tristan da Cunha courtesy of the DigitalGlobe Foundation. (B) Zoom to Edinburgh of the Seven Seas. (C) Perspective view based on drone imagery. (D) Low-altitude drone image.

### Global Low and Medium-Resolution Mapping Data Open Street Map

A topographic map of the island was available via Open Street Map (OSM). Unfortunately, little metadata about currency, source and quality were obtainable. To evaluate the accuracy of the geolocation, OSM was compared to 6 control points derived from the GCN described in section "Low-Altitude Aerial Images." Since these control points were only located around the small local harbor area, they only provide reference information to estimate the translations in Easting and Northing but no reliable rotation of the map. The OSM map is displaced by approximately 37 m in Easting and 69 m in Northing, as shown in **Figure 5**.

#### SRTM and Aster-GDEM

Shuttle Radar Topography Mission and ASTER-GDEM provide free medium resolution DEM at the Digital Terrain Elevation Data standard (DTED) level 2 standard with 1 arc-second or approximately 30 m. Both data sets were acquired via United States Geological Survey Earth Explorer (USGS, 2018). The SRTM 1 arc-second global DEM was declassified for public use in 2014 and offers almost worldwide coverage with a substantially increased level of detail (**Figures 2A,B**). Despite the promise of complete coverage with no voids, both products still show some data voids on the southern part of the island and close to the summit of Queen Mary's Peak. The ASTER-GDEM was jointly developed by Japan's Ministry of Economy, Trade, and Industry (METI) and the United States National Aeronautics and Space Administration (NASA) (Tachikawa et al., 2011). Its current version, which is available since 2011, benefits from many improvements compared to the first version, but still contains artifacts and elevation errors on a local scale (Tachikawa et al., 2011). The DEM over Tristan da Cunha shows a slightly better level of detail at 30 m GSD compared to SRTM. It does not include data voids but shows some noise (**Figure 2C**).

To increase reliability and accuracy, both models were fused together using BAE Systems' SocetGXP terrain-merge algorithm (**Figure 2D**). To increase reliability and accuracy, both models were merged. Following a statistical assessment of the co-registration between both models the Aster-GDEM was merged with the SRTM DEM which served as Reference DEM. After a fine registration with the removal of bias and outliers, both models are merged by combining and interpolating height values of corresponding raster cells. This was done using the BAE Systems SocetGXP terrain merge function.

FIGURE 2 | (A) SRTM DEM with 3<sup>00</sup> GSD. (B) SRTM DEM with 1<sup>00</sup> GSD. (C) ASTER-GDEM Version 2.0 with 1<sup>00</sup> GSD. (D) mDEM with 1<sup>00</sup> GSD.

The merged model does not include data voids and shows a reduced level of noise compared to the Aster-GDEM but an increased level detail compared to the SRTM DEM. This medium resolution DEM (mDEM) at 30 m GSD was subsequently used as so-called 'seed DEM' in the stereo matching process using the VHR imagery.

#### VHR Optical Earth Observation Data

A comprehensive stack of VHR images was provided by DigitalGlobe which included multispectral images from QuickBird (QB), WorldView2 (WV2) and WorldView3 (WV3). Overall 64 basic and 32 Standard2A image products were delivered. While Standard2A image products are already georectified using an unspecified coarse DEM, basic image products are only radiometrically corrected thus suitable for terrain extraction. A detailed overview of product specifications is given in DigitalGobe's technical manual (Digital Globe, 2014, 2016). **Table 2** provides a detailed overview of the basic satellite scenes used in this study. The scenes have spatial resolutions of up to 0.61 m for QuickBird and 0.31 m for WorldView Satellites and were acquired between February 2009 (QB) and July 2017 (WV3). Although many acquisitions include cloud coverage and atmospheric effects, a number of images were identified which provide excellent image quality over the complete extent of the island. The images have been acquired from a variety of viewing angles and are suitable for the extraction of DEMs. **Figure 3** shows the "WV2 image stack" and gives an overview of the image quality.

#### TABLE 1 | Available data sets and sources with spatial resolution, coverage and accuracy.


\* Technical Reference Digital Globe, 2016.

TABLE 2 | DigitalGlobe basic image products, scenes of the prime image triplet are marked in green, scene of the MVS image stack marked in gray and green.


FIGURE 3 | (A) WorldView 2 image stack. (B) Geometric resolution and feature definition. (C) Lack of co-registration between QB and OSM. (D) WorldView3 Scene showing clouds and haze. (E) WorldView3 scene showing shadow and haze.

#### Low-Altitude Aerial Images

An off-the-shelf DJI Phantom 3 Professional drone was used to acquire aerial photography around the harbor area and research installations. The mission aimed to collect drone imagery suitable to create detailed maps of a limited area covering the new GNSS and existing Doppler Orbitography Radiopositioning Integrated by Satellite (DORIS) stations, the new tide gauges in the harbor but also to collect panoramic photographs of the surrounding areas for documentation purposes.

Weather conditions permitted only 2 flight missions on October 8 and 18, 2017. Due to the harsh conditions with high wind speeds as well as the identified offsets in the base mapping data, both missions were executed as free flight missions in manual mode. A dense block of images was captured from a range of different flying heights covering the area of interest from different perspectives and viewing angles. Overall, 373 usable images were acquired which provide good overlapping coverage of the areas of interests as well as scenic panorama shots showing the region beyond the settlement of Edinburgh of the Seven Seas. Although not intended for 3D modeling, these images were subsequently added to the analysis. **Figure 4** gives an overview of the acquired block of images. Ground control was provided by ten well-distributed targets (**Figure 5**) which were established by the geodetic survey.

# Geodetic Ground Control

A highly accurate ground control network (GCN) with subcentimeter absolute positional accuracy, which was required for the sensor installations, was established around the area of the harbor and research installations. The geodetic survey, based on a combination of GNSS and total station observations, also provided Ground Control Points (G) to support the drone survey with sub-centimeter positional accuracy.

The continuous GNSS station TCTA (1035) and point 1003 were observed successively for 3 days. The coordinates were derived from the daily observation files using precise point positioning (PPP) and employing the Bernese GNSS Software V5.2 (Dach et al., 2015). The coordinates were then averaged over the 3 days. Subsequently Stations 1002, 1006 and 1005 were occupied for 2 h each, and their coordinates were obtained from static processing using Leica GeoOffice 8.0 in which the coordinates of TCTA were fixed to the values obtained from PPP. In all processings, precise IGS satellite orbit and clock products were used as well as individual absolute GNSS antenna calibrations.

A terrestrial survey was conducted using a high-precision total station (LeicaTS30) and a high-precision digital level (Leica DNA03). The primary aim of this survey was to establish the local GCN, which will enable the monitoring of the GNSS station TCTA, the tide gauges and the DORIS station. This highly precise GCN provided the basis to establish the control used for the aerial survey. Using the Leica TS30 total station, the coordinates of ten GCPs have been determined to an absolute accuracy of better than 1 cm in XYZ. **Figure 5** depicts the location of the GCN as well as the GCPs.

FIGURE 4 | (A) Flight missions overlaid on OSM in 2D. (B) Perspective view shows different altitudes and image perspectives. (C) Image positions after bundle block adjustment. (D) Orthophoto with overlaid image positions.

# METHODOLOGY

This study investigates the combination of archived VHR satellite and low-altitude drone imagery to generate highresolution DEMs. The main perspective was to implement and verify rigorous photogrammetric workflows based on existing commercial off-the-shelf tools. Point clouds from both data sources are merged to provide a topographic representation of the entire study site while preserving the rich details captured by the drone survey over the central area of interest.

Throughout the study, all data sets were rigorously projected into Universal Transversal Mercator (UTM) Zone 28H (South) coordinates based on the World Geodetic System (WGS) 1984 datum. With the absence of local gravimetric data, the vertical datum was defined as WGS84 ellipsoidal heights.

# Digital Terrain Extraction From Optical Satellite Imagery

Exploiting the VHR imagery, various DEMs, which cover the entire island, were extracted following photogrammetric methodologies with the intent to maximize spatial resolution, completeness and accuracy. In the scope of this study, we tested the performance of three geospatial software packages. Since the implementations of the photogrammetric workflows inside these packages differ in detail, all required processing steps are described and mapped to a general workflow shown in **Figure 6**.

The initial data preparation includes the selection and preprocessing of the satellite scene. This is followed by a spaceborne triangulation using tie points as well as ground control. The results of the triangulation provide accurately georeferenced and well-oriented satellite scenes. The exact alignment between the scenes is crucial for the following terrain extraction steps. All software packages provide several image matching algorithms to generate DSMs as rasterized models, irregular triangular networks (TIN) or point clouds. A final verification and filtering process is required as the results of the matching process always contains artifacts and noise.

Most terrain extraction algorithms support the use of seed DEMs to improve computational efficiency and reduce the number of blunders. Usually, freely available low or medium resolution DEMs (e.g., DTED level 0) are used for this purpose. As higher quality seed DEMs increase performance, the mDEM (2.1) was used as initial seed DEM whenever possible. The general

FIGURE 5 | Geodetic control network (red triangles) with photogrammetric GCP (yellow circles) overlaid on OSM; the figure also shows the displacement of OSM compared to the GCN.

terrain extraction workflow followed a stepwise procedure. Using the mDEM with 30 m GSD as seed, a DEM with 10 m GSD was matched. In the next iteration this DEM was injected as a seed to extract a DEM with higher spatial resolution. Following this gradual approach, DEMs with 10, 5, 2 and 1 m GSD were extracted.

#### Image Verification and Preprocessing

The DigitalGlobe VHR image repository provided a rich data source to derive a high-resolution DEM of the island. However, temporal decorrelation (e.g., changes in landcover between image acquisitions), atmospheric effects that influence image quality (e.g., haze), cloud coverage, shadows, unfavorable viewing angles, differences in spatial resolution as well as partial coverage, are presenting challenges for the terrain extraction process. Thus, the choice of well-suited scenes and the adequate preprocessing have a direct impact on the quality of the derived DEMs.

The initial data verification step aimed to select the most suitable satellite scenes from the stack of 64 basic image products. A triplet of satellite scenes was identified which showed excellent image quality, covered the entire island almost cloud-free and provided good geometric conditions for DEM extraction. Using a MVS configuration based on a large stack of overlapping images is expected to increase the reliability of the extracted DEMs, but is more computationally expensive and relies on good and homogenous image quality. Additional 10 images were chosen over the area of the settlement, to investigate the difference in performance and quality between the image triplet and a MVS images stack.

The selection of satellite scenes from the repository followed the subsequent criteria in the given order:


Since the selected satellite scenes showed excellent image quality, extensive preprocessing which includes atmospheric correction, cloud masking as well as pan-sharpening was not required. In this study, only the panchromatic channel was used for DEM extraction.

#### BAE Systems – SocetGXP

SocetGXP includes a powerful photogrammetric toolset, which follows rigorous procedures and allows a high degree of flexibility and control. The software package emerged from SocetSET, which was initially developed by Helava as one of the first digital photogrammetric workstations. Today, the comprehensive software suite is mainly used by the geospatial intelligence community.

Since DigitalGlobe formats and sensor models are supported in SocetGXP the import of large scenes via the TIL file format was straightforward and contained the sensor orientation via rational polynomial coefficients (RPC). The subsequent triangulation was performed using manual well distributed and automatically measured tie-points. Since the GCN only covers a small area around the harbor area, only one GCP was derived as a single point of reference which only allows to detect and correct for a shift in position. The results of the bundle block adjustment were used to update the RPC information of the satellite scene.

SocetGXP includes 3 algorithms for DEM extraction optimized for various types of data and terrain. This study utilized the classical auto terrain extraction (ATE) module based on normalized cross-correlation (NCC) as well as the so-called next generation terrain extraction (NGATE) module, which provides an advanced algorithm that "combines area-matching and edge-matching, using each approach to guide the other and thus reduces blunders" (Walker, 2007). The theory and algorithms of NGATE have been described by Zhang (2006) and Zhang et al. (2006). In both methods the merged mDEM was utilized as a seed model.

#### Hexagon Geospatial–Erdas Imagine

The widely used Erdas Imagine package is a comprehensive geospatial software suite, which integrates GIS, remote sensing and photogrammetric functionalities. Similar to SocetGXP, DigitalGlobe's sensor models and data formats are well supported.

The satellite scenes were triangulated as an image block. Since automatic tie point detection did not provide usable results, only well-distributed manually measured tie points where used in the triangulation. Similar to the analysis in SocetGXP, one GCP derived from the GCN provided control in 3 translations (X, Y, Z). Since the quality of the image alignment impacts the performance of stereo matching algorithms, the results of the triangulation were used to modify and improve the satellite orientation parameters (RPCs).

Hexagon Geospatial offers a range of terrain extraction modules optimized for different sensors or applications. This study examined the performance of the automatic terrain extraction module (ATE), the enhanced terrain extraction (eATE) and the SGM implementation by Tridicon. While the ATE module is based on an undisclosed feature-based image-matching algorithm, eATE provides a range of options to optimize matching and processing performance, reduce blunder, filter noise and classify the results. Here we used the NCC matching algorithm in combination with standard blunder detection based principal component analysis (PCA) and reverse matching. Filtering and classification were not applied and the results were provided as raster DEM and point cloud. The SGM module was initialized based on the results of the triangulation. Single stereopairs are matched to partial point clouds, which are merged subsequently to a combined point cloud.

#### PCI Geomatics – Geomatica OrthoEngine

PCI Geomatics provides a very efficient remote sensing software tool kit, which offers flexibility through an extensive library of algorithms accessible via a Python API. Using the Geomatica OrthoEngine module, digital elevation models from single stereopair were derived via different image matching algorithms, including NCC and SGM. Only one stereo-pair was selected, which was acquired on the same orbit and showed an almost perfect alignment. Epipolar images were generated as input for stereo matching using the NCC algorithm to produce a highquality DEM.

# Drone Photogrammetry

feart-08-00319 September 8, 2020 Time: 10:47 # 11

Generating 3D point clouds from image blocks using highly integrated SfM engines, is a mostly automatic process with few options for user inputs and control. In the scope of this study, we report on the results created with the drone photogrammetry packages Pix4D Mapper (4.2.25). Pix4D provides detailed data quality reports, which contain the results of the aerial triangulation, camera calibration, quality estimations of the 3D positioning as well as metadata of the derived 3D products.

To investigate the performance of the low-precession GNSS/INS sensors onboard the drone, an aerial triangulation was computed without any ground control points. Only the positions measured by the onboard GNSS receiver where used. In a first step only images of the "free flight image block" were included in the triangulation. To strengthen the geometry of the aerial image block further and to provide better coverage of occluded areas, the number of images was increased by adding the scenic images originally intended for documentation purposes in a second step. The final aerial triangulation included all 373 usable images. Subsequently the GCPs derived from the geodetic survey where introduced stepwise as checkpoints (CP) and ground control points.

Investigations on the stability of camera systems deployed on small drone systems by Cramer et al. (2017) and Gerke and Przybilla (2016) found that such lightweight cameras are often geometrically unstable and require frequent recalibration. Thus all triangulations included a self-calibration for the interior orientation of the drone camera.

# Data Quality Control and Multiscale Data Integration

The accuracy and quality of the image orientation are given by the results of spaceborne or aerial triangulation. This includes the alignment of the images relative to each other as well as the positional accuracy. These parameters are also indicators for the relative and absolute accuracy of the derived DEM. As already mentioned ground truth was only available around the harbor area via the geodetic surveyed GCN which showed subcentimeter absolute positional accuracy. For the drone survey nine well distributed photogrammetric targets provided adequate ground control. However, due to the small extent of the GCN in comparison to the overall island, the GCN could only provide a single point of reference for the spaceborne image block. One GCP and one additional CP of similar definition and positional accuracy were derived from the GCN to estimate and improve the positional accuracy of the directly georeferenced satellite scenes.

The generated DEMs were inspected for blunders, noise, completeness and, subsequently, cleaned and filtered. This was followed by cross-validation between derived spaceborne DEMs and the dense point clouds derived from the drone survey. Due to a lack of suitable reference data a thorough assessment as described by Höhle and Höhle (2009) or Höhle and Potuckova (2011) of the DEMs could not be conducted for this case study but would be highly desirable.

Finally, both models from spaceborne VHR photogrammetry (**Figure 7A**) as well as drone photogrammetry (**Figure 7B**) were merged based on the exact georeference provided by the GCN. The point clouds provided a flexible data format to merge such heterogeneous models, which differ in extent and scale. Two dense point cloud models derived by SGM at approximately 1 m and 0.5 m spacing were chosen as good representations for the entire island and merged with the rich point cloud from the drone photogrammetry. The latter covered the local area as well as the extended areas near the settlement (**Figure 7B**) for comparison. The final point cloud preserves fine local detail and provides a good overall representation of the area.

# RESULTS AND ANALYSIS

This section reports on the achieved results of terrain extraction from satellite imagery and the drone photogrammetry. It provides a cross-validation and shows the fused model as a combined point cloud from spaceborne and low-altitude drone photogrammetry.

# High-Resolution DEMs From Spaceborne Images

Following the described photogrammetric workflows highresolution DEMs were extracted from VHR imagery using the three popular software packages. All results show a high level of completeness and a low noise level after blunder removal. The absolute positioning accuracy was verified using control points derived from the GCN.

#### VHR Image Verification, Selection and Preprocessing

Through the initial selection process, three satellite scenes (image triplet) were identified from the stack of archive data which provided excellent image quality, suitable viewing angles for stereo matching, low temporal decorrelation and full coverage of the area of interest. The selected scenes also showed a good pre-alignment from direct georeference. An image pair captured in November 2015 during the same orbit provided a cloud free view of the entire island. The B/H ratio was calculated to be approximately 0.6 which indicates ideal conditions for stereo matching. The third image was captured in March 2012 from a parallel orbit with a different viewing angle.

#### Spaceborne Triangulation

The Image triangulation was performed in SocetGXP and Erdas Imagine. Since the selected image triplet showed good pre-alignment via direct georeference an initial triangulation was performed using only tie points. To avoid unintended error extrapolation and wrapping effects to the final model, only one GCP and one CP from the GCN were included in the triangulation.

FIGURE 7 | (A) Spaceborne DEM based on WV2 imagery with 2 m resolution generated in SocetGXP. (B) Point cloud from drone photogrammetry using Pix4D with a priori camera positions (blue spheres) and after bundle block adjustment using geodetic GCPs (green spheres).

#### **Triangulation in SocetGXP**

Two separate triangulations were performed. The first triangulation included the image triplet over the full extent of the island. Overall 24 well-distributed tie points, as well as the GCP and CP, were chosen and measured in a semi-automatic mode. The result of the triangulation shows RMSE of 0.48 pixels for the image residuals and a shift of 0.17 m in X, 0.05 m in Y and 0.04 m in Z. Considering a GSD of approx. 0.5 m, which limits the precision to measure ground control points in images, this result shows agreement between satellite direct orientation data and the geodetic control points. The second triangulation was based on the stack of 10 WV2 images, which focused on the local area of Edinburgh of the Seven Seas. This time 60 tie points and the same control points as in the first triangulation were measured in order to provide sufficient tie points in the image stack. The result showed a similar RMSE of 0.63 pixels for the image residuals but a higher shift of 1.02 m in X, 1.56 m in Y and 0.31 m in Z.

#### **Triangulation in Erdas imagine**

For the triangulation in Erdas Imagine, 13 manually selected and well-distributed tie points were measured. The tie points were located in similar regions and the GCP was identical as for the triangulation in SocetGXP. The triangulation shows a RMSE of 0.08 pixels for the image residuals and a shift of 0.04 m in X, 0.08 m in Y and 0.12 m in Z. Therefore, the results are similar to those from SocetGXP. The CP provides a validation of the triangulation results, and showed a shift of 0.41 m in X, 0.48 m in Y and 1.18 m in Z.

#### Terrain Extraction

**Table 3** provides an overview of the employed methods and resolution with a brief remark on the achieved results. All tested software suites include various terrain extraction algorithms. Choosing the appropriate methodology and optimizing input parameters remains a challenge and often relies on extensive user experience. Many of the terrain extraction modules do provide automatic blunder detection and noise reduction strategies, e.g., reverse matching, which show various degrees of performance and also add processing time.

#### **SocetGXP**

The terrain extraction module in SocetGXP enables flexible customization of algorithms that can be adapted to suit different data and terrain types. The presented experiments are based on standard strategies for the ATE and NGATE algorithms. Several different methods were tested whereby ATE adaptive strategy and NGATE low-contrast strategy performed best for the employed satellite imagery.

Two series of experiments were performed based on the image triplet and the stack of images over the area of Edinburgh of the Seven Seas. Series 1 covered the full island and produced DEMs with up to 1 m GSD using the ATE and NGATE strategy. ATE provided the best performance for a DEM with 2 m GSD (**Figure 8A**), which shows a good level of detail and contains less noise than the DEM generated by NGATE. A higher sampling rate at 1 m GSD based on NGATE showed an increased level of noise and artifacts but no improvements in the level of detail. Series 2 only covers the northwest area of the island around the settlement. ATE and NGATE were used to generate a DEM with 2 m and 1 m GSD based on the stack with 10 satellite scenes. No improvement could be observed. However, the computational effort rose sharply and the results showed higher levels of noise and artifacts compared to the models derived from the image triplet. For urban areas, NGATE showed better performance than ATE. The DSMs depicted in **Figure 9** have been generated using the ATE module with GSDs of 10, 5, 2 and 1 m with no filtering or corrections applied. In the 10 m resolution DSM, many artifacts can be identified at steep slopes (**Figure 9B**). The DSM with 2 m resolution (**Figure 9D**) provides a good trade-off showing low


#### TABLE 3 | Overview of Terrain matching experiments.

feart-08-00319 September 8, 2020 Time: 10:47 # 13

\* Average point sample distance.

noise level and a high level of detail. In contrast, models with 1 m GSD derived from the MVS configuration with 10 satellite images do not show improvements but appear to include more noise and artifacts (**Figures 9E,F**).

#### **Erdas imagine**

The terrain extraction modules in Erdas Imagine require a varying level of user input and customization. SGM and ATE draw on the results of the triangulation to validate the quality of the stereo-pairs and predefine suitable parameters for stereo matching. The hierarchical approach of using the mDEM as seed model while increasing the resolution stepwise proved to be a beneficial strategy in ATE and eATE regarding computational performance, blunders and noise level. The SGM algorithm did not require any seed DEM.

While the DEMs generated by ATE appeared to be noisy, the final results from eATE showed a good level of detail with very few artifacts. The eATE module includes a very flexible terrain extraction methodology but demands a detailed technical understanding in order to balance computational performance with output quality. Matching at the highest pyramidal level provided a dense DEM at 0.5 m (approx. 1 pixel) GSD but included a high level of noise and was computationally expensive. **Figure 8B** shows the DEM derived by eATE at 2m GSD. The SGM module performs a pre-assessment based on the results of the triangulation to evaluate the quality of the image alignment for SGM matching. For the image triplet only the stereo-pair acquired in November 2015 provided adequate alignment quality. Dropping the quality requirements would degrade the accuracy of the final product and increase the level of noise and artifacts. The best results from SGM was derived from the single stereo-pair (**Figure 8D**), while the combined approach using all three stereo-pairs included a significant amount of noise and blunders. Dense point clouds were generated at pyramidal level 0 (1-pixel sample distance) with approximately 280 Million points and pyramidal level 1 (2-pixel sample distance) with approximately 80 million points. Both point clouds required intensive cleaning and filtering. With an average point spacing of approx. 1 m and a lower noise level, the SGM point cloud at pyramidal level 1 provided an adequate trade-off between sampling density and performance.

#### **Geomatica OrthoEngine**

The workflow used in Geomatica OrthoEngine did not include a triangulation. This was justifiable based on the findings of the image selection process and the results of the triangulation shown in section "Spaceborne Triangulation." Only one stereopair from November 2015 was used to generate a highquality DEM at 2 m GSD using the NCC algorithm. Higher spatial resolutions, the use of the SGM algorithm or additional stereo-pairs did not provide any substantial improvements. The results show an almost complete DEM with little noise (**Figure 8C**).

### Point Clouds and DSMs From Drone Photogrammetry

An overview of the results of the drone survey is given in **Figure 7B**. The dense block of images covered the area around the harbor and a large number of perspective and panoramic views provided an almost complete overview of the surrounding mountain slopes. These additional images enabled the generation of a sparse point cloud of the neighboring areas during image matching. Due to the expected distortions and extrapolation of errors, this sparse point cloud of the mountain slopes was excluded in the final model.

**Figure 4** shows the overall configuration of the image block. The image acquisition in different flying heights and orientations, perspective views of specific objects, e.g., the harbor crane, and the panoramic shots created a geometrically stable block. **Table 4** summarizes the results of the aerial triangulation using different numbers of ground control and checkpoints.

FIGURE 8 | Raster DEM with 2m GSD generated with (A) BAE SocetGXP ATE, (B) Erdas Imagine eATE, (C) PCI Geomatica Orthoengine. (D) Point cloud generated by Erdas SGM with approx. 1 m point sampling distance.

Only CPs were introduced in the first triangulation. A RMSE of 1.09 m in Northing, 1.12 m in Easting and 191.39 m in height showed the expected positional performance in plain but indicated a significant shift in the heights between the image position recorded by the drone and the GCN. In the second triangulation, five points and in the final triangulation all ten points were used as full GCPs. Both attempts showed positional accuracies in the range of a few centimeters, whereby the height accuracy was slightly lower than for the horizontal components. The orthophoto in **Figure 10A** verifies the positional accuracy. 363 from 373 images were used for camera self-calibration and the cross-correlation matrix in **Figure 10D** shows low correlations between the camera calibration parameters.

Following the successful triangulation, the point cloud was generated using dense image matching. The quality of this point cloud greatly depends on the configuration and stability of the image block. **Figures 10A-C** give a detailed view of the harbor area. The point cloud shows a high level of detail, completeness and has a low level of noise. To avoid artifacts and distortions as seen in the top right corner of **Figure 10A,** the area of interest was restricted to well-captured areas of the harbor. Bordering regions were excluded from the final point cloud, DSM and orthophotos.

GSD derived from a MVS image stack.

TABLE 4 | Drone photogrammetry results of triangulation and camera calibration.


# Point Cloud Integration From Spaceborne and Drone Photogrammetry

Regardless of the differences in spatial resolution or point sampling, the point clouds from spaceborne and drone-based data appear to be visually well-aligned as is shown in **Figure 11A** which gives an overview of the entire island. **Figure 11B** focuses on the local area showing the merged point clouds derived from spaceborne and drone imagery. In this figure the point cloud from drone photogrammetry included the extended regions of the hill slopes (Section "Drone Photogrammetry"). A visual

investigation using cross-sections, cloud-to-cloud comparison and map overlays showed no obvious misalignments or rotations between spaceborne and drone-based models in the area covered by the drone survey **Figures 11D,E**. Although the sparse point cloud from the drone photogrammetry that covered the extended region, was immensely extrapolated, no significant distortions and wrapping effects were detected in comparison to the spaceborne model.

A closer examination of the harbor region showed a low level of remaining noise but poor object definition in the filtered models from the spaceborne photogrammetry (**Figure 11C**). No apparent differences in the level of detail between the SGM model at 1 m and 0.5 m GSD were noticed. Spaceborne and drone-based models aligned well and no apparent shift in co-registration was observed in plain (X and Y) and in height (Z). A cloud-to-cloud comparison using the ICP algorithm between the drone-based (reference) and spaceborne point clouds showed an average distance of 0.26 m with a standard deviation of 0.68 m. Despite the differences in resolutions and object definitions, this result confirms a good fit between both data sets. An additional cloud-to-cloud registration between both data sets derived a translation vector of −0.32 m in X, 0.09 m in Y and 0.25 m in Z, without any rotations. This reduced the average distance between both data sets to 0.14 m with a standard deviation of 0.53 m. Furthermore, the final comparison of the Z components shows a mean difference of 0.01 m with a standard deviation of 0.22 m. Due to its higher positional accuracy, the drone-based point cloud served as reference point cloud. Since only 3 small translations and no scales or rotations were applied, extrapolation of errors toward the full extent of the final model were avoided.

# DISCUSSION

In the absence of accurate topographic data derived from dedicated airborne mapping sensors, VHR satellite images and low-altitude drone images provide alternative data sources to generate topographic mapping data. For the presented study, a heterogeneous data set consistent of VHR imagery, low-altitude drone imagery as well as local ground control was assembled over the remote island of Tristan da Cunha with the aim to produce a surface representation with high fidelity and geometric resolution. Conventional photogrammetric workflows implemented on widely available software tools ensure reproducible results and allows to adopt this workflow to other study sites. Due to the differences in image scale and extent, both data sets where processed separately and subsequently merged into a combined point cloud. In the presented case, geodetic ground control provided sub-centimeter absolute positional accuracy for the central area of interest. This ground control was used to georeference and assess both datasets. In the absence of ground control of such accurracy, the dataset with the expected highest positional accuracy should be used to reference the other datasets.

FIGURE 11 | (A) Overview merged point clouds. (B) Area around Edinburgh of the Seven Seas. (C) Perspective view harbor area point cloud from SGM. (D) Top view of the final cloud around the harbor area. (E) Perspective view harbor area blended point clouds.

# Satellite Photogrammetry

Satellite operators have collected large archives of VHR satellite images with spatial resolution better than 1 m. These data provide an inexpensive source of imagery to generate high-resolution DEMs with 1–2 m GSD over most regions of our planet. The latest generation of highly agile VHR image satellites, e.g., WorldView3, are capable of collecting accurately georefenced images with up to 0.30 m nominal GSD in the panchromatic channel. The image archive of DigitalGlobe/Maxar comprised 100 Petabyte of images in 2016 with its fleet of satellites currently adding approximately 80 Terra Byte of image data every day. With many new satellite constellations in preparation, imagery will be available from a number of suppliers. Nevertheless, the selection and quality assessment of suitable image pairs and preprocessing is a crucial step for a successful terrain extraction. In this case study an image triplet was identified, which showed ideal conditions for the terrain extraction. The selected images provided full coverage of the island, good geometric conditions with favorable base/height ratios, high radiometric image quality without apparent atmospheric effects, no cloud coverage and low temporal decorrelation, e.g., changes in land cover.

While we detected noteworthy misalignments and shifts in the orientation of older satellite scenes, e.g., for QuickBird2, the error in the direct georeference of recent WorldView2 and WorldView3 images used in this study was detected to be below 0.5 m in Easting and Northing. This was much better as in the specification provided in DigitalGlobe/Maxar technical references, which state the absolute geolocation accuracy for a nadir image to be at 5 m circular error in 90% of cases (Digital Globe, 2014, 2016), but they are also in line with the evaluation by Poli et al. (2015). Following a consequent photogrammetric workflow, an aerial triangulation was performed to optimize the alignment between the satellite scenes. The relative alignment between stereo-pairs has a significant impact on the performance of the stereo matching algorithms and thus on the quality of the extracted DEMs.

The use of ground control needs to be considered carefully. If the reliability and accuracy of GCPs is questionable and a high a priori quality of the image orientation via direct georeference is expected, it is advisable to perform a triangulation based on tie points only. The expected accuracy of the direct georeference depends on the origin and age of the image data and should be provided as metadata. Furthermore, many scientific Calibration/Validation studies may also provide an independent assessment of VHR imagery (Poli et al., 2009; Poli and Toutin, 2010; Reinartz et al., 2010; Dowman et al., 2012; Perko et al., 2018). Since ground truth information was limited to a small area in the center of the study site, only one GCP was introduced to the spaceborne triangulation. This verified the accuracy of the direct georeference of the VHR imagery and allowed for small horizontal (or planimetric) as well as vertical corrections. The

geometric and radiometric quality of the VHR images acquired by the WorldView constellations was investigated on established test sites by Angiuli et al. (2010) and Poli et al. (2015). The study confirms the accuracy of the direct georeference found in previous studies on a remote location, which adds confidence in the quality of the data source.

The assessment of the different terrain extraction modules showed that all software packages include suitable algorithms to derive high-quality DEMs. Especially the use of the NCC algorithm produced high-quality DEMs at 2 m GSD with low noise level and few blunders. Furthermore, the use of an improved 'seed' DEMs increased the efficiency of the terrain extraction process. The Tridicon SGM implementation was more computationally efficient and produced a denser point cloud at approximately 1 m point spacing, but also included a higher level of noise and outliers. A brute force approach based on a MVS configuration with 10 VHR images showed a higher level of noise but no noteworthy improvements at a much higher computational expense. This can be related to temporal decorrelation, lower radiometric quality, atmospheric effects, less favorable image geometry and shortcomings in the image alignment. In a successful study Loghin et al. (2020) investigated the potential and limitations of small object presentations in DSMs derived from purposely collected VHR imagery.

# Drone Photogrammetry

Capabilities and endurance of small and portable multirotor drones as used in this study are still limited but are expanding all the time. The latest generation of semi-professional systems have improved positional sensors, cameras, extended battery capacity and use more sophisticated flight management software, which enable the capture of significantly larger areas. Flight management and planning software enable reliable automated missions in predefined standard configuration comparable to conventional photogrammetric blocks. Free flight missions allow the capture of bespoke image blocks but are considered difficult and do require skilled operators. The limited stability of small cameras (Gerke and Przybilla, 2016) does require frequent recalibration or self-calibration as a part of the aerial triangulation. Following the recommendation of Fraser (2018), a suitable photogrammetric block requires redundant images captured with significant overlaps, different perspective views, different flight heights, as well as a significant number of welldistributed GCPs.

In this study, it was challenging to deploy the DJI Phantom 3 drone under the harsh weather conditions. Only 2 flight missions in manual flight mode could be performed during a two weeks mission. Still, the drone survey provided a geometrically stable block of images that captured the main area of interest from a range of viewing angles and provided perspective views of the overall scene.

The absolute positioning accuracy using only the low-grade GNSS sensors onboard the drone showed good performance in Easting and Northing but had a significant bias in heights. Such height offsets are not uncommon issues when using consumergrade drones. Inaccurate or misinterpreted information from image metadata (EXIF headers) seems to be a trivial but frequent problem, which has been reported in many discussions and investigations (Griffiths and Burningham, 2018; Bertin et al., 2020). Thus, accurate ground control information is still indispensable in order to prevent such errors. Furthermore, no boresight calibration between the positional sensor and camera position was applied. Many studies have investigated the influence of the accuracy and distribution of GCPs (Tonkin and Midgley, 2016; James et al., 2017; Martínez-Carricondo et al., 2018; Rangel et al., 2018) which highlight the importance of well-distributed and accurately determined GCPs. Ten evenly distributed GCPs were derived from the geodetic control network assure the accuracy of the model both in absolute positioning and the avoidance of distortions of the 3D model.

An increasing number of computer vision-based processing software suits offer highly automated post-processing with the generation of dense point clouds or mesh models. Many of these operate as 'black box' tools without the possibility to modify or customize workflows. Moreover, they provide limited information about the quality of the generated model. Software tools should provide detailed processing logs and reports which describe the results and with quantitative quality indicators. The drone photogrammetry package Pix4D used in this study provided a comprehensive quality report of the aerial triangulation, camera self-calibration estimations on the achieved relative and absolute positional accuracy of the model. Based on a strong block geometry with images captured from several flying heights and perspectives, stable and uncorrelated camera calibration parameters were obtained for both flights.

The dense point cloud from the drone photogrammetry only covered a small area of interest whereby a sparse point cloud generated based on the perspective images, provided a reasonable representation of the nearby mountain slopes which was unexpected and unreported. Based on poor geometric conditions, such point clouds are likely to include large distortions through error extrapolations. They are useful for visualizations but should be discarded for analytical applications.

# Point Cloud Integration and Validation

Since the spaceborne and the low-altitude drone-based models were accurately georeferenced using the geodetic ground control, both models were merged based on their geolocation without the use of natural targets or cloud to cloud registration. Due to the difference in spatial resolution and extent, the application of a best fit cloud-to-cloud registration via an iterative closest point (ICP) algorithm requires additional thoughts to avoid the extrapolation of errors on from the local area to larger extent one hand or degradation of the positional accuracy on the other hand. Detailed metadata about the data acquisition methods, spatial resolution, relative accuracy as well as absolute positional accuracies are required to guide this process.

Filtering and cleaning remain to be an essential part of the DEM extraction process and have a significant impact on the results of any subsequent data analysis. Ideally, this step should be performed by the informed user of the DEM, who is best aware of the requirements of the intended application.

# CONCLUSION

feart-08-00319 September 8, 2020 Time: 10:47 # 19

In this study, we used archived VHR satellite images and low altitude drone imagery to produce a high-resolution DEM over the remote island of Tristan da Cunha. The DEMs derived from VHR imagery cover the whole island and were generated following a rigorous photogrammetric workflow. We assessed the performance of the terrain extraction and dense image matching algorithms included in three widely available geospatial software packages. For the lowaltitude drone survey, we used the drone photogrammetry package Pix4D to derive a dense point cloud of the harbor area and research installations with up to 500 points/m<sup>2</sup> . A geodetic survey provided ground control with sub-centimeter absolute positional accuracy in the same area. A cloudto-cloud comparison of spaceborne and drone-based point clouds showed a tight fit based on their georeference. The co-registration of both models was optimized via small translations in plain and height and subsequently merged to a single point cloud.

We have barely touched the potential of the use of archived VHR satellite images. The high-resolution DEMs were derived from well know conventional stereophotogrammetric methods. However, image selection and preprocessing are crucial. Good-quality DEMs can be generated if suitable stereo-pairs or image triplets are available from archived data. We demonstrated the application of best practice photogrammetric workflows on widely available geospatial software suits. Multiview image matching using selected combinations of stereo-pairs should further improve the completeness and accuracy of the models in occluded regions and help to eliminate blunders. Considering a large number of available images, multi-temporal change detection, e.g., for landslides, should be possible.

Although the positional and navigational sensors of small drones have improved considerably, the use of well-distributed GCPs is still indispensable to assure geometric fidelity and the accurate georeference. Drone surveys using low-cost semiprofessional drones should follow best practice recommendations for a stable photogrammetric block.

By applying a rigorous georeferenced workflow with well established and accurate control points, merging the DEMs in point clouds representation was straight forward despite the differences in geometric resolution and extent. In the absence of accurate ground control, the dataset with the expected highest positional accuracy should provide the reference to coregister all other datasets. The lack of suitable reference data prevented a thorough quantitative assessment, especially of the spaceborne DEMs. However, this is not uncommon for projects of this kind. The collection of extended ground control and reference data will be highly desirable tasks for future missions.

The study showed that archived VHR satellite imagery and low-altitude drone imagery provides two alternative data sources to generate high-resolution DEMs over regions where airborne mapping data are not available. In the presented study, we combined the strength and weaknesses of both data sources to produce a high-resolution point cloud that covers the full extent of the study site and preserves fine detail of a local area of interest. DTMs as well as orthophotos, can easily be derived from the final point cloud.

The growing archive of VHR images collected by current and future satellite fleets will allow the generation of multi-temporal high-resolution DEMs over most parts of the world and enable dedicated geo-monitoring tasks. Low-cost and simple to use drones are already a useful tool to collect high-resolution data sets over small areas. We demonstrated that both data sources combined can provide high-quality DEMs over remote regions with poor mapping.

# DATA AVAILABILITY STATEMENT

The extended datasets for this study are available on request to the corresponding author. The final 3D point clouds as well as 2.5 raster DEMs are available on the Pangea data repository with the article's doi: 10.3389/feart.2020.00319.

# AUTHOR CONTRIBUTIONS

DB was responsible for the analysis of both the spaceborne and aerial photogrammetric data sets over Tristan da Cunha while FT was responsible for all installations on the island, the geodetic surveys and their analysis, the drone imagery acquisition, and the project coordination overall. Both authors contributed to the article and approved the submitted version.

# FUNDING

DB is funded by the University of Luxembourg. The DigitalGlobe Foundation provided the satellite imagery under the umbrella of Image Grand. The installations and survey work on Tristan da Cunha were funded through projects at the University of Luxembourg (in particular SGSL) and the United Kingdom National Oceanography Centre.

# ACKNOWLEDGMENTS

We acknowledge Hexagon Geospatial for providing access to and advice on Erdas Imagine as well as additional software licenses for the automatic terrain extraction modules, PCI Geomatics for providing Geomatica OrthoEngine software as well as professional expertise and guidance, BAE Systems for temporary access to SocketGXP, and the help and support from the community, various contract workers and local administration of Tristan da Cunha, the contributions from the IGN, France, with respect to the DORIS station and the IGS for data and satellite products. Finally we would like to thank the reviewers for their constructive input to this study.

# REFERENCES

feart-08-00319 September 8, 2020 Time: 10:47 # 20


Part 1 International Archives of Photogrammetry and Remote Sensing - ISSN: 1682-1777, Vol. XXXVIII, Calgary, AB.


**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Backes and Teferle. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.