• If we just consider a straight segment V1-V2 we'll wrongly assume that the elevation profile has a constant slope.
  (as shown by case a) on the side figure)
• But the actual elevation profile of the underlaying surface could easily be completely different and much more articulated.
  (as shown by cases b) and c) on the side figure)
• Conclusion: we definitely need a more sophisticated approach when draping LINESTRINGs and POLYGONs.
  One safely ensuring to be absolutely sure to never suppress any relevant information about Elevations available from the raster grid.

• We just need to preliminary densify our LINESTRING or POLYGON Geometries by interpolating more Vertices at regular intervals.
  In other words, we must always call ST_Segmentize() before draping our geometries over the DEM/DTM grid.
• Such interpolated Vertices should ideally ensure that every relevant cell on the raster grid will be effectively overlapped.
  And accordingly to the Shannon sampling theorem the optimal length limit to be passed to ST_Segmentize() is half of the cell size.

Warning

• Huge Geometries presenting a very high number of Vertices are usually very slow to be evaluated by SQL Functions such as ST_Intersects etc.
• So we must necessarily simplify all draped geometries in the post-processing steps.
  This is a classic two-steps approach:
  • First of all we'll always start by appropriately densifying all Geometries before draping.
  • Once draping is done, we must simplify in order to get our Geometries as lean as possible:
    • Only the most relevant interpolated vertices will be preserved.
      i.e. the ones corresponding to some relevant variation in Elevation.
    • Any other interpolated vertex will be simply discarded, because it's not really useful.
  • It's a reasonable and effective compromise: our draped Geometries will certainly put on some extra weight.
    But no more of what's strictly indispensable for avoiding any serious loss of information and degraded accuracy.

step #1.a

SELECT ImportSHP( 'C:\vanuatu\toscana\Iternet_Marzo_2019\shp\civici' , 
       'civici' , 'CP1252' , 3003 , 'Geometry' , 'PK_UID' , 'POINT' , 0 , 0 , 1 , 1 , 'LOWERCASE' );

step #1.b

SELECT ImportSHP( 'C:\vanuatu\toscana\Iternet_Marzo_2019\shp\elem_strad' ,
       'elem_strad' , 'CP1252' , 3003 , 'Geometry' , 'PK_UID' , 'LINESTRING' , 0 , 0 , 1 , 1 , 'LOWERCASE' );

step #1.c

SELECT RL2_CreateRasterCoverage ( 'dtmoro' , 'FLOAT' , 'DATAGRID' , 1 , 'LZMA' , 100 , 512 , 512 , 3003 , 10.0 , 10.0 , 
       RL2_SetPixelValue ( RL2_CreatePixel ( 'FLOAT' , 'DATAGRID' , 1 ) , 0, -9999 ) , 1 , 0 , 1 , 1 , 1 , 1 );

step #1.d

SELECT RL2_LoadRaster ( 'dtmoro' , '/home/sandro/Scaricati/draping/dtmoro.asc' , 0 , 3003 , 1 , 1 );

Important notice for Windows users

• The most recent and accurate Lidar Coverages should be processed with the highest priority
• Then the less recent and accurate Lidar Coverage should be processed in decreasing priority order.
• Finally the medium-.resolution 10x10 dtmoro should be processed as a last choice for filling any remaining gap uncovered by Lidar DTM/DSM.

CREATE TEMPORARY TABLE draping_aux (
    progr INTEGER NOT NULL,
    db_prefix TEXT,
    coverage_name TEXT NOT NULL);

• the name of the helper table is not relevant, but the names of the columns must exactly be as shown in the above example.
• progr is the relative priority index (the lower the value, the higher the priority).
• db_prefix is the prefix of the attached DB containing the Raster Coverage.
  When it's NULL the MAIN DB will be always intended.
• coverage_name is the name of a Raster Coverage to be processed for Draping Geometries.

--
-- inserting all Lidar Raster Coverages in decreasing priority order
--
...
INSERT INTO draping_aux ( coverage_name , progr ) VALUES ( 'lidar_green_2011' , 10 );
INSERT INTO draping_aux ( coverage_name , progr ) VALUES ( 'lidar_green_2010' , 11 );
INSERT INTO draping_aux ( coverage_name , progr ) VALUES ( 'lidar_red_2009' , 12 );
INSERT INTO draping_aux ( coverage_name , progr ) VALUES ( 'lidar_red_2008' , 13 );
...

--
-- inserting the 10x10 dtmoro Raster Coverage with lowest priority
--
INSERT INTO draping_aux ( coverage_name , progr ) VALUES ( 'dtmoro' , 99 );

SELECT RL2_DrapeGeometries( NULL , NULL , 'draping_aux' , 'civici' , 'geometry' , 'geom3d' , -9999 , 5.0 , 5.0 , 0 );
---------------------------------------
1



SELECT RL2_DrapeGeometries( NULL , NULL , 'draping_aux' , 'elem_strad' , 'geometry' , 'geom3d' , -9999 , 5.0 , 5.0 , 0 );
---------------------------------------
1

• the second argument (Raster Coverage name) is now NULL, because we are now processing multiple Raster Coverages.
• the third argument (that was previously NULL) now specifies the name of the helper table containing the list of all Raster Coverages to be processed in the given priority order.
  Note: these two arguments are mutually exclusive: you can define the one or the other, but not both at the same time. And you must specify one of them.

Useful hint

SELECT db_prefix, coverage_name
FROM draping_aux
ORDER BY progr;

1. RL2_DrapeGeometries() will always start by creating the 3D geometry column.
2. then it will initially populate the 3D geometry column by inserting into it all geometries found within the 2D geometry column:
   • all X and Y (and eventually M) coordinates will be preserved as they are.
   • all Z coordinates will be initially set to the NO-DATA value.
   • if a densification distance has been specified as many vertices as required will be interpolated into each 3D geometry.
3. now RL2_DrapeGeometry() will start processing the first Raster Coverage (or the unique one if a just single Raster Coverage has been specified).
4. the Raster Coverage will be processed by extracting a single Tile at each time.
   All Geometries intersecting the current Tile will be processed (and this explains why the Spatial Table must be supported by a Spatial Index).
   All points/vertices intersecting the Tile and still presenting NO-DATA values will assume the elevation value taken from the corresponding cells of the Tile.
   Any point/vertex presenting a Z-value different from NO-DATA will never change in a second time, and this explains why the execution priority of each Raster Coverage is critically relevant.
5. when all Tiles of the current Raster Coverage will be finally processed, then the next Raster Coverage (if any) will be processed.
   The cycle will resume from step 3) until consuming all Raster Coverages.
6. the final step will be the one to simplify all 3D geometries accordingly to the given simplification distance.

SELECT RL2_DrapeGeometries( NULL , 'dtmoro' , NULL , 'civici' , 'geometry' , 'geom3d' , -9999 , 5.0 , 5.0 , 1 );
---------------------------------------
1

• the tenth and last argument accepted by RL2_DrapeGeometries() is update_m, and it's of the boolean type:
  • if set to FALSE (this is the default setting) only the Z-values will be affected by Draping, and all M-values will be preserved untouched.
  • if set to TRUE only the M-values will be affected by Draping, and all Z-values will be preserved untouched.