We would like to thank the reviewers for their constructive review. Special acknowledgement to Ann Grant for her valuable remarks and corrections concerning English version.
1Digital Elevation Models (DEMs) are produced from various techniques such as photogrammetry from aerial photos and satellite images (Ayhan et al., 2006; Paul et al., 2004; Casson et al., 2003), radar images (AbuBakr et al., 2013; Hong et al., 2010), laser scanning using airborne Light Detection and Ranging (LiDAR; Allen et al., 2012) and classical ground survey methods (Nelson et al., 2009; Hirt et al., 2010). Their precision and accuracy will depend on the sensor used, technique, scale and type of terrain (Florinsky, 2012). Detailed information about the earth’s surface is increasingly required in order to perform simulations and calculations from the terrain that respond to problems such as landslides and floods. In Mexico, this detailed information is not yet available for all the country, the number of check point data bases is limited, and most of the bare earth surface models are generated by the use of contour line interpolations such as in ARCGIS, IDRISI and ILWIS software (Hutchinson, 1988; Douglas, 2000). DEMs provided by INEGI (Instituto Nacional de Estadística Geografía e Informática) follow this procedure and cover most parts of the country. Vector data can be used to generate DEMs, but the nature of the terrain of fluvial and coastal plains requires other types of interventions to produce detailed results.
2Generally, the interpolation of contour lines, a phenomenon inherent to raster contour line drawing produces an overestimation of altitude values corresponding to these contour lines in the resulting images (Taud et al., 1999). This problem can be partially solved by changing the hypsometric scale. A β spline smoothing applied to altitude values in millimetres attempts to diminish this overestimation (Parrot, 2014a). Furthermore, this millimetre hypsometric scale eliminates the staircase effect observed in wide zones that present a slight slope. This improvement is not so necessary in mountainous regions because of their high number of closely spaced contour lines, but in coastal plains the problem remains and needs another type of solution. For instance, in the Coatzacoalcos lower basin (Veracruz State, Mexico) used here as a test area (fig. 1), 2889.027 km2 lies at between 0 and 10 m, i.e. for 35.5% of the lower basin there are no details; 92% of the 0-10 m zone has a slope of 0°.
Fig. 1 – Area of Veracruz State, Mexico, used as the basis of a flood simulation.
Fig. 1 – Zone d’étude dans l'Etat de Véracruz, Mexique, utilisée comme base de la simulation hydrologique.
DEM interpolation (Parrot, 2012a) from dxf contour lines (INEGI, 2012).
Interpolation du Modèle Numérique de Terrain (Parrot, 2012a) à partir d’un fichier dxf contenant les courbes de niveau (INEGI, 2012).
3Here, we propose an alternative method for producing a DEM for coastal and flood plains by means of vector data. First, the treatments used to generate this DEM are described. Then, the result is validated by the root mean square roughness (RMSR) and the root mean square error (RMSE) using the Felicísimo (1994) method, a LiDAR Digital Terrain Model (DTM) and a regional flood simulation to compare the resulting areas and volumes due to flooding since this region is regularly affected by flooding. The importance of this product lies in its potential for use as a valuable tool for civil protection.
4The test area used to illustrate the method is in southeastern Mexico, in the lower basin of the Coatzacoalcos River fluvial plain in Veracruz State, where the continent narrows into the Isthmus of Tehuantepec (Oaxaca and Veracruz States). The river rises in the Sierra Niltepec, Oaxaca State, and flows through 325 km from south to north until arriving at the Gulf of Mexico. The basin has an area of 17,369 km2 and the runoff is estimated as 28,093 Mm3 (Comisión Nacional del Agua, 2010). This area is regularly affected by floods and is, therefore, ideal for simulation of such an event as a validation tool.
5The vector format is a tagged data representation of all the information contained in a digital drawing file; each drawing element in the file is preceded by an integer number that represents a group code that can be modelled directly using a GIS. Even if nowadays numerical DEMs can be generated by stereoscopic digital images (Baily et al., 2003; Schiefer and Gilbert, 2007; Publighe and Fava, 2013), DEM building of raster type from vector contour lines remains a low-cost treatment that provides good results and that can be improved as we show in this paper. The contour lines used in vector format were digitized from stereo pairs of aerial photographs so that the available scales of these data are 1:50,000 or 1:250,000; the DEMs in raster format were obtained mainly from the interpolation of vector contour lines providing resolutions of 15, 30, 60, 90 and 120 m (INEGI, 2014). The resulting DEMs were verified by level banks and geodetic points at the national scale. The resulting sets of contour lines provided by INEGI (2012) do not take into account the contour line corresponding to sea level. This value has to be sought in another data type that contains information about water bodies. The other water bodies, such as lagoons, lakes and dams, have no altitude value. These perimeters need to be integrated in the set of contour lines to attribute, in a following stage, to each pixel that belongs to the sea level or water body an altitude value before final interpolation. The present DEM generation method requires six stages (fig. 3) as follows:
Fig. 2 – Altitude ranges.
Fig. 2 – Classes d'altitude.
A: Superimposed areas on the shadowed DEM: dark gray, 0-10 m; light gray, 11-20 m. B: Surface percentage of hypsometric classes.
A : Tranches d’altitude reportées sur un estompage du Modèle Numérique de Terrain : gris foncé : 0-10 m, gris clair : 11-20 m. B : Histogramme du pourcentage des surfaces correspondant aux tranches d’altitude.
Fig. 3 – Flow chart for DEM production for coastal and fluvial plains.
Fig. 3 – Organigramme de la production de MNT dans les plaines côtières et plaines alluviales.
Software and algorithms used: Parrot (1998, 2005, 2006, 2011, 2012b-g, 2013, 2014a, 2014b) and Parrot and Ramírez-Núñez (2012a-c).
Logiciels et algorithmes utilisés : Parrot (1998, 2011, 2012b-g, 2013, 2014b) and Parrot and Ramírez-Núñez (2012a-c).
6As oceans are classed as water bodies, the first vector contour line for coastal areas generally starts at 10 m altitude in INEGI data. The vector water bodies file contains labelled objects such as river banks and perimeter of the sea surface considered as polygons within the topographic maps. This file also contains channels and aqueducts but none of these features has an altitude value. A first algorithm (Sum_cn_ha_dxf) attributes the altitude value 0 to coastlines and introduces this feature into the vector contour line file. As the river banks, even if they are inland, have the same label, the altitude of all river banks is 0. Moreover, oxbow and lagoon borders also have this altitude. When the vector contour line file has been completed in this way, a second algorithm generates a raster image of 8 bits (Transf_dxf_v2) where a grey tone value is given to each contour line and water-body border according to its altitude (Fig. 4A). Artifacts inherent to raster format such as hiatus and 4-neighbour pixels are removed (Net_curve2 and Hiatus).
Fig. 4 – Basic raster images for DEM generation.
Fig. 4 – Images raster de base utilisées pour engendrer le MNT.
A: Contour lines and water bodies (8 bits). B: Surface of the rivers. C: River banks. D: Water bodies. E: Perimeter of the water bodies. F: Remaining surface inside the flood plain (< 10 m). Note: Figures A, C and E correspond to a dilation (one iteration) to emphasize the contours.
A : Courbes de niveau et étendues d’eau (8 bits). B : Surface couverte par les rivières. C : Périmètre des berges. D : Plans d'eau. E : Périmètre des plans d'eau. F : Surface restante à l'intérieur de la plaine inondable (< 10 m). Note : les contours des images A, C et D ont été dilatés pour des raisons de lisibilité.
7All the water body limits are codified with the value zero and this value is assigned manually to the inner surface. This manual operation allows all the water bodies to be labelled. At this stage, different treatments can be used specifically for water bodies that do not correspond to river surfaces, such as oxbows, lagoons and dams. It is possible to define an altitude value for each labelled water body in order to generate the DEM, or to define two different groups: the first corresponds to the water bodies that belong to the fluvial plain (an altitude lower than 10 m) and the second group concerns water bodies with an altitude above 10 m. At the end of these pre-treatments, five binary images are obtained (Parrot 2011, 2013), images that represent the base for DEM generation (see flowchart in figure 3). These images are as follows: river surface (fig. 4B), river banks (fig. 4C), water bodies in the fluvial plain (fig. 4D), water body perimeters (fig. 4E) and the remaining surface inside the flood plain considered as zones with an altitude lower than 10 m (fig. 4F).
8The raster file that contains the grey-tone contour lines (8 bits per pixel) and its corresponding table of grey tone/altitude are used to produce an image that contains contour lines employing a 4-bytes codification of its constitutive pixels in such a way that it is possible to register altitudinal values in metres, decimetres, centimetres or millimetres (Brod4_mx and Crear_tabla). This image will be used to generate the DEM at the final stage after integrating the results provided by the treatment described as Stage 4. The contour lines are reported in an integer image where the background corresponds to the null value -999999, and this provisional image is named PI_C4 (fig. 3).
9This process is based on the first four binary images: river surface, river banks, water bodies inside the flood plain and the image of the water body perimeters. The main purpose of the developed algorithm is to calculate the altitude values of all the water surfaces by means of two different approaches.
10The first approach needs to separate the isolated water bodies and river surfaces in order to do a specific treatment that concerns the rivers. In that case, the altitude values of the river banks are calculated separately before a river surface is interpolated. The procedure consists in following each river bank from the lower to the upper point and in assigning to each pixel that belongs to the river banks its altitude value according to the number of steps and the altitude value of the two river bank extremities. Regarding the isolated water bodies, the above-mentioned labelling allows attribution separately to each water surface its own value of altitude (Ramírez-Núñez and Parrot, 2014).
11In practice, this procedure is time consuming and therefore limited for a wider use on large areas; it requires many specific and local interventions in order to eliminate some incoherent results especially when confluents are encountered. For this reason, a simpler treatment has been defined and employed. In this case, the regional water level representing the flood expanse is considered as a tilted plan used as a hypsometric reference to calculate the elevation of the river surface and the altitude of the isolated water body surfaces. The number of artifacts generated by this method depends on the configuration of the drainage network and the lack of preferred orientation. In the test area, the drainage network runs globally from south to north, and only a few meanders may cause a slight change in its orientation. Regarding the water surface isolated, the developed algorithm (RiverBodies) allows definition of a reference point used to calculate the different water body surfaces that can be either the centre of gravity of the surface, the point of the highest elevation or the point of the lowest altitude.
12This partial result (fig. 5) is an image called water surface altitude [WSA] codified by using 4-byte pixels.
Fig. 5 – Recalculated altitudes of the surfaces of the water bodies.
Fig. 5 – Altitude recalculée des étendues d’eau.
13The data used to generate the riverbed altitudes are the lower altitude contour line image that in a general way limits the floodplain (fig. 4A), the image corresponding to the water surface (WSA), the image of the river banks (fig. 4C), the binary image of the water body perimeters (fig. 4E) and the binary surface of the flood plain between the altitudes 0 and 10 m (fig. 4F). As for the final interpolation, the treatment is as follows. For each pixel located in the flood plain, the algorithm New_fast_cauce (Parrot and Ramírez-Núñez, 2012b) searches the smaller distance ds between this pixel and a pixel of the contour line that corresponds to the upper limit of the flood plain (here, the altitude As of this contour line is 10 metres, 100 decimetres, 1000 centimetres or 10000 millimetres) and the smaller distance di between this pixel and the river bank or a water body perimeter (Fig. 6). The corresponding altitude Ai depends on the altitude of the water body reached at this point. Then the altitude value AP of the studied pixel is derived from Ap = Ai + [(As – Ai) x (di/d)] where d = di + ds.
14The result of the floodplain altitudes (FPA) is a file using 4-byte pixels (fig. 7).
Fig. 6 – Example of interpolation.
Fig. 6 – Exemple d'interpolation.
As is the upper altitude (here always 10 m); Ai is the lower altitude of the nearest water body perimeter; ds is the distance As-Ap; di is the distance Ap-Ai; and Ap is the resulting altitude of the pixel.
As correspond à l’altitude supérieure (dans le cas présent celle-ci est toujours égale à 10 m) ; Ai correspond à l'altitude inférieure dont la valeur dépend de celle du périmètre de la masse d’eau la plus proche ; ds correspond à la distance jusqu’à As et di à la distance jusqu’à Ai ; Ap est l'altitude résultant de l’interpolation pour le pixel étudié.
Fig. 7 – Recalculated altitudes on the flood plain between 0 and 10 m without taking into account the altitude of the water bodies.
Fig. 7 – Altitude recalculée de la plaine inondable entre 0 et 10 m, sans tenir compte de l’altitude des étendues d’eau.
Grey-tones image of the FPA integer 4 file.
L’image FPA est une image de 32 bits en teintes de gris.
15The water surface DEM (WSA) that corresponds to the interpolated rivers and water bodies, together with the former FPA, is superposed (Superpos_rio) on the provisional image PI_C4 that contains the contour lines. The altitude AP of the pixels corresponding to the remaining null values is calculated by the linear interpolation defined above, but in this case the algorithm computes it in each altitude layer (Newmiel, Miel4_mx). An altitude layer corresponds to the hypsometric interval between two consecutive contour lines that are used to define ds and di as well as Ai and As.
16Finally, a smoothing treatment is applied outside the water surfaces used as a mask. The smoothing integrated in the program Suav_mask (Parrot, 2012f) is based on a β spline function. The resulting DEM has a 5 metre resolution and a hypsometric scale in centimetres.
17The accuracy of a DEM depends on factors such as the type of data, interpolation, resolution and sensor used to generate it. Low- and medium-frequency errors are related to DEMs produced by means of contour line interpolation (Rieger, 1996; Florinsky, 1998, 2012; Aguilar et al., 2005; Ghilani and Wolf, 2008). The RMSE of elevation is a general way to estimate the accuracy of a DEM, but this evaluation depends on the data from which it was generated (Weschler, 1999). Many estimators use reference altitude points (Ivanov and Kruzhkov, 1992; Bolstad and Stowe, 1994; Wechsler, 1999) or morphometric variables (Young, 1978; Evans, 1979). Generally, validation criteria consist of a comparison between the resulting DEM and another more accurate surface (Wood, 1996; Wechsler, 2000). Nevertheless, the number of reference altitude points is generally limited and it is difficult to consider the reference DEM as the real surface (Florinsky, 2012).
18In this research, in a first approach, we assess the accuracy of our DEM by calculating the RMSR and defining a criterion of homogeneity (H). This index corresponds to the difference between 100 and the quotient multiplied by 100 of the average of the differences between the values of the roughness computed by lines, columns and whole image and the average of these values (see later).
19The RMSE is calculated according to the Felicísimo (1994) procedure, which compares each pixel of the DEM with a local surface provided by its environment. According to Florinsky (2012), this procedure is an elegant validation approach that does not require field data, the point number of which is in any way limited (Wood, 1996; Brasington et al., 2000; Heritage et al., 2009). The results for the RMSE and the RMSR appear in Table 1. Furthermore, in order to overcome the lack of actual field data, the validation is also based on a comparison involving a LiDAR DTM which has the same horizontal and vertical resolution as the DEM generated by the technique described above.
20The last validation corresponds to an application of a regional flood simulation using the procedure described above and applied to the resulting DEM and a new LiDAR product that gives detailed information at a local level.
21The RMSR is the square root of the sum of the squares of the difference of altitude between each pixel (i, j) of an image of size m × n and the mean value μ calculated as follows.
22The RMSR can be calculated according to the lines (Rl), columns (Rc) or image (Rt). The standard deviation of the whole image is calculated and a coefficient of Homogeneity (H) is proposed and corresponds to:
23The program takes into account the absolute value of the differences between Rl, Rc and Rt. The same calculation can be applied to a morphologic variable such as the slope (tab. 1).
Tab. 1 – RMSR calculated for the generated DEM and the morphometric variable slope, taking into account the Felicísimo (1994) approach.
Tab. 1 – Calcul du RMSR appliqué au MNT et sur une variable morphométrique (pente) calculée en tenant compte du traitement proposé par Felicísimo (1994).
24As mentioned above, methods of accuracy estimation are based on the analysis of the differences observed between altitude values provided by two different sources. The estimator RMSE is calculated as follows:
25where yi is an elevation point from the resulting DEM, yj the value of the corresponding point on the “reference” surface and N the number of sample points. Felicísimo (1994) proposed that a “reference” surface be generated by considering the hypsometric values of the four cardinal points of each studied pixel. As this process is applied to all the DEM pixels, it is also possible to calculate the arithmetic mean and the standard deviation of these differences.
26As the last validation concerns a comparison between the simulations of the flooding expansion calculated using respectively the generated DEM and a Digital Terrain Model provided by LiDAR data (INEGI, 2013), the same treatment has been applied to both models. It is also possible to compare the two roughness results in a treatment done line by line (fig. 8). The RMSE and the standard deviation of the generated DEM are lower (0.0074 and 0.0059) than for the LiDAR DTM (0.0945 and 0.0793), showing that the latter model has a greater level of errors in relation to local manual corrections.
27Finally, a comparison between the two digital models needs to use the fast Fourier transform (Baudemont, 1999) in order to smooth the meander scrollbars (fig. 9).
Fig. 8 – RMSR comparison between the generated DEM and the LiDAR DTM.
Fig. 8 – Résultats pour le RMSR (root mean square roughness) des deux modèles numériques de terrain.
Fig. 9 – Comparison between the generated DEM and the result of the application of the fast Fourier transform to the LiDAR DTM.
Fig. 9 – Comparaison entre MNT et la transformée de Fourier appliquée au MNT LiDAR.
A: LiDAR DTM. B: FFT application (frequency 5% elimination). C: Generated DEM.
A : MNT LiDAR. B : Traitement de la transformée rapide de Fourier. C : MNT.
28A regional flood simulation has taken into account the resulting DEM and a LiDAR DTM in order to compare and so validate the accuracy in terms of volume and surface. In both cases DEMs have a 5 m pixel resolution and a hypsometric scale in centimetres. In accordance with the height reached during a regional flood event, it is possible to calculate the water surface and the related water volume.
29In a first approach, the flood is simulated from the intersection between the DEM surface and the upper surface of a water sheet. This water sheet is defined by considering the local increasing value of the water surface due to the studied or simulated event upstream and downstream of the main drainage network of the training area. It has also been calculated by use of a local Gaussian function in order to follow the displacement of the general water wave (program Gaussian_Lateral_Flooding; Parrot, 2014c).
30The area and volume were calculated considering a maximum water level rise of 2 m above the river surface; this corresponds to field measurements done after a flood in 2010 and the data reported by the inhabitants about former floods.
31The area calculated from our DEM of the flood plain was 716.63 km2, whereas the LiDAR DTM was from 787.82 km2. Figure 10 shows the overlap between the two DEMs as well as those areas calculated by the individual models.
Fig. 10 – Flooded areas in the lower basin of the Coatzacoalcos River for the two types of DEM used.
Fig. 10 – Les zones inondées dans le bassin inférieur de la fleuve Coatzacoalcos pour les deux types de DEM utilisés.
Black: LiDAR calculation; Light gray: area solely covered by the generated DEM; Dark gray: area common to the two calculations, LiDAR DTM and DEM.
Noir : surface seulement couverte par LiDAR DTM ; gris clair : surface uniquement couverte par le MNT provenant du traitement ; gris foncé : calcul de la surface commune aux deux modèles digitaux (MNT créé et LiDAR DTM).
32Such surfaces correspond respectively to a volume of 1.02 km3 for the DEM and of 0.86 km3 for the LiDAR DTM. This difference is clearly associated with the higher level of roughness of the LiDAR DTM due to the presence of many scrollbars associated with meandering, whereas the algorithm used to generate the DEM smoothes the riverbed bottom.
33As floods represent one of the most frequent risks in Mexico, and detailed DEMs such as LiDAR for flood plains and coastal areas are not yet available for all parts of the country, we propose the production of an adaptive raster DEM taking into account vector data at the 1:50,000 scale. The multidirectional interpolation method used here has been reported to give good results (Parrot, 1998; Pérez-Vega and Mas, 2009). This constructive DEM can be used to provide a quick response for governmental agencies in case of emergency.
34Until the generation of LiDAR DTMs that will cover all Mexican regions, especially in flood and coastal plains, the DEM obtained by the proposed method provides reliable results in general calculations such as simulation of regional flooding, for which the detailed scale of LiDAR is not necessary.
35Moreover, for the moment, LiDAR data for flood-plain and coastal areas need corrections that concern mainly the river surface. LiDAR systems use two frequencies, and water depth is determined by measuring the time delay to receive the return signal from the seafloor and water surfaces. When the infrared is reflected from the sea surface the higher-frequency green laser penetrates through the water depending on its clarity. The accuracy of the result depends on the correct management of these two responses. Actually, in the studied region abnormal measurement concerning the river surface indicates altitudes lower than the real values, for instance -0.54 m to the south of Minatitlan City which is at 38 km from the river mouth. A special treatment is required for the surface of a river under a forest canopy (Maune, 2007).
36These observations show clearly that the method developed here to achieve an accurate representation of the relief on coastal plains is an efficient tool for use at a regional scale.