Keywords

1 Introduction

Land planning and management requires an accurate and rigorous knowledge of the area [1,2,3]. The complexity of urban or natural scenarios is an issue that is increasingly driving the use of a 3D approach over the traditional 2D approach [4].

3D vision technologies in recent years are becoming more and more popular, the 3D representation of the territory and the elements that feature it is a very debated issue in science and is very actual [5]. The demand comes from different sectors: from industry to government for state and local administrations [6]. 3D maps of the territory offer significant advantages over a traditional 2D representation; the optimized and detailed three-dimensional representation of the territory, infrastructure and urban areas can be used for the analysis of morpho-evolutionary phenomena and for hydraulic analysis to assess the impact on anthropic activities, in addition to the provision of a higher quality visualization of information for planning applications [7, 8].

Three-dimensional reconstruction of land and anthropized territory requires the use of advanced and integrated processing techniques applied to many different types of geospatial data [9,10,11], which must be sourced by experienced users and involve quite significant investment, besides time consuming, especially if the areas to be modeled are large [12].

In recent years some government agencies have made available, free of charge, a few geospatial data. The resolution to make these data available arose to meet mainly public but also private needs. Among these are the LiDAR (Light Detection and Ranging) data acquired by aerial platform, high density data whose altimetric accuracy is on the order of the decimeter. The European Commission has recently published a report with the aim of disclosing all useful information on non-commercial LiDAR data made available by eighteen European Union countries, including Italy [13]. In Italy, LiDAR data have been acquired on behalf of the “Ministero dell'Ambiente e della Tutela del Territorio e del Mare” (MATTM) as part of the “Piano Straordinario di Telerilevamento” (PST-A) (Law 179 of 31/07/2002 art. 27) and are provided on demand [14].

The availability of data has led to a sudden increase in experimentation by both public and private entities for the development of new methods for the 3D modeling of the territory and the anthropic environment starting from data coming from remote sensing techniques [2, 15].

The methods proposed in the scientific community for the 3D reconstruction of the built environment can be grouped into two different macro-groups [11]: (i) methods that involve the integration of multiple data sources, typically both Airborne Laser Scanner (ALS) data and data acquired from 2D vector and raster cartography [16, 17]; (ii) methods that involve the use of only point clouds from ALS [18]. These latter methods are mainly constrained by the density of the point cloud, which mainly affects the extraction and vectorization of each single structure outlines. Differently, the methods based on the use of multiple data sources are more versatile and, in some respects, more accurate, since they are not influenced by the density of the LiDAR point cloud; the planimetric shape of the object is inferred directly from the 2D vector or raster cartographic information and subsequently extruded through the altimetric information coming from the data.

The methods based on the integration of different datasets for 3D modeling needs meeting strict requirements: the different spatial data, in addition to being framed in the same reference system, must also be processed to be aligned.

The main difficulties come up if the data are referred to different datums and there is no certain information about which is the reference datum and which is the target one; in these scenarios the error in the transformation may not result as acceptable [19]. These incongruences are mainly due to the different nature of the data, derived by means of different techniques, as well as their accuracy [20, 21].

Another example of open-source geo-spatial data are the National Topographic Data Base (NTDB) derived from a base cartography. The NTDB, produced by Public Administrations for the purposes of planning and land management according to the Content Specifications of the Catalogue of Spatial Data annexed to the Ministerial Decree - November 10, 2011 “Technical rules for the definition of the specifications of the Geotopographic database content”, follows the General Guidelines stated by the INSPIRE Directive (Directive 2007/2/EC of March 14, 2007).

A NTDB enables the structuring of spatial information into layers and classes linked topologically, this interaction should ensure a detailed and digital view of the territory [22]. At the beginning, NTDBs were designed for 2D use, also because there were no suitable tools (hardware and software) for modeling and managing 3D data, so the object-oriented approach and the extension to 3D has not been fully considered [9]. For up-to-date 3D models of the territory, data from different sources must be made congruent with each other and must be integrated. A very significant derivative application of NTDBs relevant to all countries that, like Italy, are affected by the occurrence of large areas subjected to high hydrogeological risk is their use for environmental protection, slope susceptibility analysis, risk assessment, and the resilience of road infrastructure to natural and anthropic effects [23, 24]. Since these aspects were generally not focused when NTDBs were originally developed, when they are used for these purposes the products may show some critical issues, not only in terms of accuracy inherent at the project scale, but also in relation to their reliability with respect to the reality of the territory.

Having this goal in mind, we have developed a procedure to implement the integration into a GIS environment of 2D and 3D data coming from two different sources and of different nature to model the 3D anthropized territory. The procedure was then applied on an area of great interest as it is strongly affected by hydrogeological phenomena that have a strong impact on the viability.

2 Test Area and Data

The area being considered, located in the municipality of Salerno (Fig. 1), is especially vulnerable in terms of hydrogeological risk and is characterized by the simultaneous presence of multiple transport infrastructures concentrated in an urban area of high socio-economic value.

Fig. 1.
figure 1

Test area. (a) Location of the area; (b) Google Earth imagery of the area, the yellow polygon delimits the area being analyzed. (Color figure online)

Full size image

The area analyzed belongs to the strip of land on the slopes of Monte San Liberatore, between Vietri sul Mare and Salerno, which represents a Strategic Infrastructure Corridor for the co-presence of stretches of major international infrastructure: the E45 freeway, the ex-state road SS18, the Naples-Salerno railway line, other minor roads (provincial and municipal) and the commercial port of Salerno. The E45 freeway from Salerno to Cava de’ Tirreni is classified for various stretches at high and very high “rischio frana” and “pericolosità da frana” level in the basin planning in force. The peculiar morphological configuration, characterized by the rocky spur of San Liberatore and the presence of various altimetrically staggered infrastructural networks, makes it ideal for validating the methodology proposed.

On the test area, the available data are the LiDAR from MATTM and the NTDB derived from the “Numerical Regional Technical Cartography” at the scale 1:5,000. The large scale implies the need for periodic updating to monitor evolutionary dynamics within the area, otherwise the NTDB becomes less and less reliable.

The LiDAR data on the area were acquired on April 2013 using an Optech ALTM 3100 laser scanner by Optech, mounted on an aircraft flying at an altitude between 1,500 and 1,800 m height AGL (Above Ground Level). The direction of flight is parallel to the coastline, the maximum scan angle used is 25 deg and the scan frequency has been set equal to 100 kHz.

Detailed information about the acquisition and the dataset can be found in [25]. The LiDAR data were processed and filtered using TerraScan software from Terrasolid. The specifications of the data, as stated by the distributing authority, are as follows: point density greater than 1.5 points per square meter, planimetric accuracy (2σ) of 30 cm, altimetric accuracy (1σ) of 15 cm.

LiDAR data are processed at various levels:

  • XYZ: Point cloud (*.xyz). Each line in the text file contains ellipsoidal coordinates longitude and latitude (DD) and elevation (m), the value of intensity, a code that defines whether the point belongs to the ground (2, “ground”) or not (1, “no ground”). Reference system: WGS84 (EPSG: 4979), ellipsoidal WGS84 height.

  • DSM First, DSM Last: Digital Surface Model, raster file format (*.asc), first-pulse and last-pulse respectively; GSD of 1 m inland and 2 m on the coast. 3D grid cells, planimetric size 10–5 DD. Reference system: EPSG: 4979.

  • DTM: Digital Terrain Model, raster file format (*.asc); GSD of 1 m inland and 2 m on the coast. 3D grid cells, planimetric size 10–5 DD. Reference system: EPSG: 4979. The DTM is derived by TIN interpolation of points classified as “ground”.

  • Intensity: raster file format (*.asc) containing the reflection values of the laser beam. 2D grid cells (10–5 DD). Reference system: EPSG: 4326.

Analyses carried out on some NTDB strata have highlighted anomalies relating to the elevation component, on the polygons of the “Immobili e Antropizzazioni” layer, for the “Opere delle infrastrutture di trasporto” and “Strade” themes. The errors in elevation are presumably due to editing processes that are not appropriate for the 3D geometry: during the editing phases, to some vertices of the polygons of the roads, which overlap with other layers, for example the contour line layer, were attributed incorrect elevation values, as Figs. 2a–b shows.

For LiDAR, starting from the point cloud, given in geographic coordinates φ, λ, first a TIN was formed, a structure particularly suitable for representing complex landforms and highlighting orographic discontinuities, then this TIN was rasterized to have a matrix of cells more easily handled in common GIS software. The coordinate system φ, λ is not an isothermal coordinate system on the ellipsoid, so that an equal angular increment of latitude or longitude corresponds to an arc element of different length.

This affects the raster output, where bands formed by rectangular cells can be observed (Fig. 2c). The other problem related to DTMs is the poor rendering of areas in correspondence with infrastructures or certain geomorphological details, including areas modified by anthropic processes (slopes, roadways, etc.).

Fig. 2.
figure 2

a, b) NTDB: 3D excerpt view of polygons (theme “Strade”) overlaid on contour lines; c) LiDAR data from MATTM: DTM excerpt.

Full size image

3 Methods

Since both input data do not meet the requirements of most experts and users of cartography - the DTM from LiDAR provided by the Ministry has some scaling, and some NTDB polygons have aberrations in the elevation component - to optimize the outputs and obtain a 3D representation of the territory and built environment giving an accurate rendering of geomorphological and anthropic features with their typical shape, an automated tool was implemented in ArcMap environment.

The process is implemented in Python using the ArcGIS library (ArcPy) and is based on the generation of a DTM and contour lines from the LiDAR data using planimetric information from the NTDB. The main steps of the process consist in: (i) homogenization of the different reference systems, (ii) extraction of the height from the LiDAR point cloud and attribution of the same to the 2D polygons implemented in the NTDB, (iii) generation of the DTM from the LiDAR data (iv), generation and filtering of the contour lines.

3.1 Homogenization of Coordinate Reference System

The files provided by the Ministry containing the LiDAR data (*.xyz) are imported into the GIS environment; the function implemented allows to import a single file or a set of files contained in a folder. The output consists of a file named Point Features (*.shp file), in the attributes Table the plano-altimetric coordinates are reported. Only the points associated with the “Ground” class (identified by scalar 2) are imported.

The next step is to associate the data with the correct datum. This is done using the Define Projection tool; for the planimetric coordinates it is EPSG: 4326. Since the reference systems of the input data are different from each other, coordinate/datum transformations that meet the nationally mandated specifications are required, also to enable interoperability between regions.

In Europe, so also in Italy, the system chosen is the ETRS89, realization ETRF2000 (epoch 2008.0), which is mandatory for the Public Administration, stated in the Ministerial Decree of November 10, 2011, as well as being specified in the European Directive INSPIRE (Technical Guidelines Annex I - D2.8.I.1). In Italy ETRF2000 is materialized by the RDN2008 (EPSG:6706).

To transform the coordinates of the LiDAR point cloud into ETRF00, we have used the grids provided by “Istituto Geografico Militare” (IGM), spaced at 7.50’ in longitude and 5’ in latitude. They cover the extension of an element of the map of Italy at a scale of 1:50,000. Each grid node will be associated with a corrective value for longitude (Δλ) and latitude (Δφ) contained in the *.GK2 file [26].

The algorithm implemented to transform the planimetric coordinate of the LiDAR point cloud into the ETRS89/ETRF00 (epoch 2008) involves the following steps:

  1. 1.

    Generation of two 6 × 6 grids (one for each coordinate) having 36 nodes with spacing 7.50’ in longitude and 5’ in latitude (DD). The coordinates of the origin of the grid (in the south-west node) are given in ETRS89/ETRF89 (input system);

  2. 2.

    Computation, through bilinear interpolation of the corrective values given in the *.GK2 file, of the corrections ΔλPi and ΔφPi to be assigned to the generic Pi point of the LiDAR cloud;

  3. 3.

    Conversion of ΔλPi and ΔφPi corrections from DMS into DD;

  4. 4.

    Application of ΔλPi and ΔφPi corrections to LiDAR point coordinates to accomplish the transformation into ETRF2000 (EPSG: 6706), according to the formula:

    $$ \begin{gathered} \lambda_{ETRF00} = \lambda_{ETRF89} + \Delta \lambda \hfill \\ \varphi_{ETRF00} = \varphi_{ETRF89} + \Delta \varphi \hfill \\ \end{gathered} $$
    (1)
  5. 5.

    Conversion from geographic coordinates (EPSG:6706) into UTM/ETRF2000 map coordinates (EPSG:7792).

The *.GK2 file, in addition to the corrective parameters for the transformation of planimetric coordinates, also contains the values of the geoid undulation N (ITALGEO2005) to transform the heights from ellipsoidal h to orthometric H. The values, spaced at 2’ DD in planimetry, are contained in the rows 267–406.

The conversion of the elevation values follows the following steps:

  1. 1.

    Creation of a grid with spacing 2’ DD, made of 14 columns and 10 rows, the planimetric coordinates of the origin of the grid (in the south-west node) are given in ETRS89/ETRF89;

  2. 2.

    Bilinear interpolation of the undulation values, given at the grid nodes, to compute the value for each point in the LiDAR cloud;

  3. 3.

    Computation of orthometric heights H with the known formula:

    $$ H = h - N $$
    (2)

3.2 3D Modeling of the Roads Within the NTDB

The next step is to create the 3D polygons of the roads, with a twofold purpose:

  1. 1.

    To correct the errors on the elevation values of the polygons for the themes “Opere delle infrastrutture di trasporto”, “Ferrovie” and “Strade” in the NTDB;

  2. 2.

    To generate the 3D breaklines to correctly represent the roads, which are elements of discontinuity, in the DTM.

To this end, a specific procedure has been developed to extract altimetric values from the LiDAR point cloud and then attribute them to the 2D polygons of the roads within the NTDB. This will also provide a 3D modeling of the road system.

The process has been implemented in ArcMap environment, with some algorithms developed apart in MATLAB and Python, and it is based on the following steps:

  1. 1.

    Extraction from LiDAR data of points belonging to the road surface;

  2. 2.

    Extraction from NTDB of 2D polygons belonging to the layer “Immobili e Antropizzazioni” (Fig. 3);

  3. 3.

    Division of the polygons, along the longitudinal direction, into sub-polygons of assigned length; the length must be congruent with the resolution of the DTM to be produced;

  4. 4.

    Join process: a column is added to the polygon attributes table; the column will contain the average elevation values of the points within the sub-polygons;

  5. 5.

    Transformation of the 2D sub-polygons into 3D sub-polygons by attributing the elevation values contained within the attribute table, by running the tool Feature to 3D By Attribute.

The process of step #1 has been carried out using an algorithm implemented specifically for Mobile Laser Scanner (MLS) data [27], based on the extraction of ground points, belonging to the road surface, from the MLS trajectory. Unlike MLS data, in which it is possible to extract the trajectory given the scanning angle and the GPS time, in ALS data this input is not available, as the geometry and the point cloud acquisition method are different. Hence, the required information has been extracted from the NTDB; the layer related to the viability contains the polylines, at the axis of the carriageway, of the roads (Fig. 3, about urban and suburban roads). Here the polylines have been used as an alternative to the MLS trajectory for the application of the algorithm. The process requires as input the roadway polyline and the width of the sub-polygons.

Fig. 3.
figure 3

Polygons extracted from the NTDB on the test site: vertices of the polyline of the road network (blue points) overlapped by the polygons of the infrastructure one (purple polygons). (Color figure online)

Full size image

3.3 3D Modeling of the Buildings Within the NTDB

The reconstruction of the 3D model is carried out by joining elementary volumes, means solid that are generated by the extrusion along the vertical of a surface, called the extrusion surface, up to a given height, called the extrusion height. Since the extrusion heights are absolute values, the direction of extrusion can be up or down.

The polygons of the buildings contained in the NTDB are lacking in altimetric values; therefore, they are used as a footprint for the extrusion process, that is based on the following steps:

  1. 1.

    Extraction of the 2D polygons relative to the “Edificato” theme of the “Immobili e Antropizzazioni” layer;

  2. 2.

    Selection of the LiDAR points within each polygon using the Selection by Location function;

  3. 3.

    Removal of the points belonging to the vertical walls of the buildings through the computation of the slope;

  4. 4.

    Selection of the points belonging to the single planes of the pitch by means of the aspect computation;

  5. 5.

    Semi-automatic generation of the pitch planes, through the implementation of the RANdom SAmple Consensus (RANSAC) algorithm, using the primitive plane;

  6. 6.

    Computation of the negative extrusion height (from the pitch towards the ground) for each interpolated plane by computing the difference in elevation between the DTM (average height of the site area) and the minimum height of the pitch plane.

The removal of the points belonging to the vertical walls is done by computing the slope on the TIN surface created from the points selected through the polygons of the buildings [28], the function used is Surface Slope. The function results in a polygonal feature, having the same structure as the mesh. To each polygon is associated the value of the slope, the slope is computed with respect to the horizontal plane. Only the points within the polygons which are associated with slope values lower than the chosen threshold will be extracted. The identification of the points belonging to the single pitch plane, in the case of hipped roofs, is done by computing the aspect on the TIN surface interpolated with the points resulting from the previous process, the function used is Surface Aspect. The function returns a polygonal feature, as for the Surface Slope function, with associated aspect values.

To have the number of classes equal to the number of slopes, a classification according to the Natural Breaks method has been carried out. The points belonging to each single class will be interpolated with the RANSAC algorithm proposed by [29], using the plane primitive.

3.4 DTM Construction

To better represent the land surface, we chose a TIN-type DTM rather than GRID one. To be effective, it must be assumed that the input data are free of outliers.

The implemented script creates the TIN with the Create TIN tool, which respects the Delaunay algorithm. The input data are the LiDAR point cloud (*.xyz) transformed into the reference system EPSG: 7792, with orthometric heights, limited to the ground class.

The tool allows the insertion of breaklines in the form of 3D polygonal surfaces; in our application we used as input the 3D polygons of the roads, obtained from the application of the processes detailed in Sect. 3.2.

After structuring the points into a TIN, the edge triangles were removed using the Delineate TIN tool. The last step of the process consists in the conversion of the 3D TIN into a two-dimensional image (raster), where each pixel is associated with an elevation value, performed by implementing the tool TIN to Raster. The cell size must be chosen according to the representation scale, which in turn must be congruent with the accuracy and density of the input data. The interpolation method used is linear; each triangle (face) belonging to the TIN identifies a plane in space, passing through the three vertices of the triangle itself. Using this interpolating plane, the elevation values of the corresponding raster cell are computed. Generally, as the resolution increases, the output raster will more accurately represent the TIN surface. Contour lines should represent small geomorphological details with their typical shape, making them immediately perceptible. This sometimes requires the enhancement of some terrain shapes, which are generally filtered and smoothed in standard processes [2].

The implemented algorithm is not limited to the generation of contour lines but includes a geometric filtering and smoothing process. The contour lines are generated from the raster DTM by applying the Contour function. New fields representing the geometric characteristics of each contour are added to the contour layer attribute table. For each individual contour line, the following are computed: (i) length; (ii) vertices coordinates.

Through the combination of the computed parameters and the height for each individual contour, the implemented algorithm removes/corrects all possible aberrations or elements that should not be present at a given representation scale. These include bull's eyes, altimetric inconsistencies, non-coincident edges, non-closed areas. The smoothing algorithm used is based on a local polynomial interpolation. The algorithm smooths the lines according to an assigned tolerance; the tolerance controls the maximum distance between the interpolated polyline and the spline vertex. The smaller the length, more detail will be preserved and the longer the processing time. The smoothing process was implemented using the Smooth line tool. The smoothing algorithm used is PAEK (Polynomial Approximation with Exponential Kernel).

4 Results

The accuracy of the coordinate transformation was checked by comparing our results with those resulting from the use of the official IGM conversion software (Verto3k). The average difference in the coordinates between 1,000 random samples was 0.3 mm ± 0.8 mm in East and ± 0.6 mm in North (Fig. 4) whereas for geographical coordinates the average differences was of the order of 1–10−8 DD. The differences are not significant, they are well below the accuracy of the LiDAR data, so the process can be deemed verified. As for the elevations, the comparison provided zero differences.

The construction of the breaklines corresponding to the roads was carried out by the extraction of the polygons relative to the “Immobili e Antropizzazioni” layer, only those relative to the themes of “Strade” and “Ferrovie” belonging to the domain level were selected (Fig. 5a).

The extraction of LiDAR points belonging to the road surface was made by applying the algorithm briefly described in Sect. 3.2 [27]. The maxDistance parameter was set equal to 0.25 m to account for the road cross slope while for the other parameters the defaults were used.

Fig. 4.
figure 4

BoxPlot of the differences between UTM coordinates from Verto3K and those resulting from the implemented process.

Full size image

The polygons belonging to the domain layer were divided into sub-polygons of 0.5 m along the travel direction (Fig. 5b), less would not allow to contain enough points, due to the limited density of the LiDAR data. Each sub-polygon is associated with an elevation that is the average of those of the points, which are selected by a spatial join process. Therefore, each sub-polygon is horizontal and has that elevation value (Fig. 5c); such simplification is justified since the elevation differences between the center and the roadside, due to cross slopes, are lower than the accuracy of the elevation component of LiDAR data.

For the process of building extrusion, the 2D polygons related theme “Edificato” of the “Immobili e Antropizzazioni” layer were extracted. For each polygon, LiDAR points were selected (Fig. 5d), and they were then organized into TINs and the dihedral angle of each triangle face with respect to the horizontal plane was computed using the Surface Slope function. Triangle faces with an angle greater than 45 deg and their vertices were then deleted to remove all those points belonging to the vertical walls of buildings. The Surface Aspect function was applied to the remaining surfaces, and the values obtained were classified according to the Natural Breaks method. This type of classification allowed obtaining as many classes as the number of pitches. The pitch planes were constructed by RANSAC interpolation on the points belonging to the individual classes, using the single primitive plane (Fig. 5e).

The LiDAR point cloud transformed into EPSG: 7792, with orthometric heights, was run through spatial interpolation processes to generate the DTM. The TIN method was used, according to the Delaunay criterion and with the inclusion of breaklines derived from 3D infrastructure polygons. The variable size of the triangles results in a model with variable resolution. In addition, all edge triangles with a length greater than 30 m were deleted. Figure 6 shows the shaded relief maps of the DTM of the MATTM (Fig. 6a) and the one generated through the implemented process (Fig. 6b). In detail, in Fig. 6a it is possible to notice the poor rendering of the roadways, especially on the highway with staggered carriageways and on the shoulders of the viaducts.

In addition, it is also evident the scaling already discussed in Sect. 2. Figure 6b shows the DTM resulting from the optimization process in which it is possible to note a clear improvement in the definition of some geomorphological details, including areas modified by anthropic processes (slopes, roadways, etc.). The TIN is rasterized using a 20 cm pixel size. The effect of the scaling present in the MATTM data is absent, this is due to the interpolation being done in the cartographic system.

Fig. 5.
figure 5

Procedure for constructing breaklines. (a) Extraction of the infrastructures belonging to the domain level; (b) Perspective view of the polygons extracted with the infrastructure point cloud (blue points) overlaid; (c) Zoom-in of the 3D sub-polygons with the point cloud overlaid. d) LiDAR points extracted from 2D polygons of the NTDB; e) Pitch plans interpolated by RANSAC algorithm. (Color figure online)

Full size image
Fig. 6.
figure 6

Shaded relief maps of the test area. (a) DTM provided by the Ministry; (b) optimized DTM.

Full size image

The contour lines are derived from the TIN raster with 1 m spacing. Figure 7 shows a comparison between contour lines derived from the DTM provided by the MATTM (Fig. 7a) and those derived from the optimization process implemented in GIS (Fig. 7b). The contour lines shown in Fig. 7a are featured by discontinuities, essentially related to the DTM interpolation and visible in the shaded relief map of Fig. 6a.

Aberrations have been removed with the automated filtering process, run by setting the tolerance value at 2 m because the tests carried out have proved the effectiveness of setting a tolerance equal to twice the equidistance so to enhance the landform without overly smoothing the edges.

Figure 8 shows the 3D model in ArcScene environment, produced from extrusion of building and infrastructure polygons extracted from NTDB. The DTM was optimized using Breaklines.

Fig. 7.
figure 7

Excerpt of contour maps with 1m spacing equidistance. (a) Contour lines from MATTM-provided DTM; (b) Contour lines from optimized DTM.

Full size image
Fig. 8.
figure 8

a, b, c, d) perspective views of the 3D model of “Vallone e Viadotto Olivieri” created from NTDB-LiDAR integration.

Full size image

5 Conclusions

The work highlights the benefits of an integrated use of regional NTDB and LiDAR data, through automatic procedures implemented using tools in ArcGIS Desktop environment. The integrated use of the two data sets (LiDAR and NTDB) requires data framed in the National Geodetic System and thus needs specific procedures for datum transformation and elevation homogenization.

The automated process allowed to optimize the DTM derived from LiDAR data, through correction and attribution of the altimetric component to the geometries contained in the NTDB, namely to the polygons of the roads, and to generate the 3D breaklines used in the DTM interpolation. The proposed methodology has also allowed to produce a DTM free of the aberrations present in the MATTM one.

The algorithm has been implemented in Python, using the libraries available in ArcMap. Some of the process for the identification of the pitch planes and for the generation of the polygons of the 3D buildings has been developed in a computing environment outside the GIS, given the complexity of the algorithms used and not included in the ArcMap libraries. The implemented methodology offers, moreover, an optimized and detailed 3D representation of the land, infrastructures, and urbanization, useful for the identification and analysis of morpho-evolutionary phenomena that may impact on human activities.