# [ncl-talk] Computing Aspect

Michael Notaro mnotaro at wisc.edu
Mon Mar 12 15:27:38 MDT 2018

```I used the following to compute terrain aspect based

on elevation.  But it seems the values are possibly not

quite right, since they disagree with another aspect

dataset I have.  Does anyone see something wrong below?

Michael

ny=dims(0)

nx=dims(1)

aspect=new((/ny,nx/),"float")

aspect at _FillValue=-999.

aspect=-999.

; a b c

; d e f

; g h i

do ii=1,nx-2

do jj=1,ny-2

if (.not.ismissing(elev(jj,ii))) then

a=elevation(jj+1,ii-1)

b=elevation(jj+1,ii)

c=elevation(jj+1,ii+1)

d=elevation(jj,ii-1)

e=elevation(jj,ii)

f=elevation(jj,ii+1)

g=elevation(jj-1,ii-1)

h=elevation(jj-1,ii)

i=elevation(jj-1,ii+1)

dzdx=((a+2.*d+g)-(c+2.*f+i))/(8.*800.) ; 800 m PRISM

dzdy=((g+2.*h+i)-(a+2.*b+c))/(8.*800.)

aspect(jj,ii)=57.29578*atan2(dzdy,dzdx)

if (aspect(jj,ii).le.0.) then

aspect(jj,ii)=90.-aspect(jj,ii)

else if (aspect(jj,ii).gt.90.) then

aspect(jj,ii)=360.-aspect(jj,ii)+90.

else

aspect(jj,ii)=90.-aspect(jj,ii)

end if

end if

end if

end do

end do

Michael Notaro
Associate Director
Nelson Institute Center for Climatic Research