[ncl-talk] Fwd: Computing Aspect

[Dennis Shea]

Sent from my iPhone

Begin forwarded message:

> From: Michael Notaro <mnotaro at wisc.edu>
> Date: March 14, 2018 at 2:11:53 PM MDT
> To: Dennis Shea <shea at ucar.edu>
> Subject: Re: [ncl-talk] Computing Aspect
> 
> Dennis, thanks for your feedback and suggestions.
> I believe the code below works.  I tried different
> configurations of simple models, facing N, NE, E, SE, S,
> SW, W, NW, and it looked right each time.  
> Michael
> 
> ; a b c
> ; d e f
> ; g h i
> 
> a=5000.
> b=3000.
> c=1000.
> d=5000.
> e=3000.
> f=1000.
> g=5000.
> h=3000.
> i=1000.
> 
> dx=800.
> dy=800.
> 
> dzdy=((a+2.*b+c)-(g+2.*h+i))/(2.*dy)
> dzdx=((c+2.*f+i)-(a+2.*d+g))/(2.*dx)
> aspect=180.+57.29578*atan2(dzdx,dzdy)
> 
> print(dzdy)
> print(dzdx)
> print(aspect)
> 
> 
> Michael Notaro
> Associate Director
> Nelson Institute Center for Climatic Research
> University of Wisconsin-Madison
> Phone: (608) 261-1503
> Email: mnotaro at wisc.edu
> 
> 
>  
> From: Dennis Shea <shea at ucar.edu>
> Sent: Wednesday, March 14, 2018 9:29 AM
> To: Michael Notaro
> Subject: Re: [ncl-talk] Computing Aspect
>  
> Sigh ...  offline ...
> 
> atan2 has angles  -180 to 180 degrees with those 2 values relative to the x- axis (mathematical)
> 
> In arcmap
> aspect is < 90 then .... direction=90-aspect         which is an azimuth with respect to North from 0-360
> 
> aspect is > 90 then .... direction = 360-aspect     ditto
> 
> 
> You may have to play with this.
> 
> 
> ---
> 
> When you get a function that works, send it to me and I'll make a funcytion out of it.
> 
> 
> Cheers
> 
> D
> 
> On Tue, Mar 13, 2018 at 11:27 PM, Dennis Shea <shea at ucar.edu> wrote:
> 
> 
> I think this:
> 
> dzdx = ( (c + 2*f + i) - (a + 2*d + g))/(2*dx)   ; right - left
> ----
> dzdy = ( (a + 2*b + c) - (g + 2*h + i))/(2*dy)   ; top - bottom
> ----
> aspect = 57.29578*atan2(dzdy,dzdx)   
> 
> if (aspect.lt.0) then
>     aspect = 90 - aspect
> end if
> 
> ---
> This link has a formula that uses atan ,,, not atan2
> 
> http://www.spatialanalysisonline.com/HTML/index.html?gradient__slope_and_aspect.htm
> 
> D
> 
> On Tue, Mar 13, 2018 at 10:51 PM, Dennis Shea <shea at ucar.edu> wrote:
> 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
> ncl-talk Info Page - University Corporation for ...
> mailman.ucar.edu
> This email list is for NCL users to exchange information, ask questions, and report bugs on NCL or other related software if appropriate. It is also for the ...
> 
> 
> 
> 
> 
> 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180314/5e7bf5f5/attachment.html>