[ncl-talk] Computing Aspect

[Michael Notaro]

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
University of Wisconsin-Madison
Phone: (608) 261-1503
Email: mnotaro at wisc.edu
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180312/534d9d7a/attachment.html>