Hello all,
For those interested in generating Local Relief Models based on the methodology of Hesse (2010), I have made a simple batch file to automate the task:
Ralf Hesse (2010), LiDAR derived Local Relief Models – a new tool for archaeological prospection, Archaeological Prospection 17(2):67 - 72. DOI:10.1002/arp.374
Script (Windows):
REM 1. Low pass filter original DTM/DEM
REM 2. Original DEM (SX57NE.grd) - low-pass DEM -> Difference grid
REM 3. Extract zero contours from difference grid
REM 4. Extract point coordinates/z-values from zero contour data applied to original DEM
REM 5. Grid surface from zero-contour extracted points
REM 6. Local relief model = Original DEM - Purged simplified model DEM
grdinfo SX57NE.grd -C | gawk “{xmin=$2; xmax=$3; ymin=$4; ymax=$5; nx=$10; ny=$11}; {print xmin, xmax, ymin, ymax, 4 * nx, 4 * ny}” > xycoord.txt
@ for /f “tokens=1” %%i in (xycoord.txt) do set xmin=%%i
@ for /f “tokens=2” %%i in (xycoord.txt) do set xmax=%%i
@ for /f “tokens=3” %%i in (xycoord.txt) do set ymin=%%i
@ for /f “tokens=4” %%i in (xycoord.txt) do set ymax=%%i
@ for /f “tokens=5” %%i in (xycoord.txt) do set nx=%%i
@ for /f “tokens=6” %%i in (xycoord.txt) do set ny=%%i
REM 1. Low-pass filter (test values for best result) = DEM_LP where crop.grd=original DEM
grdfft SX57NE.grd -GLP_DEM.grd -F-/50 -N%nx%/%ny% -V
REM 2. DIFF_DEM = DEM - DEM_LP Simple Local Relief Model (approximation)
grdmath SX57NE.grd LP_DEM.grd SUB -V = DIFF_DEM.grd
REM 3. ZERO_CONTOURS from DIFF_DEM extract DIFF_POINTS
grdcontour DIFF_DEM.grd -C0, -DDIFF_DEM_Zero-cont.dat -V
REM 4. Extract DEM_POINTS_ORIGINAL using DIFF_POINTS as input
grdtrack DIFF_DEM_Zero-cont.dat -GSX57NE.grd -V > Simplified_DEM.dat
REM Pre-process data for Surface
blockmean Simplified_DEM.dat -R%xmin%/%xmax%/%ymin%/%ymax% -i0,1,3 -I5 -V > Simplified_DEM_mean.dat
REM 5. Grid DEM_POINTS_ORIGINAL to SIMPLIFIED_DEM XYZ (Minimum curvature)
surface Simplified_DEM_mean.dat -R%xmin%/%xmax%/%ymin%/%ymax% -I1 -GSimplified_DEM.grd -T0 -C0.1 -V
REM 6. LRM_DEM = ORIGINAL_DEM - SIMPLIFIED_DEM
grdmath SX57NE.grd Simplified_DEM.grd SUB -V = Local_Relief_Model_revised.grd
Note: the low-pass filter would need some experimenting with to achieve the best result for the LiDar data resolution.
Hope this is useful.
Lester