Merging DEMs

merging DEMs

This page described Cshell scripts and ARC/INFO macros they create to merge a large number of 7.5' digital elevation models. This is written primarily for Dan Saul, and for my own records. The scripts have been tested under the Cshell under solaris 8 unix. If they misbehave on your system, or if you just want to see how they work, you can cut out pieces of interesting lines and paste them into your command line.

I have hundreds of DEMs posted at http://gis.ess.washington.edu/data/raster/ This examples shows how I import and mosaic the 216 10-meter quads currently in stock for zone 11 (east of 120°) Washington State. I have a directory for each 30'x60' quadrangle, containing up to 64 files for 7.5' quads. (Actually, the northern quads also contain over-edge files above 49°.) The .zip files expand into single files with a .dem extension. For example, /job/www/tenmeter/byquad/ritzville/q941.zip expands into q941.dem. The naming scheme follows the DNR convention: this is the ninth row down from 49° and the 41st column east from 125°. This is a standard ASCII USGS file of the type that is gradually being supplanted by the SDTS format. I curently have no 10-meter data for Pendleton or Spokane quads, so I pulled them out of the list. It is necessary that this be a nested loop, or the list within parentheses would be too long. I call the following file unzipcom, and I launch it from the csh prompt with "source unzipcom". You can also copy and paste it to the command line. This script also creates a file named import.aml.

cat /dev/null >! import.aml
#foreach dir (okanagan ritzville wallawalla pendleton sandpoint spokane pullman grange)
foreach dir (okanagan ritzville wallawalla sandpoint pullman grange)
  foreach fil (/job/www/tenmeter/byquad/$dir/*.zip)
    unzip $fil
    set quad=$fil:t
    echo $quad:r = demgrid\($quad:r.dem\) >> import.aml
  end
end

This will unzip all the .dem files into the current directory. It creates import.aml, which looks like

q141 = demgrid(q141.dem)
q142 = demgrid(q142.dem)
...

The problem with import.aml is that ARC/INFO fails to close certain files, so it will crash if it is not restarted several times. We therefore use the split command to break import.aml into pieces 53 lines long. We create a script named callarccom and source it to call ARC/INFO repeatedly to import the DEMs. I have chosen to use the demgrid function, which imports these files in integer decimeters. The latticedem command would import them as floating-point meters.

split -53 import.aml import
echo "echo importing .dem files to integer grids" >! callarccom
foreach aml (importa*)
  echo quit >> $aml
  echo arc \'grid  $aml\' >> callarccom
end
echo calling callarccom
source callarccom

In my case, the automatically-created file callarccom looks like

arc 'grid importaa'
arc 'grid importab'
arc 'grid importac'
arc 'grid importad'

This leaves us with a large number (216 in my case) DEMs, but ARC/INFO can merge only 12 at a time. The following script creates amls to create a hierarchy of intermediate grids.

ls *.dem | sed 's/[^0-9]//g' | sort -n | sed 's/^/q/' | paste - - - - - - - - - - - - | tr '\011' ',' | cat -n | awk '{print "tmp"$1" = merge("$2")"}' | sed 's/,,*)/)/' > mergea.aml
echo You can delete the .dem files now
while (1 < `cat mergea.aml | wc -l`)
  cut -d ' ' -f1 mergea.aml | paste - - - - - - - - - - - - | tr '\011' ',' | cat -n | awk '{print "tmpb"$1" = merge("$2")"}' | sed 's/,,*)/)/' >! mergeb.aml
  break
end
while (1 < `cat mergeb.aml | wc -l`)
  cut -d ' ' -f1 mergeb.aml | paste - - - - - - - - - - - - | tr '\011' ',' | cat -n | awk '{print "tmpc"$1" = merge("$2")"}'| sed 's/,,*)/)/' >! mergec.aml
  break
end
echo "/*merging amls" >! merge.aml
echo grid >> merge.aml
echo "&run mergea.aml" >> merge.aml
if (-e mergeb.aml)then
  echo "&run mergeb.aml" >> merge.aml
  if (-e mergec.aml)then
    echo "&run mergec.aml" >> merge.aml
    echo "rename tmpc1 mergeddem" >> merge.aml
  else
    echo "rename tmpa1 mergeddem" >> merge.aml
  endif
else
  echo "rename tmp1 mergeddem" >> merge.aml
endif
arc "&run merge.aml"

Aren't those while loops strange? I couldn't find any way to get the Cshell to do something like
if(`cat mergeb.aml | wc -l` > 1)
In my case, mergea.aml looks like

tmp1 = merge(q42,q43,q44,q45,q46,q47,q48,q49,q50,q51,q52,q53)
tmp2 = merge(q54,q55,q56,q57,q58,q59,q60,q61,q62,q141,q142,q143)
tmp3 = merge(q144,q145,q146,q147,q148,q149,q150,q151,q152,q153,q154,q155)
tmp4 = merge(q156,q157,q158,q159,q160,q161,q162,q163,q164,q241,q242,q243)
tmp5 = merge(q244,q245,q246,q247,q248,q249,q250,q251,q252,q253,q254,q255)
tmp6 = merge(q256,q257,q258,q259,q260,q261,q262,q263,q264,q341,q342,q343)
tmp7 = merge(q344,q345,q346,q347,q348,q349,q350,q351,q352,q353,q354,q355)
tmp8 = merge(q356,q357,q358,q359,q360,q361,q362,q363,q364,q441,q442,q443)
tmp9 = merge(q444,q445,q446,q447,q448,q449,q450,q451,q452,q453,q454,q455)
tmp10 = merge(q456,q457,q458,q459,q460,q461,q462,q463,q464,q541,q542,q543)
tmp11 = merge(q544,q545,q546,q547,q548,q549,q550,q551,q552,q553,q554,q555)
tmp12 = merge(q556,q557,q558,q559,q560,q561,q562,q563,q564,q641,q642,q643)
tmp13 = merge(q644,q645,q646,q647,q648,q649,q650,q651,q652,q653,q654,q655)
tmp14 = merge(q656,q657,q658,q659,q660,q661,q662,q663,q664,q741,q742,q743)
tmp15 = merge(q759,q760,q761,q762,q763,q764,q841,q842,q843,q941,q1041,q1741)
tmp16 = merge(q1841,q1846,q1847,q1941,q1946,q1947,q1948,q1949,q1950,q2041,q2258,q2259)
tmp17 = merge(q2260,q2261,q2262,q2358,q2359,q2360,q2361,q2362,q2363,q2456,q2457,q2458)
tmp18 = merge(q2459,q2460,q2461,q2462,q2463,q2464,q2465,q2559,q2560,q2562,q2564,q2565)

Files are processed by number, so the no-data aras are not too exteme. The script is written for 2 levels of the hierachy, so it can handle up to 1728 components, more than any reasonable-sized state contains. Cleanup involves deleting all the .dem files and killing all the tmp grids.

/usr/bin/ls -d tmp* q* | sed 's/^/kill /' >! delete.aml
/pre>
Now we can check for no_data slivers between quads.


tmpmask = con(^ isnull(mergeddem),1)
testpoly = gridpoly(tmpmask)
linepen 0.1 ROUND
arcs testpoly

This will show a dark line at any slivers. To fix trouble areas:

setwindow * mergeddem
patch = focalmean(mergeddem) /* the default params should cover any gap
setwindow * mergeddem
patch2  = focalmean(mergeddem)  /* and so forth
setwindow mergeddem
rename mergeddem uffda
mergeddem = merge(uffda,patch1,patch2)

Web Curator: Harvey Greenberg
about this site

merging DEMs