/*   WRR paper is based on vd5r.aml. Prior history in vd5r.aml

/*

/*   mrvbf6a - vd5r modified to incorporate min (fuzzy and) instead of 

/*	multiplication when combining flatness and lowness. The transformation

/*	following combination of flatness and lowness is no longer used, so

/*	tvf and rvf layers are not computed, vf and rf are used directly.

/*

/*   mrvbf6b - mrvbf6a modified to use hard thresholds when recombining

/*	resolutions, rather than the smooth approach that resulted in

/*	values bleeding into lower resolution classes

/*	

/*	-2 changed threshold for converting pctl to low from 0.4 to 0.5

/*

/*   mrvbf6c - tried using standard deviation - not continued

/*

/*   mrvbf6d - another attempt to remove small hilltops from inclusion in

/*	valley bottom classes by including a fraction of the low indicator

/*	from the previous scale

/*

/*	-2 changed weighting of previous low/high from 0.3 to 0.5

/*

/*   mrvbf6f - (ignoring the 6e experiment)

/*	-1 Following results from Linda Gregory's transcription of the

/*	   method to Python in ArcGIS, removed the use of nibble to

/*	   extend the DEM beyond its original borders. Seems to be 

/*	   unnecessary.

/*	   Included extension of grids by one pixel at each step

/*

/*   mrvbf6g - all of 6f, plus std dev

/*         Use sd thresholds with _3 applying to level 3, _9 to level 4 etc 

/*	-2 Reduce standard deviation threshold by a factor of 2, so first

/*	   threshold is 4 instead of 8. This follows analysis of -1 that

/*	   shows that ridge tops are still being left in that should be

/*	   eliminated.

/*	-3 Reduce pctl radius from 6 cells to 3 at all steps. Radius is the

/*	   same for first two steps.

/*	-7 As for -3, but change carrying forward of low indicator to use 

/*	   raw pctl from previous scale, not low weight. This means that the

/*	   information only comes from one scale finer, not all finer scales.

/*	   This is to attempt to remedy the defects introduced by the -3 version

/*	   while retaining its ability to remove some problems.

/*	-10 As for -8, but pctl radius at 6 cells from step 2 onwards.

/*

/*	-a2 Added code to automatically detect geographic coords

/*	    Calculates sdthres from resolution, same as slopethres

/*	-a3 Incorporate mssd calculations directly in the AML, using 

/*	    correct methods developed in SRTM project; changed name of 

/*	    sd layers to match slope naming

&args dem

/* dem: name of dem to work from

&sv keep 1

/* higher keep values result in more intermediate results being kept

&describe %dem%

&sv mindim = [min %grd$ncols% %grd$nrows%]

&sv maxresmult = 1

&sv nsteps = 2

&do &while [calc %mindim% / %maxresmult%] > [calc 5 * 3]

  &sv maxresmult = [calc %maxresmult% * 3]

  &sv nsteps = [calc %nsteps% + 1]

&end

&type Minimum DEM dimension = %mindim%, max resolution = %maxresmult%

&type %nsteps% steps

&if [null %prj$units%] &then

  &return &error No horizontal units defined for DEM %dem%

&select %prj$units%

  &when METERS

    &do

      &type DEM horizontal units is metres, good

      &sv geographic = .false.

    &end

  &when DD

    &do

      &type DEM horizontal units is decimal degrees, scale

      &sv geographic = .true.

    &end

  &otherwise

      &return &error Unrecognised horizontal units for DEM %dem%: %prj$units%

&end

/* set up scaling factors for geographic data

&if %geographic% &then &do

  &type Calculating geographic scaling factors

  &sv xyfactor [calc 110000 * [cos [calc ( %grd$ymin% + %grd$ymax% ) / 2 / 57.2958] ]  ]

  &sv zfactor [calc 1 / %xyfactor% ]

  &end

&else &do

  &sv xyfactor 1

  &sv zfactor 1

&end

&sv baseres = %grd$dx%

&sv resmetres = [calc %baseres% * %xyfactor%]

&type DEM resolution = %resmetres% metres

&sv minvbf = 1

&sv vbfres = 4  /* maximum resolution for first vbf = 1

&sv slopethres = 64  /* 64% = 33 degrees, about steepest for any deposit

&sv sdthres = 0.375  /* developed with 1.5 at 25 m

&sv pctlradius = [calc %vbfres% * 3]

&do &while %vbfres% < %resmetres%

  &sv minvbf = [calc %minvbf% + 1]

  &sv vbfres = [calc %vbfres% * 3]

  &sv slopethres = [calc %slopethres% / 2]

  &sv sdthres = [calc %sdthres% * 2]

  &sv pctlradius = [calc %pctlradius% * 3]

&end

&type CHOICE setting pctlradius = 3 * dem resolution

&sv pctlradius = [calc %resmetres% * 3.02]

/* Modified 18/10/2006 to overcome numerical issues in pctl's cell inclusion calcs

&type TEMPORARY setting minvbf = 1

&sv minvbf = 1

&type VBF for base resolution = %minvbf%

&type Slope threshold = %slopethres%, percentile radius = %pctlradius%

setwindow %dem%

setcell minof

/* set weighting for pctl value from previous step

&sv prevpctlwt = 0.5

&sv l = 0

&sv v = %minvbf%

&sv l = 0

&sv l1 = 0

&sv res = %baseres%

&do &while %l% < %nsteps%

  &if %l% = 0 &then &sv l = 1

    &else &sv l = [calc %l% + 1]

  &type === Step %l% ====

  &if %l% > 2 &then &do

    &sv res = [calc %res% * 3]

    &sv resmetres = [calc %res% * %xyfactor%]

  &end

  &type Resolution = %res% (%resmetres% metres)

  &type Slope threshold = %slopethres%, pctl radius = %pctlradius%

  &type Std dev threshold = %sdthres%

  &type Checkpoint 1

  &if %l% > 2 &then &do

    &type Checkpoint 1.1

    /* create a smoothed version of the DEM from the previous resolution

    &if [exists dem_%l%_%l1% -grid] &then kill dem_%l%_%l1%

    &describe %dem%

    setwindow [calc %grd$xmin% - %res%] [calc %grd$ymin% - %res%] [calc %grd$xmax% + %res%] [calc %grd$ymax% + %res%]

    &if %l% = 3 &then &do

      dem_%l%_%l1% = focalsum(%dem%, weight, gauss3, data) / focalsum(con(%dem% > -999, 1, 0), weight, gauss3, data)

    &end

    &else &do

      dem_%l%_%l1% = focalsum(dem_%l1%_%l1%, weight, gauss3, data) / focalsum(con(dem_%l1%_%l1% > -999, 1, 0), weight, gauss3, data)

    &end

    &type Checkpoint 1.2

    /* calculate standard deviation

    &if [exists sd_%l% -grid] &then kill sd_%l%

    &if %l% = 3 &then &do

      sd_%l% = resample(blockstd(%dem%, rectangle, 3, 3), %res%, nearest)

    &end

    &else &do

      &if [exists sswg_%l% -grid] &then kill sswg_%l%

      &if [exists ssbg_%l% -grid] &then kill ssbg_%l%

      sswg_%l% = aggregate(pow(sd_%l1%, 2), 3, sum)

      ssbg_%l% = 9 * pow(resample(blockstd(dem_%l1%_%l1%, rectangle, 3, 3), %res%, nearest), 2)

      sd_%l% = sqrt((sswg_%l% + ssbg_%l%) / 9)

      &if %keep% <= 4 &then &do

        kill ssbg_%l%

        kill sswg_%l%

      &end

    &end

    &type Checkpoint 1.3

    &if %keep% <= 2 &then &do

      /* get rid of DEM from previous step

      &if %l1% > 2 &then kill dem_%l1%_%l1%

      &if %l1% > 2 &then kill sd_%l1%

    &end

  &end

  &type Checkpoint 2

  &if [exists s_%l%_1 -grid] &then kill s_%l%_1

  &if %l% = 1 &then &do

    /* first step, compute slope from base resolution

    &type Compute slope from DEM

    s_1_1 = slope(%dem%, %zfactor%, percentrise)

  &end 

  &else &if %l% = 2 &then &do

    /* second step, slope is same as from first step

    &if %keep% <= 4 &then rename s_1_1 s_2_1

    &else copy s_1_1 s_2_1

  &end

  &else &do

    /* compute slope from smoothed DEM

    &if [exists s_%l%_%l1% -grid] &then kill s_%l%_%l1%

    s_%l%_%l1% = slope(dem_%l%_%l1%, %zfactor%, percentrise)

    &if %l% > 2 &then &do

      /* resample back to base resolution

      s_%l%_1 = resample(s_%l%_%l1%, %baseres%, bilinear)

    &end

    &if %keep% <= 5 &then kill s_%l%_%l1%

    &if [exists dem_%l%_%l% -grid] &then kill dem_%l%_%l%

    dem_%l%_%l% = resample(dem_%l%_%l1%, %res%)

    &if %keep% <= 5 &then kill dem_%l%_%l1%

  &end

  &if [exists f_%l% -grid] &then kill f_%l%

  f_%l% = 1 / (1 + pow(s_%l%_1 / %slopethres%, 4))

  &if %keep% <= 4 &then &do

    /* get rid of s_%l%_1, except when l = 1 and we need it for s_2_1

    &if %l% ^= 1 &then kill s_%l%_1

  &end

  &type Checkpoint 3

  &if [exists pctl_%l%_1 -grid] &then kill pctl_%l%_1

  &if [exists cpctl_%l%_1 -grid] &then kill cpctl_%l%_1

  &if [exists cf_%l% -grid] &then kill cf_%l%

  &if [exists low_%l% -grid] &then kill low_%l%

  &if [exists vf_%l% -grid] &then kill vf_%l%

  &if [exists high_%l% -grid] &then kill high_%l%

  &if [exists rf_%l% -grid] &then kill rf_%l%

  &if %l% <= 2 &then &do

    &r chunkproc pctl %dem% pctl_%l%_1 %pctlradius% %geographic%

    &if %l% = 1 &then &do

      low_%l% = 1 / (1 + pow(pctl_%l%_1 / 0.5, 3))

    &end

    &else &do

      cpctl_%l%_1 = [calc 1 - %prevpctlwt%] * pctl_%l%_1 + %prevpctlwt% * pctl_%l1%_1

      low_%l% = 1 / (1 + pow(cpctl_%l%_1 / 0.5, 3))

    &end

    vf_%l% = min(f_%l%, low_%l%)

    &if %l% = 1 &then &do

      high_%l% = 1 / (1 + pow((1 - pctl_%l%_1) / 0.5, 3))

    &end

    &else &do

      high_%l% = 1 / (1 + pow((1 - cpctl_%l%_1) / 0.5, 3))

    &end

    rf_%l% = min(f_%l%, high_%l%)

    &if %l% = 2 &then &do

      cf_2 = f_2 * f_1

      &if %keep% <= 3 &then kill f_1

    &end

  &end

  &else &do

    &if [exists pctl_%l%_%l% -grid] &then kill pctl_%l%_%l%

    &r chunkproc pctl dem_%l%_%l% pctl_%l%_%l% %pctlradius% %geographic%

    pctl_%l%_1 = resample(pctl_%l%_%l%, %baseres%, bilinear)

    cpctl_%l%_1 = [calc 1 - %prevpctlwt%] * pctl_%l%_1 + %prevpctlwt% * pctl_%l1%_1

    &if %keep% <= 5 &then kill pctl_%l%_%l%

    cf_%l% = f_%l% * cf_%l1%

    &if %keep% <= 3 &then kill cf_%l1%

    low_%l% = 1 / (1 + pow(cpctl_%l%_1 / 0.5, 3))

    &if [exists lowrelief_%l% -grid] &then kill lowrelief_%l%

    lowrelief_%l% = 1 / (1 + pow(sd_%l% / %sdthres%, 4))

    vf_%l% = min(cf_%l%, low_%l%, lowrelief_%l%)

    high_%l% = 1 / (1 + pow((1 - cpctl_%l%_1) / 0.5, 3))

    rf_%l% = min(cf_%l%, high_%l%, lowrelief_%l%)

    &if %keep% <= 3 &then kill lowrelief_%l%

  &end

  &type Checkpoint 4

  &if %l% > 1 &then &do

    &if %keep% <= 4 &then kill pctl_%l1%_1

    &if %keep% <= 4 &then kill cpctl_%l%_1

    &if %keep% <= 3 &then kill f_%l%

  &end

  &if [exists mrvbf_%l% -grid] &then kill mrvbf_%l%

  &if [exists mrrtf_%l% -grid] &then kill mrrtf_%l%

  &if %l% = 2 &then &do

    mrvbf_%l% = con(vf_2 >= 0.5, vf_2 + %v% - 1, vf_1)

    mrrtf_%l% = con(rf_2 >= 0.5, rf_2 + %v% - 1, rf_1)

    &if %keep% <= 3 &then kill vf_1

    &if %keep% <= 3 &then kill rf_1

  &end

  &else &if %l% > 2 &then &do

      mrvbf_%l% = con(vf_%l% >= 0.5, vf_%l% + %v% - 1, mrvbf_%l1%)

      mrrtf_%l% = con(rf_%l% >= 0.5, rf_%l% + %v% - 1, mrrtf_%l1%)

     &if %keep% <= 3 &then kill mrvbf_%l1%

     &if %keep% <= 3 &then kill mrrtf_%l1%

  &end

  &if %l% >= 2 &then &do

    &if %keep% <= 3 &then kill vf_%l%

    &if %keep% <= 3 &then kill rf_%l%

  &end

  &type Checkpoint 5

  &sv v = [calc %v% + 1]

  &sv l1 = %l%

  &sv slopethres = [calc %slopethres% / 2]

  &if %l% > 2 &then &do

    &sv sdthres = [calc %sdthres% * 2]

  &end

  &if %l% = 1 &then &sv pctlradius = [calc %pctlradius% * 2]

  &else &sv pctlradius = [calc %pctlradius% * 3]

  &type Checkpoint 6

&end

&if %keep% <= 4 &then kill cf_%l%

&if %keep% <= 4 &then kill pctl_%l%_1

&if %keep% <= 2 &then kill dem_%l%_%l%

&if [exists mrvbf6g-a -grid] &then kill mrvbf6g-a

&if [exists mrrtf6g-a -grid] &then kill mrrtf6g-a

mrvbf6g-a = con(%dem% > -999, con(mrvbf_%l% > [calc %minvbf% - 0.5], mrvbf_%l%, 0))

mrrtf6g-a = con(%dem% > -999, con(mrrtf_%l% > [calc %minvbf% - 0.5], mrrtf_%l%, 0))

&if %keep% <= 2 &then kill mrvbf_%l%

&if %keep% <= 2 &then kill mrrtf_%l%