export REPO="."

# first migrate the whole thing to new commons area

mv $ORGANISMS/topo $LAYERS/commons/

# consolidate first set of dirs

mkdir $LAYERS/commons/topo/v1

mv -i $LAYERS/commons/topo/{experimental,scripts,tiles} \

      $LAYERS/commons/topo/v1/

# consolidate second set of dirs

mv -i $LAYERS/commons/topo/experimental_2/topo \

      $LAYERS/commons/topo/v2

rmdir $LAYERS/commons/topo/experimental_2

# remove ArcToolbox temporary files

rm $LAYERS/commons/topo/v2/topo/xx00003024.s

rm $LAYERS/commons/topo/v2/experimental/xx00013916.x

rm $LAYERS/commons/topo/v2/experimental/xx00003916.s

rm $LAYERS/commons/topo/v2/experimental/xx00013916.s

# consolidate first set of scripts into git repo

mv -i $LAYERS/commons/topo/v1/experimental/layers.aml \

      $REPO/terrain/research/oregon/arcgis/v1/

mv -i $LAYERS/commons/topo/v1/experimental/aggregate.aml \

      $REPO/terrain/research/oregon/arcgis/v1/

mv -i $LAYERS/commons/topo/v1/experimental/multiscalesmooth9a_clean.aml \

      $REPO/terrain/research/oregon/arcgis/v1/

mv -i $LAYERS/commons/topo/v1/experimental/tilemerge_oregon1.aml \

      $REPO/terrain/research/oregon/arcgis/v1/

mv -i $LAYERS/commons/topo/v1/scripts/unpacktiles.aml \

      $REPO/terrain/research/oregon/arcgis/v1/

mv -i $LAYERS/commons/topo/v1/scripts/unziptiles.py \

      $REPO/terrain/research/oregon/arcgis/v1/

mv -i $LAYERS/commons/topo/v1/scripts/mrvbf/mrvbf6g-a3.aml \

      $REPO/terrain/research/oregon/arcgis/v1/

mv -i $LAYERS/commons/topo/v1/scripts/mrvbf/pctl.aml \

      $REPO/terrain/research/oregon/arcgis/v1/

mv -i $LAYERS/commons/topo/v1/scripts/mrvbf/chunkproc.aml \

      $REPO/terrain/research/oregon/arcgis/v1/

mv -i $LAYERS/commons/topo/v1/scripts/mrvbf/pctl-limits.aml \

      $REPO/terrain/research/oregon/arcgis/v1/

mv -i $LAYERS/commons/topo/v1/scripts/mrvbf/gauss3 \

      $REPO/terrain/research/oregon/arcgis/v1/

# for now keep this mysterioud Windows binary around...

mv -i $LAYERS/commons/topo/v1/scripts/mrvbf/pctl.exe \

      $LAYERS/commons/topo/v1/

# remove now-empty scripts directory

rmdir $LAYERS/commons/topo/v1/scripts/mrvbf

rmdir $LAYERS/commons/topo/v1/scripts

# consolidate second set of scripts into git repo

mv -i $LAYERS/commons/topo/v2/experimental/layers.aml \

      $REPO/terrain/research/oregon/arcgis/v2/

mv -i $LAYERS/commons/topo/v2/experimental/aggregate.aml \

      $REPO/terrain/research/oregon/arcgis/v2/

mv -i $LAYERS/commons/topo/v2/experimental/multiscalesmooth9a_clean.aml \

      $REPO/terrain/research/oregon/arcgis/v2/

mv -i $LAYERS/commons/topo/v2/experimental/tilemerge_oregon1.aml \

      $REPO/terrain/research/oregon/arcgis/v2/

mv -i $LAYERS/commons/topo/v2/scripts/unpacktiles.aml \

      $REPO/terrain/research/oregon/arcgis/v2/

mv -i $LAYERS/commons/topo/v2/scripts/unziptiles.py \

      $REPO/terrain/research/oregon/arcgis/v2/

mv -i $LAYERS/commons/topo/v2/topo/mrvbf6g-a3.aml \

      $REPO/terrain/research/oregon/arcgis/v2/

mv -i $LAYERS/commons/topo/v2/topo/pctl.aml \

      $REPO/terrain/research/oregon/arcgis/v2/

mv -i $LAYERS/commons/topo/v2/topo/chunkproc.aml \

      $REPO/terrain/research/oregon/arcgis/v2/

mv -i $LAYERS/commons/topo/v2/topo/pctl-limits.aml \

      $REPO/terrain/research/oregon/arcgis/v2/

mv -i $LAYERS/commons/topo/v2/topo/gauss3 \

      $REPO/terrain/research/oregon/arcgis/v2/

/* need to ensure that the aggregation boxes are consistent across the whole set of layers

/* so might need a setwindow here to enforce the correct location

/* Question about 1km Behrman version - should that be projected from 30" stats, or

/* should data be projected to Behrman and aggregated over 1km cell. Far from the equator

/* these will give quite different results

/* elevation - derive mean, min, max, std. dev, skew

elev_mean = aggregate(srtmv41_smth, 10, mean)

elev_min = aggregate(srtmv41_smth, 10, min)

elev_max = aggregate(srtmv41_smth, 10, max)

elev_range = elev_max - elev_min

/* stupid arc doesn't have a standard deviation function for aggregation

elev_sd = sqrt(aggregate(sqr(srtmv41_smth - elev_mean), 10, sum) / 100)

/* slope - just mean and standard deviation

slopedeg_mean = aggregate(slopedeg, 10, mean)

slopedeg_sd = sqrt(aggregate(sqr(slopdeg - slopedeg_mean), 10, sum) / 100)

/* wetness index - mean and standard deviation

/* twi not done yet

/* Invoke another AML on a grid to produce a grid

/* Other AML takes arguments: gridname outgridname otherargs

&args targetaml gridname outgridname otherargs:REST

/* targetaml-limits.aml takes same arguments as target.aml (except for

/* output file) and sets maximum number of cells and overlap required

&run %targetaml%-limits %gridname% [unquote %otherargs%]

&describe %gridname%

&sv gridxmin = %grd$xmin%

&sv gridxmax = %grd$xmax%

&sv gridymin = %grd$ymin%

&sv gridymax = %grd$ymax%

&sv gridncols = %grd$ncols%

&sv gridnrows = %grd$nrows%

&sv cellsize = %grd$dx%

/****** TEMPORARY *******/

/* setcell %cellsize%

/**** END TEMPORARY *****/

&sv ncells = %grd$ncols% * %grd$nrows%

&type %ncells% cells in %gridname%

&if %ncells% <= %.maxcells% &then &do

  /* no need to do chunks, just invoke target AML

  &type Process as single chunk

  &run %targetaml% %gridname% %outgridname% [unquote %otherargs%]

  &return

&end

/* need to do chunks

&type Process in multiple chunks

/* &type X: %gridxmin% %gridxmax% (%gridncols%) 

/* &type Y: %gridymin% %gridymax% (%gridnrows%) 

/* &type %ncells% cells, %.maxcells% max, %.border% border

&sv nchunks = [truncate [calc %ncells% / %.maxcells%]] + 1

&type Need at least %nchunks% chunks

&sv chunkdim = [sqrt %.maxcells%] + 2 * %.border%

&if %gridncols% < %chunkdim% &then

  &do

    &sv colchunks = 1

    &sv rowchunks = %nchunks%

    &sv chunkncols = %gridncols%

    &sv chunknrows = [truncate [calc %gridnrows% / %rowchunks%]] + 1

  &end

&else &if %gridnrows% < %chunkdim% &then

  &do

    &sv rowchunks = 1

    &sv colchunks = %nchunks%

    &sv chunkncols = [truncate [calc %gridncols% / %colchunks%]] + 1

  &end

&else 

  &do

    &sv rowchunks = [truncate [calc %gridnrows% / %chunkdim%]] + 1

    &sv colchunks = [truncate [calc %gridncols% / %chunkdim%]] + 1

    &sv chunkncols = [truncate [calc %gridncols% / %colchunks%]] + 1

    &sv chunknrows = [truncate [calc %gridnrows% / %rowchunks%]] + 1

  &end

&type %colchunks% colchunks, %rowchunks% rowchunks

&type Chunk size %chunkncols% cols by %chunknrows% rows

&sv qtrcell = %cellsize% / 4

&sv halfcell = %cellsize% / 2

&sv colist

&sv chunki = 1

&do &while %chunki% <= %colchunks%

  &sv chunkj = 1

  &sv colist_col

  &do &while %chunkj% <= %rowchunks%

    &type Chunk %chunki% %chunkj%

    &sv coxmin = %gridxmin% + ( %chunki% - 1 ) * %chunkncols% * %cellsize% - %qtrcell%

    &sv coxmax = %coxmin% + %chunkncols% * %cellsize% + %halfcell%

    &sv coymin = %gridymin% + ( %chunkj% - 1 ) * %chunknrows% * %cellsize% - %qtrcell%

    &sv coymax = %coymin% + %chunknrows% * %cellsize% + %halfcell%

    /* &type Output X: %coxmin% %coxmax%

    /* &type Output Y: %coymin% %coymax%

    &if %chunki% = 1 &then &sv chunkxmin = %coxmin%

    &else &sv chunkxmin = %coxmin% - %.border% * %cellsize%

    &if %chunki% = %colchunks% &then &sv chunkxmax = %coxmax%

    &else &sv chunkxmax = %coxmax% + %.border% * %cellsize%

    &if %chunkj% = 1 &then &sv chunkymin = %coymin%

    &else &sv chunkymin = %coymin% - %.border% * %cellsize%

    &if %chunkj% = %rowchunks% &then &sv chunkymax = %coymax%

    &else &sv chunkymax = %coymax% + %.border% * %cellsize%

    /* &type Chunk X: %chunkxmin% %chunkxmax%

    /* &type Chunk Y: %chunkymin% %chunkymax%

    /* clip input file to chunk (with borders)

    setwindow %chunkxmin% %chunkymin% %chunkxmax% %chunkymax% %gridname%

    &if [exists chunk -grid] &then kill chunk

    chunk = %gridname%  

    &if [exists co_tmp -grid] &then kill co_tmp

    &run %targetaml% chunk co_tmp [unquote %otherargs%]

    kill chunk

    /* clip output file to chunk (without borders)

    setwindow %coxmin% %coymin% %coxmax% %coymax% %gridname%

    &if [exists co_%chunki%_%chunkj% -grid] &then kill co_%chunki%_%chunkj%

    co_%chunki%_%chunkj% = co_tmp

    kill co_tmp

    &if %chunkj% = 1 &then &sv colist_col = co_%chunki%_%chunkj%

    &else &sv colist_col = %colist_col%, co_%chunki%_%chunkj%

    &sv chunkj = %chunkj% + 1

  &end

  setwindow maxof

  &if [exists co_%chunki% -grid] &then kill co_%chunki%

  co_%chunki% = merge(%colist_col%)

  kill (!%colist_col%!)

  &if %chunki% = 1 &then &sv colist = co_%chunki%

  &else &sv colist = %colist%, co_%chunki%

  &sv chunki = %chunki% + 1

&end

setwindow maxof

&if [exists %outgridname% -grid] &then kill %outgridname%

%outgridname% = merge(%colist%)

kill (!%colist%!)

11 11

0.017 0.046 0.100 0.174 0.242 0.271 0.242 0.174 0.100 0.046 0.017 

0.046 0.124 0.271 0.472 0.659 0.736 0.659 0.472 0.271 0.124 0.046 

0.100 0.271 0.590 1.028 1.434 1.603 1.434 1.028 0.590 0.271 0.100 

0.174 0.472 1.028 1.791 2.500 2.793 2.500 1.791 1.028 0.472 0.174 

0.242 0.659 1.434 2.500 3.488 3.898 3.488 2.500 1.434 0.659 0.242 

0.271 0.736 1.603 2.793 3.898 4.356 3.898 2.793 1.603 0.736 0.271 

0.242 0.659 1.434 2.500 3.488 3.898 3.488 2.500 1.434 0.659 0.242 

0.174 0.472 1.028 1.791 2.500 2.793 2.500 1.791 1.028 0.472 0.174 

0.100 0.271 0.590 1.028 1.434 1.603 1.434 1.028 0.590 0.271 0.100 

0.046 0.124 0.271 0.472 0.659 0.736 0.659 0.472 0.271 0.124 0.046 

0.017 0.046 0.100 0.174 0.242 0.271 0.242 0.174 0.100 0.046 0.017 

/* using srtmv41 3x3 degree patch smoothed to srtmv41_smth

&describe srtmv41_smth

&sv cellsize %grd$dx% * 111000

dzdx = (srtmv41_smth(1,0) - srtmv41_smth(-1,0)) / (2 * %cellsize% * cos($$ymap / deg))

dzdy = (srtmv41_smth(0,-1) - srtmv41_smth(0,1)) / (2 * %cellsize%)

slope = sqrt(sqr(dzdx) + sqr(dzdy))

slopedeg = atan(slope) * deg

/* use scaling factor that is average of x and y scaling factors, 0.707 and 1 ~= 0.85

/* and 0.85 / 111000 ~ 0.0000107

arcslope = slope(srtmv41_smth, 0.0000107, degree)

hillshade = hillshade(srtmv41_smth, 300, 30, shade, 0.00001)

aspect1 = atan2(0 - dzdy, 0 - dzdx) * deg

aspect = con(aspect1 > 90, 450 - aspect1, 90 - aspect1)

kill aspect1

/*use elevation 

flowdir = flowdirection(srtmv41_smth)

/*cell area to account for latitudinal area effects, smallest area near pole

setwindow srtmv41_smth

setcell srtmv41_smth 

cellarea = sqr(%cellsize%) * cos($$ymap / deg)

flowacc = flowaccumulation(flowdir, cellarea)

/*account for sinks

/*change # defaulto z limit

/*caution fill is a command not function

/*fill srtmv41_smth srtmv41_fill sink # flowdir_fill

flowacc_fill = flowaccumulation(flowdir_fill, cellarea)

/*scatch_fill is specific catchment area in m 

scatch_fill = flowacc_fill / sqrt(cellarea) + sqrt(cellarea) / 2

/*setting the min slope to 0.1

adjslope = tan(max(slopedeg, 0.1) / deg)

twi_fill = ln(scatch_fill / adjslope)

/*   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%

&args ingrid sd prob bbox:rest

/* 9a - limit to 4 steps, see if that has any significant deterioration of smoothing performance. Should fix

/* the problem with islands and headlands - turns out also need to remove water except

/* for one cell adjacent to land, and give that a higher uncertainty. See cluster_multiscalesmooth9a_clean_216

/* Version 9:

/* Focus on implementing identical algorithm to directsmooth2 using multiscale method i.e. aggregating by factor of 3 

/* from already aggregated data, rather than from original resolution each time.

/* Version 8:

/* WARNING: have not yet checked that the additional weighting of the gaussian smoothing is not messing with

/* the calculations of variance etc.

/* Replaced simple 3x3 aggregation with gaussian smoothing

/* Kernel is chosen to give appropriate level of smoothing and produce a fairly accurate approximation of the smoothed

/* surface by interpolation of the coarser scale smoothed values

/* Details in gaussian.xls on john's laptop

/* Version 7:

/* Further reworking of how the final values are selected - a mean is a candidate if its associated variance

/* is less than the mean sample uncertainty, and the mean with the lowest variance among the candidates is the chosen one.

/* Implement that in the nested sense by taking the lowest group variance divided by the chi^2 value, and its associated mean variance,

/* and if that is lower than the data point variance the 

/* approximate critical value of chi^2/N with N degrees of freedom at 5% level as 1 + 2.45/sqrt(N) + 0.55/N

/* for 1% level use 1 + 3.4/sqrt(N) + 2.4/N

/* Version 6:

/* Done from scratch after careful working through of theory.

/* ingrid is the (potentially sparse) grid of data to be smoothed and

/* interpolated, which can be of different extent to the output

/* (resolution is assumed to be the same, could adapt this to relax that)

/*

/* var can be a constant or a grid

/*

/* bbox can be either a grid name or the 'xmin ymin xmax ymax' parameters

/* for setwindow

&type NB - using standard deviation as noise specification now, not variance!

/* set up chisq parameters

&sv chisqa = [calc 2.807 - 0.6422 * [log10 %prob% ] - 3.410 * %prob% ** 0.3411 ]

&sv chisqb = [calc -5.871 - 3.675 * [log10 %prob% ] + 4.690 * %prob% ** 0.3377 ]

&type chisq parameters %chisqa% %chisqb%

setcell %ingrid%

/* work out maximum of ingrid and bbox extents

setwindow [unquote %bbox%] %ingrid%

bboxgrid = 1

setwindow maxof

workgrid = %ingrid% + bboxgrid

setwindow workgrid

kill workgrid

/* naming:

/*  h - the value being smoothed/interpolated

/*  vg - total variance of group of data, or of individual measurement

/*  v_bg - variance between groups

/*  v_wg - variance within groups

/*  wa - weighting for aggregation, based on total variance

/*  vm - variance of the calculated mean

/*  mv - mean of finer scale variances

/*  n - effective number of measurements

/* NB - only calculating sample variances here, not variances of estimated means.

/* Also note that v0_bg is an uncertainty, not a sample variance

/* and v1_bg is total variances, but both are labelled as "between-group" to simplify the smoothing

h0 = %ingrid%

v0 = con(^ isnull(h0), sqr(%sd%))

vg0 = v0

w0 = con(isnull(v0), 0, 1.0 / v0)

wsq0 = sqr(w0)

n0 = con(^ isnull(h0), 1, 0)

&describe v0

&sv bigvar %grd$zmax%

setcell minof

/* aggregate to broader scales

&sv i 1

&sv done .false.

&describe h0

&do &until %done%

  &sv j [calc %i% - 1]

  &type Aggregate from %j% to %i%

  &describe h%j%

  &sv cell3 [calc %grd$dx% * 3]

  &describe h0

  &sv nx0 [round [calc %grd$xmin% / %cell3% - 0.5]]

  &sv ny0 [round [calc %grd$ymin% / %cell3% - 0.5]]

  &sv nx1 [round [calc %grd$xmax% / %cell3% + 0.5]]

  &sv ny1 [round [calc %grd$ymax% / %cell3% + 0.5]]

  &sv x0 [calc ( %nx0% - 0.5 ) * %cell3%]

  &sv y0 [calc ( %ny0% - 0.5 ) * %cell3%]

  &sv x1 [calc ( %nx1% + 0.5 ) * %cell3%]

  &sv y1 [calc ( %ny1% + 0.5 ) * %cell3%]

  setwindow %x0% %y0% %x1% %y1%

  w%i% = aggregate(w%j%, 3, sum)

  wsq%i% = aggregate(wsq%j%, 3, sum)

  n%i% = aggregate(n%j%, 3, sum)

  neff%i% = w%i% * w%i% / wsq%i%

  h%i% = aggregate(w%j% * h%j%, 3, sum) / w%i%

  vbg%i% = aggregate(w%j% * sqr(h%j% - h%i%), 3, sum) / w%i%

  &if %i% eq 1 &then vwg%i% = n%i% - n%i% /* zero, but with window and cell size set for us

  &else vwg%i% = aggregate(w%j% * vg%j%, 3, sum) / w%i%

  vg%i% = vbg%i% + vwg%i%

  vm%i% = 1.0 / w%i%

  mv%i% = n%i% / w%i%

  chisq%i% = 1 + %chisqa% / sqrt(neff%i% - 1) + %chisqb% / (neff%i% - 1)

  v%i% = con(vg%i% / chisq%i% < mv%i%, vm%i%, vg%i%)

  /* remove everything except h%i% and v%i%

  kill w%j%

  kill wsq%j%

  kill n%j%

  kill neff%i%

  kill vbg%i%

  kill vwg%i%

  kill vg%j%

  kill vm%i%

  kill mv%i%

  kill chisq%i%

  &sv done %i% eq 4

  &sv i [calc %i% + 1]

&end

&sv maxstep [calc %i% - 1]

&sv bigvar [calc %bigvar% * 10]

kill w%maxstep%

kill wsq%maxstep%

kill n%maxstep%

kill vg%maxstep%

/* smooth, refine and combine each layer in turn

copy h%maxstep% hs%maxstep%

copy v%maxstep% vs%maxstep%

kill h%maxstep%

kill v%maxstep%

setcell hs%maxstep%

setwindow hs%maxstep%

&do j := %maxstep% &to 1 &by -1

  &sv i [calc %j% - 1]

  &type Refine from %j% to %i%

  /* for the first stage where the coarser grid is refined and smoothed, set window to the coarse grid

  setcell h%i%

  setwindow maxof

  /* create smoothed higher resolution versions of h and v_bg, hopefully with no nulls!

  hs%j%_%i% = focalmean(hs%j%, circle, 2)

  vs%j%_%i% = focalmean(vs%j%, circle, 2)

  setcell h%i%

  &describe h%i%

  &sv cellsize %grd$dx%

  &describe bboxgrid

  setwindow [calc %grd$xmin% - 4 * %cellsize%] [calc %grd$ymin% - 4 * %cellsize%] [calc %grd$xmax% + 4 * %cellsize%] [calc %grd$ymax% + 4 * %cellsize%] h%i%

  /* create no-null version of finer h and v

  h%i%_c = con(isnull(h%i%), 0, h%i%)

  v%i%_c = con(isnull(v%i%), %bigvar%, v%i%)

  /* combine two values using least variance

  hs%i% = (h%i%_c / v%i%_c + hs%j%_%i% / vs%j%_%i% ) / (1.0 / v%i%_c + 1.0 / vs%j%_%i%)

  vs%i% = 1 / (1.0 / v%i%_c + 1.0 / vs%j%_%i%)

  kill v%i%_c  

  kill h%i%_c

  kill v%i%

  kill h%i%

  kill vs%j%_%i%

  kill hs%j%_%i%

  kill hs%j%

  kill vs%j%

&end

/* result is hs0, with variance vs0

kill bboxgrid

&args gridname radius geographic

&describe %gridname%

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

&else &sv xyfactor 1

&sv radiuscells = [calc %radius% / ( %grd$dx% * %xyfactor% )]

&type radius = %radiuscells% cells

&sv .maxcells 50000000

/* &sv .maxcells 500000 /* Testing on Kyeamba 

&sv .border [calc %radiuscells% + 1]

&args dem pctlname radius geographic

&type Running pctl %dem% %pctlname% %radius% %geographic%

&if [exists demtmp.flt -file] &then [delete demtmp.flt -file]

demtmp.flt = gridfloat(%dem%)

&if %geographic% &then &sys pctl -ud %radius% demtmp.flt pctl.flt

&else &sys pctl %radius% demtmp.flt pctl.flt

&if [delete demtmp.flt -file] ^= 0 &then &type ERROR deleting demtmp.flt

&if [delete demtmp.hdr -file] ^= 0 &then &type ERROR deleting demtmp.hdr

&if [delete demtmp.prj -file] ^= 0 &then &type ERROR deleting demtmp.prj

&if [exists %pctlname% -grid] &then kill %pctlname%

%pctlname% = floatgrid(pctl.flt)

&if [delete pctl.flt -file] ^= 0 &then &type ERROR deleting pctl.flt

&if [delete pctl.hdr -file] ^= 0 &then &type ERROR deleting pctl.hdr

/* merge 9 tiles in Oregon around Portland for testing methods

&if %:program%_ eq ARC_ &then grid

&if %:program%_ ne GRID_ &then

&return This program must be run from ARC or GRID.

&sv tiledir \\jupiter\\organisms\topo\tiles

&sv outdir \\jupiter\\organisms\topo\experimental

&sv name srtmv41

&sv vernum 1

/* This can be extended to larger areas (beyond 3x3 used here) but at some stage the

/* tilelist variable is going to get too long - AML has a limited line length

/* At this point, it would be better to tile by longitude band i.e. do a merge for

/* each value of w, then merge the resulting strips as a separate step

&sv tilelist

&do w = 122 &to 124

  &do n = 44 &to 46

    &sv tilelist %tilelist%, %tiledir%\w%w%\n%n%\%vernum%\%name%\w%w%n%n%

  &end

&end

%outdir%\srtmv41 = merge(%tilelist%)

:q

/* look for all .asc files in input dir

/* load them and split into 1 degree tiles

&if %:program%_ eq ARC_ &then grid

&if %:program%_ ne GRID_ &then

&return This program must be run from ARC or GRID.

&sv indir \\jupiter\\organisms\topo\incoming\srtmv41

&sv tiledir \\jupiter\\organisms\topo\tiles

&sv name srtmv41

&sv vernum 1

&do x = 19 &to 72

  &if %x% > 9 &then &sv xn %x%; &else &sv xn 0%x%

  &do y = 1 &to 36

    &if %y% > 9 &then &sv yn %y%; &else &sv yn 0%y%

    &sv ascfile %indir%\srtm_%xn%_%yn%.asc

    &if [exists %ascfile% -file] &then

    &do

      &type Got %ascfile%

      /* import the ASCII file 

      &sv tmpgrd [scratchname -dir -full]

      /* &sv tmpgrd c:\workspace\xx00000 - was using for testing

      /* &type Using %tmpgrd% as temporary grid

      %tmpgrd% = asciigrid(%ascfile%)

      /* determine x (long) and y (lat) range

      &describe %tmpgrd%

      &sv xmin %grd$xmin%

      &sv ymin %grd$ymin%

      &sv xmax %xmin% + 5

      &sv ymax %ymin% + 5

      &sv cell %grd$dx%

      &type %xmin% %ymin% %xmax% %ymax%

      &do x = %xmin% &to %xmax%

        &if %x% < 0 &then &sv xtile w[round [abs %x%]]; &else &sv xtile e[round %x%]

        &sv x1 %x% + 1

        &do y = %ymin% &to %ymax%

          &if %y% < 0 &then &sv ytile s[round [abs %y%]]; &else &sv ytile n[round %y%]

          &sv y1 %y% + 1

          setwindow %x% %y% %x1% %y1%

          &type Tile %x% %y%

          show setwindow

          &sv tmptile [scratchname -dir -full]

          &type Create tile %xtile%%ytile%

          %tmptile% = %tmpgrd%

          &describe %tmptile%

          &if %grd$stdv% > 0 &then &do

            /* grid is not all nodata

            &sv dirname %tiledir%\%xtile%

            &if ^ [exists %dirname% -directory] &then &sys mkdir %dirname%

            &sv dirname %tiledir%\%xtile%\%ytile%

            &if ^ [exists %dirname% -directory] &then &sys mkdir %dirname%

            &sv dirname %tiledir%\%xtile%\%ytile%\%vernum%

            &if ^ [exists %dirname% -directory] &then &sys mkdir %dirname%

            &sv dirname %tiledir%\%xtile%\%ytile%\%vernum%\%name%

            &if ^ [exists %dirname% -directory] &then &sys mkdir %dirname%

            &sv gridname %dirname%\%xtile%%ytile%

            &if [exists %gridname% -grid] &then &do

              &type WARNING: did not expect to find %gridname%!

            &end

            &else &do

              copy %tmptile% %dirname%\%xtile%%ytile%

            &end

          &end

          kill %tmptile%

        &end

      &end

      kill %tmpgrd%

    &end

  &end 

&end

import os, fnmatch, zipfile, datetime

outpath = '//jupiter/organisms/topo/incoming/srtmv41'

inpath = '//jupiter/organisms/SRTM_90m_ASCII_v4.1'

logname = outpath + '/unzip.log'

logfile = open(logname, 'a')

logfile.write('Start unziptiles.py at ' + str(datetime.datetime.today()) + '\n')

nprocessed = 0

nerrors = 0

try:

    for filename in os.listdir(inpath):

        fullname = inpath + '/' + filename

        if fnmatch.fnmatch(fullname, '*.zip'):

            print fullname

            try:

                zip = zipfile.ZipFile(fullname, 'r')

                for zipi in zip.infolist():

                    name = zipi.filename

                    if name != 'readme.txt':

                        outfilename = outpath + '/' + name

                        print outfilename

                        if os.path.exists(outfilename):

                            print 'done'

                        else:

                            print 'Process', outfilename

                            outfile = open(outfilename, 'w')

                            outfile.write(zip.read(name))

                            outfile.close()

                            nprocessed = nprocessed + 1

                zip.close()

            except Exception, e:

                print str(e)

                logfile.write(str(e) + '\n')

                nerrors = nerrors + 1

except Exception, e:

    print "Exception:", e.str()

    logfile.write(str(e) + '\n')

logfile.write('Processed ' + str(nprocessed) + '\n')

logfile.write('Errors: ' + str(nerrors) + '\n')

logfile.write('End unziptiles.py at ' + str(datetime.datetime.today()) + '\n')

logfile.close()

/* need to ensure that the aggregation boxes are consistent across the whole set of layers

/* so might need a setwindow here to enforce the correct location

/* Question about 1km Behrman version - should that be projected from 30" stats, or

/* should data be projected to Behrman and aggregated over 1km cell. Far from the equator

/* these will give quite different results

/* elevation - derive mean, min, max, std. dev, skew

elev_mean = aggregate(srtm_or_smth, 10, mean)

elev_min = aggregate(srtm_or_smth, 10, min)

elev_max = aggregate(srtm_or_smth, 10, max)

elev_range = elev_max - elev_min

/* stupid arc doesn't have a standard deviation function for aggregation

elev_sd = sqrt(aggregate(sqr(srtm_or_smth - elev_mean), 10, sum) / 100)

/* slope - just mean and standard deviation

slopedeg_mean = aggregate(slopedeg, 10, mean)

slopedeg_sd = sqrt(aggregate(sqr(slopdeg - slopedeg_mean), 10, sum) / 100)

/* wetness index - mean and standard deviation

/* twi not done yet

twi_fill_mean = aggregate(twi_fill, 10, mean)

twi_fill_sd = sqrt(aggregate(sqr(twi_fill - twi_fill_mean), 10, sum) / 100)

/* Invoke another AML on a grid to produce a grid

/* Other AML takes arguments: gridname outgridname otherargs

&args targetaml gridname outgridname otherargs:REST

/* targetaml-limits.aml takes same arguments as target.aml (except for

/* output file) and sets maximum number of cells and overlap required

&run %targetaml%-limits %gridname% [unquote %otherargs%]

&describe %gridname%

&sv gridxmin = %grd$xmin%

&sv gridxmax = %grd$xmax%

&sv gridymin = %grd$ymin%

&sv gridymax = %grd$ymax%

&sv gridncols = %grd$ncols%

&sv gridnrows = %grd$nrows%

&sv cellsize = %grd$dx%

/****** TEMPORARY *******/

/* setcell %cellsize%

/**** END TEMPORARY *****/

&sv ncells = %grd$ncols% * %grd$nrows%

&type %ncells% cells in %gridname%

&if %ncells% <= %.maxcells% &then &do

  /* no need to do chunks, just invoke target AML

  &type Process as single chunk

  &run %targetaml% %gridname% %outgridname% [unquote %otherargs%]

  &return

&end

/* need to do chunks

&type Process in multiple chunks

/* &type X: %gridxmin% %gridxmax% (%gridncols%) 

/* &type Y: %gridymin% %gridymax% (%gridnrows%) 

/* &type %ncells% cells, %.maxcells% max, %.border% border

&sv nchunks = [truncate [calc %ncells% / %.maxcells%]] + 1

&type Need at least %nchunks% chunks

&sv chunkdim = [sqrt %.maxcells%] + 2 * %.border%

&if %gridncols% < %chunkdim% &then

  &do

    &sv colchunks = 1

    &sv rowchunks = %nchunks%

    &sv chunkncols = %gridncols%

    &sv chunknrows = [truncate [calc %gridnrows% / %rowchunks%]] + 1

  &end

&else &if %gridnrows% < %chunkdim% &then

  &do

    &sv rowchunks = 1

    &sv colchunks = %nchunks%

    &sv chunkncols = [truncate [calc %gridncols% / %colchunks%]] + 1

  &end

&else 

  &do

    &sv rowchunks = [truncate [calc %gridnrows% / %chunkdim%]] + 1

    &sv colchunks = [truncate [calc %gridncols% / %chunkdim%]] + 1

    &sv chunkncols = [truncate [calc %gridncols% / %colchunks%]] + 1

    &sv chunknrows = [truncate [calc %gridnrows% / %rowchunks%]] + 1

  &end

&type %colchunks% colchunks, %rowchunks% rowchunks

&type Chunk size %chunkncols% cols by %chunknrows% rows

&sv qtrcell = %cellsize% / 4

&sv halfcell = %cellsize% / 2

&sv colist

&sv chunki = 1

&do &while %chunki% <= %colchunks%

  &sv chunkj = 1

  &sv colist_col

  &do &while %chunkj% <= %rowchunks%

    &type Chunk %chunki% %chunkj%

    &sv coxmin = %gridxmin% + ( %chunki% - 1 ) * %chunkncols% * %cellsize% - %qtrcell%

    &sv coxmax = %coxmin% + %chunkncols% * %cellsize% + %halfcell%

    &sv coymin = %gridymin% + ( %chunkj% - 1 ) * %chunknrows% * %cellsize% - %qtrcell%

    &sv coymax = %coymin% + %chunknrows% * %cellsize% + %halfcell%

    /* &type Output X: %coxmin% %coxmax%

    /* &type Output Y: %coymin% %coymax%

    &if %chunki% = 1 &then &sv chunkxmin = %coxmin%

    &else &sv chunkxmin = %coxmin% - %.border% * %cellsize%

    &if %chunki% = %colchunks% &then &sv chunkxmax = %coxmax%

    &else &sv chunkxmax = %coxmax% + %.border% * %cellsize%

    &if %chunkj% = 1 &then &sv chunkymin = %coymin%

    &else &sv chunkymin = %coymin% - %.border% * %cellsize%

    &if %chunkj% = %rowchunks% &then &sv chunkymax = %coymax%

    &else &sv chunkymax = %coymax% + %.border% * %cellsize%

    /* &type Chunk X: %chunkxmin% %chunkxmax%

    /* &type Chunk Y: %chunkymin% %chunkymax%

    /* clip input file to chunk (with borders)

    setwindow %chunkxmin% %chunkymin% %chunkxmax% %chunkymax% %gridname%

    &if [exists chunk -grid] &then kill chunk

    chunk = %gridname%  

    &if [exists co_tmp -grid] &then kill co_tmp

    &run %targetaml% chunk co_tmp [unquote %otherargs%]

    kill chunk

    /* clip output file to chunk (without borders)

    setwindow %coxmin% %coymin% %coxmax% %coymax% %gridname%

    &if [exists co_%chunki%_%chunkj% -grid] &then kill co_%chunki%_%chunkj%

    co_%chunki%_%chunkj% = co_tmp

    kill co_tmp

    &if %chunkj% = 1 &then &sv colist_col = co_%chunki%_%chunkj%

    &else &sv colist_col = %colist_col%, co_%chunki%_%chunkj%

    &sv chunkj = %chunkj% + 1

  &end

  setwindow maxof

  &if [exists co_%chunki% -grid] &then kill co_%chunki%

  co_%chunki% = merge(%colist_col%)

  kill (!%colist_col%!)

  &if %chunki% = 1 &then &sv colist = co_%chunki%

  &else &sv colist = %colist%, co_%chunki%

  &sv chunki = %chunki% + 1

&end

setwindow maxof

&if [exists %outgridname% -grid] &then kill %outgridname%

%outgridname% = merge(%colist%)

kill (!%colist%!)

11 11

0.017 0.046 0.100 0.174 0.242 0.271 0.242 0.174 0.100 0.046 0.017 

0.046 0.124 0.271 0.472 0.659 0.736 0.659 0.472 0.271 0.124 0.046 

0.100 0.271 0.590 1.028 1.434 1.603 1.434 1.028 0.590 0.271 0.100 

0.174 0.472 1.028 1.791 2.500 2.793 2.500 1.791 1.028 0.472 0.174 

0.242 0.659 1.434 2.500 3.488 3.898 3.488 2.500 1.434 0.659 0.242 

0.271 0.736 1.603 2.793 3.898 4.356 3.898 2.793 1.603 0.736 0.271 

0.242 0.659 1.434 2.500 3.488 3.898 3.488 2.500 1.434 0.659 0.242 

0.174 0.472 1.028 1.791 2.500 2.793 2.500 1.791 1.028 0.472 0.174 

0.100 0.271 0.590 1.028 1.434 1.603 1.434 1.028 0.590 0.271 0.100 

0.046 0.124 0.271 0.472 0.659 0.736 0.659 0.472 0.271 0.124 0.046 

0.017 0.046 0.100 0.174 0.242 0.271 0.242 0.174 0.100 0.046 0.017 

/* using srtmv41 3x3 degree patch smoothed to srtmv41_smth

/* need to run GRID

&describe srtm_or_smth

&sv cellsize %grd$dx% * 111000

dzdx = (srtm_or_smth(1,0) - srtm_or_smth(-1,0)) / (2 * %cellsize% * cos($$ymap / deg))

dzdy = (srtm_or_smth(0,-1) - srtm_or_smth(0,1)) / (2 * %cellsize%)

slope = sqrt(sqr(dzdx) + sqr(dzdy))

slopedeg = atan(slope) * deg

/* use scaling factor that is average of x and y scaling factors, 0.707 and 1 ~= 0.85

/* and 0.85 / 111000 ~ 0.0000107

arcslope = slope(srtm_or_smth, 0.0000107, degree)

hillshade = hillshade(srtm_or_smth, 300, 30, shade, 0.00001)

aspect1 = atan2(0 - dzdy, 0 - dzdx) * deg

aspect = con(aspect1 > 90, 450 - aspect1, 90 - aspect1)

kill aspect1

/*use elevation 

flowdir = flowdirection(srtm_or_smth)

/*cell area to account for latitudinal area effects, smallest area near pole

setwindow srtm_or_smth

setcell srtm_or_smth 

cellarea = sqr(%cellsize%) * cos($$ymap / deg)

flowacc = flowaccumulation(flowdir, cellarea)

/*account for sinks

/*change # defaulto z limit

/*caution fill is a command not function

/*fill srtm_or_smth srtm_or_fill sink # flowdir_fill

flowacc_fill = flowaccumulation(flowdir_fill, cellarea)

/*scatch_fill is specific catchment area in m

/*wrong filling? flowdir_fill produces artefacts

scatch_fill = flowacc_fill / sqrt(cellarea) + sqrt(cellarea) / 2

/*setting the min slope to 0.1

adjslope = tan(max(slopedeg, 0.1) / deg)

twi_fill = ln(scatch_fill / adjslope)

/*   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