[ncl-talk] Computing Aspect

[Dennis Shea]

Hi Mike,

Maybe a minus sign in the 2nd argument of atan2?

Given your:

; a b c

; d e f

; g h i


dzdx = ( (c + 2*f + i) - (a + 2*d + g))/dx

----

dzdy = ( (a + 2*b + c) - (g + 2*h + i))/dy


or, maybe, not sure of sign convention


dzdy = ( (g + 2*h + i) - (a + 2*b + c))/dy

----

Your code may need a minus sign before the 2nd argument:


aspect = 57.29578*atan2(dzdy,-dzdx)   ; <=== minus sign ???


if (aspect.lt.0) then

    aspect = 90 - aspect

end if


------------------


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

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

end if

--


If you use a minus sign, maybe you do not need other 'if' blocks?


On Mon, Mar 12, 2018 at 3:27 PM, Michael Notaro <mnotaro at wisc.edu> wrote:

> 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
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180313/e0480d06/attachment.html>