/* Program: dem2hillsh.aml /* /* Purpose: Creates ARC elevation grid and hillshade grid from a DEM file. /* /* Usage: &r dem2hillsh /* /* {clip_cover} /* {datum} {projection} {horiz_units} {UTM_Zone} /* Ex: &r dem2hillsh gol.dem elev hillshade 10 315 35 5 0 1.682 goln ~ /* nad27 UTM meters 13 /* &r dem2hillsh mor.dem morz hillsh 30 310 40 4 10 -0.65 # NAD83 utm /* /* Output: Outputs are 1) a raw elevation grid of the same resolution as the /* input DEM, and 2) a hillshaded elevation grid with cellsize /* of , clipped to the extent of , rotated an /* angle of (positive angle = clockwise, negative = counter), /* with shaded relief parameters: altitude of , azimuth of , /* and vertical exaggeration of . /* /* Notes: Run this program without arguments to see the usage and examples. /* /* The is the size, in cells, of the square /* filter which will be passed over the hillshade grid (use 0 for no /* smoothing, use 10 for significant smoothing, any integer okay). /* If an angle of 0 is entered, the hillshade grid will not be rotated. /* Use the normal "#" character to skip optional arguments. /* If no UTM zone is passed in, no zone information will exist /* in the projection file. /* /* History: JKosovich 03/05/97 Coding for abandoned mine lands, based on /* AMLs image2grid, makeshade, drg2grid, etc. /* JKosovich 09/10/97 Renamed from old dem2grid.aml /* JKosovich 10/10/97 Cleans up leftover Grid trash (xyxy001...) /* JKosovich 01/13/98 Turned watch off, loop to clean up grid junk /* JKosovich 02/04/98 Added argument /* /******************************************************************************* /* &args indem outelevgrd outshgrid outshres az alt ve xy angle clipcov datum ~ proj units zone &severity &error &routine bailout &s prog = [ENTRYNAME %aml$file% -NOEXT -FILE] /* current AML being run &watch %prog%.watch /* /***** Check for correct usage and correct argument types ***** /* &if [NULL %angle%] &then &do &s error_case = 1 &call error_cond &end &if [LENGTH %outelevgrd%] > 12 | [LENGTH %outshgrid%] > 12 &then &do &s error_case = 2 &call error_cond &end &do xx &list outshres angle ve /* these can be int or real &if [TYPE [VALUE %xx%]] > 0 &then &do &s error_case = 3 &call error_cond &end &end &do xx &list alt az xy zone /* these can only be integer &if ^ [NULL [VALUE %xx%]] &then &do &if [TYPE [VALUE %xx%]] ^= -1 &then &do &s error_case = 4 &call error_cond &end &end &end /* /***** Set toggle options ***** /* &s smooth = .TRUE. /* Don't smooth if 0 xy &if %xy% = 0 &then &s smooth = .FALSE. &s rotate = .TRUE. /* Don't rotate if 0 ang &if %angle% = 0 &then &s rotate = .FALSE. &if %angle% = 0.0 &then &s rotate = .FALSE. &do xx &list clipcov datum proj units zone /* if #, set to null &if ^ [NULL [VALUE %xx%]] &then &do &if [TYPE [VALUE %xx%]] = 1 &then /* if character &do &if [VALUE %xx%] = # &then &s %xx% = &end &end &end &s projdef = .TRUE. &if [NULL %proj%] | [NULL %datum%] &then &s projdef = .FALSE. &s clip = .TRUE. &if [NULL %clipcov%] &then &s clip = .FALSE. /* /***** Check for existence of input file and cover ***** /* &if ^ [EXISTS %indem% -FILE] &then &do &s error_case = 5 &call error_cond &end &if %clip% &then &do &if ^ [EXISTS %clipcov% -COVER] &then &do &s error_case = 6 &call error_cond &end &else &if ^ [EXISTS %clipcov% -POLYGON] &then &do &s error_case = 7 &call error_cond &end &end /* /***** Set terminal ***** /* &if [NULL [EXTRACT 1 [SHOW &terminal]]] &then &terminal 9999 /* /***** Convert DEM to Arc grid and convert to integer if necessary ***** /* &if [EXISTS %outelevgrd% -GRID] &then KILL %outelevgrd% ALL &if [EXISTS xxgrd1 -GRID] &then KILL xxgrd1 ALL DEMLATTICE %indem% %outelevgrd% USGS &describe %outelevgrd% &s indemres = %grd$dx% /* resolution of incoming DEM = latt res &s gridtype = %grd$type% /* type of grid: 1 = int, 2 = floating pt. &if %gridtype% = 2 &then &do DISPLAY 0 GRID xxgrd1 = INT(%outelevgrd%) /* Convert to integer, much smaller file size! QUIT &end &if [EXISTS xxgrd1 -GRID] &then &do KILL %outelevgrd% ALL RENAME xxgrd1 %outelevgrd% &end /* /***** Define projection for grid ***** /* &if %projdef% &then &do PROJECTDEFINE GRID %outelevgrd% PROJECTION %proj% DATUM %datum% &if ^ [NULL %units%] &then UNITS %units% &if ^ [NULL %zone%] &then ZONE %zone% PARAMETERS &if [UPCASE [SHOW program]] = PRJDEF &then END &end /* /***** Resample grid ***** /* &if [EXISTS xxgrd2 -GRID] &then KILL xxgrd2 ALL &if %indemres% ^= %outshres% &then &do DISPLAY 0 GRID xxgrd2 = RESAMPLE(%outelevgrd%,%outshres%) QUIT &end &if [EXISTS xxgrd2 -GRID] &then PROJECTCOPY GRID %outelevgrd% GRID xxgrd2 &else COPY %outelevgrd% xxgrd2 /* /***** Clip grid to %clipcov% extent ***** /* &if [EXISTS xxgrd3 -GRID] &then KILL xxgrd3 ALL &if %clip% &then LATTICECLIP xxgrd2 %clipcov% xxgrd3 MINIMUM 1 &if [EXISTS xxgrd3 -GRID] &then &do PROJECTCOPY GRID xxgrd2 GRID xxgrd3 KILL xxgrd2 ALL &end &else RENAME xxgrd2 xxgrd3 /* /***** Rotate grid to due N-S orientation for plotting ***** /* &if [EXISTS xxgrd4 -GRID] &then KILL xxgrd4 ALL &if %rotate% &then GRIDROTATE xxgrd3 xxgrd4 %angle% &if [EXISTS xxgrd4 -GRID] &then &do PROJECTCOPY GRID xxgrd3 GRID xxgrd4 KILL xxgrd3 ALL &end &else RENAME xxgrd3 xxgrd4 /* /***** Apply HILLSHADE parameters ***** /* &if [EXISTS %outshgrid% -GRID] &then KILL %outshgrid% ALL DISPLAY 0 GRID &if %smooth% &then %outshgrid% = HILLSHADE(FOCALMEAN(xxgrd4,RECTANGLE,%xy%,%xy%),%az%,%alt%,SHADE,%ve%) &else %outshgrid% = HILLSHADE(xxgrd4,%az%,%alt%,SHADE,%ve%) QUIT &if [EXISTS %outshgrid% -GRID] &then &do PROJECTCOPY GRID xxgrd4 GRID %outshgrid% KILL xxgrd4 ALL &type &type &type *** %prog% successful. Shaded relief grid is %outshgrid% &type &type &end &else &do &type &type &type *** %prog% ERROR: Shaded relief grid %outshgrid% not created!!! &type *** Last temporary grid made was: xxgrd4. &type &type &end /* /***** Clean up Grid junk ***** /* &messages &off &all &do n = 1 &to 9 &if [EXISTS xyxy00%n% -GRID] &then KILL xyxy00%n% ALL &end &messages &on /* /* &watch &off &return /* /* /*------------------------------------------------------------------------------ /*------------------------------ ROUTINE LIST -------------------------------- /*------------------------------------------------------------------------------ /* /* /*------------------ &routine ERROR_COND /*------------------ /* &severity &error &routine bailout &if [UPCASE [SHOW PROGRAM]] = TABLES &then QUIT &messages &on &type &type *************** PROGRAM [TRANSLATE %prog%] HALTED ***************** /* &if %error_case% = 1 &then /***** Incorrect usage &do &type USAGE: &r %prog% &type &type {clip_cover} &type {datum} {projection} {horiz_units} {UTM_Zone} &type &type Ex: &r %prog% gol.dem elev hillsh 10 315 35 5 0 1.682 goln &type nad27 UTM meters 13 &type &r %prog% abc.dem abcelv hillsh 30 310 40 4 9 -0.65 # NAD83 utm &type &type OUTPUT ELEV and HILLSHADE GRID NAMES must be 12 characters or less. &type OUTPUT HILLSHADE RESOLUTION is in ground units (eg. meters). &type AZIMUTH is whole degrees clockwise from North, North being 0. &type ALTITUDE is whole degrees above the horizon, horizon being 0. &type VERTICAL EXAGGERATION can be integer or whole, typically 5 or less. &type SMOOTHING KERNEL: use 0 for no hillshade smoothing, any integer OK. &type Use positive ANGLE to rotate grid clockwise to grid north, &type Use negative ANGLE to rotate grid counterclockwise. &type Use ANGLE of 0 for no rotation. &type If used, CLIP COVER must have POLYGON topology. &type If used, DATUM is the input data's datum (upper/lower case is okay). &type If used, PROJECTION is the input data projection. &type If used, UNITS is the the input data horizontal units. &type If used, ZONE is the input data projection zone. &end &if %error_case% = 2 &then /***** If output grid name too long &type Output grid name is too long - must be 12 characters or less!!! &if %error_case% = 3 &then /***** Argument must be int or real &type Argument %xx% must be an integer or real value!!! &if %error_case% = 4 &then /***** Argument must be integer &type Argument %xx% must be an integer value!!! &if %error_case% = 5 &then /***** If dem file not found &type DEM file %indem% does not exist!!! &if %error_case% = 6 &then /***** If clip cover not found &type Clip coverage %clipcov% does not exist!!! &if %error_case% = 7 &then /***** If clip cover not a POLY cover &type Clip coverage %clipcov% does not have POLYGON topology!!! /* &type Process halted [DATE -DOW] [DATE -CAL] [DATE -AMPM] &type ********************************************************* &type /* /* /***** Clean up Grid junk ***** /* &messages &off &all &do n = 1 &to 9 &if [EXISTS xyxy00%n% -GRID] &then KILL xyxy00%n% ALL &end &messages &on &watch &off &return; &return &error /* /* /*------------------ &routine BAILOUT /*------------------ /* &severity &error &ignore /* /***** Clean up Grid junk ***** /* &messages &off &all &do n = 1 &to 9 &if [EXISTS xyxy00%n% -GRID] &then KILL xyxy00%n% ALL &end &messages &on &type &type *** Error occurred at line %aml$errorline% in %aml$errorfile% &watch &off &return; &return &error /* /*