[ncl-talk] Computing Aspect

[Dennis Shea]

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
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180313/b0aa98ca/attachment.html>