; This program was developed by Camilo D. Rennó and Lise C. Banon using IDL/ENVI language. This code can obtained in http://urlib.net/sid.inpe.br/mtc-m19/2012/09.03.01.24 ; Version: 1.0 (sep 2012) ; This program calculates attributes based on morphometric and flow direction from a digital elevation model (DEM) ; Input: DEM and LDD grids (LDD grid is coded according to Rennó et al, 2008) ; Output: grid with 23 attributes ; References: ; Wood, J.D. The geomorphological characterisation of digital elevation models. 1996. PhD Thesis, University of Leicester, UK, http://www.soi.city.ac.uk/~jwo/phd ; Shary et al. Fundamental quantitative methods of land surface analysis. Geoderma, 107: 1-32, 2002. ; Rennó et al. HAND, a new terrain descriptor using SRTM-DEM: mapping terra-firme rainforest environments in Amazonia. Remote Sensing of Environment, 112: 3469-3481, 2008. PRO DEM_attributes Envi_Select,Title='Input DEM Grid',Fid=fid,Dims=dims,Pos=pos,/Band_Only if fid eq -1 then return ns=dims[2]-dims[1]+1 nl=dims[4]-dims[3]+1 Envi_Select,Title='Input LDD Grid',Fid=fidldd,Dims=dimsldd,Pos=posldd,/Band_Only if fidldd eq -1 then return Base=Widget_Auto_Base(Title='Morphometric Attributes') Base1=Widget_Base(Base,/Column,/Frame) OutFM=Widget_OutFM(Base1,Uvalue='outf',Func='Envi_Out_Check',/Auto) result=Auto_Wid_Mng(Base) if result.accept eq 0 then return z5=float(Envi_Get_Data(fid=fid,dims=dims,pos=pos)) temp=envi_get_projection(fid=fid,units=units,pixel_size=pixel_size) pixel_size=float(pixel_size) temp=strupcase(envi_translate_projection_units(units)) if (temp eq 'DEGREES') or (temp eq 'SECONDS') then begin if temp eq 'SECONDS' then pixel_size/=3600d envi_convert_file_coordinates, fid,lindgen(ns,nl) mod ns,lindgen(ns,nl)/ns,xps,yps,/to_map xps=pixel_size[0]*111320.7d *cos(abs(yps)/!radeg) yps=replicate(pixel_size[1]*111320.7d,ns,nl) endif else begin xps=replicate(pixel_size[0],ns,nl) yps=replicate(pixel_size[1],ns,nl) endelse z1=shift(z5,1,1) z2=shift(z5,0,1) z3=shift(z5,-1,1) z4=shift(z5,1,0) z6=shift(z5,-1,0) z7=shift(z5,1,-1) z8=shift(z5,0,-1) z9=shift(z5,-1,-1) p=(z3+z6+z9-z1-z4-z7)/(6*xps) q=(z1+z2+z3-z7-z8-z9)/(6*yps) r=((z1+z3+z4+z6+z7+z9)-2*(z2+z5+z8))/(3*xps*xps) t=((z1+z2+z3+z7+z8+z9)-2*(z4+z5+z6))/(3*yps*yps) s=(z3+z7-z1-z9)/(4*xps*yps) slope=sqrt(p^2+q^2)*100 zeroed=where(slope eq 0) res=temporary(slope) bnames='slope' ;% gradfactor=sqrt(p^2+q^2)/sqrt(1+p^2+q^2) res=[[[res]],[[100*temporary(gradfactor)]]] bnames=[bnames,'gradfactor'] ;m/m (*100) ; aspect=(atan(p,q)+!pi)*180/!pi ;-90*(1-sign(q))*(1-abs(sign(p)))+180*(1+sign(p))-180*sign(p)*acos(-q/sqrt(p^2+q^2))/!pi ;graus ; res=[[[res]],[[temporary(aspect)]]] ; bnames=[bnames,'aspect'] ;degrees horcurv=-(r*q^2-2*p*q*s+t*p^2)/((p^2+q^2)*sqrt(1+p^2+q^2)) if zeroed[0] ne -1 then horcurv[zeroed]=0 res=[[[res]],[[1000*temporary(horcurv)]]] bnames=[bnames,'horcurv'] ;1/m (*1000) plancurv=-(r*q^2-2*p*q*s+t*p^2)/((p^2+q^2)^1.5) if zeroed[0] ne -1 then plancurv[zeroed]=0 res=[[[res]],[[1000*temporary(plancurv)]]] bnames=[bnames,'plancurv'] ;1/m (*1000) vertcurv=-(r*p^2+2*p*q*s+t*q^2)/((p^2+q^2)*((1+p^2+q^2)^1.5)) if zeroed[0] ne -1 then vertcurv[zeroed]=0 res=[[[res]],[[1000*temporary(vertcurv)]]] bnames=[bnames,'vertcurv'] ;1/m (*1000) meancurv=-((1+q^2)*r-2*p*q*s+(1+p^2)*t)/(2*(1+p^2+1^2)^1.5) if zeroed[0] ne -1 then meancurv[zeroed]=0 res=[[[res]],[[1000*meancurv]]] bnames=[bnames,'meancurv'] ;1/m (*1000) unspher=sqrt(((r*sqrt((1+q^2)/(1+p^2))-t/sqrt((1+q^2)/(1+p^2)))^2)/(1+p^2+q^2)+(p*q*r*sqrt((1+q^2)/(1-p^2))-2*s*sqrt((1+q^2)/(1-p^2))+p*q*t/sqrt((1+q^2)/(1+p^2)))^2)/(2*(1+p^2+q^2)^1.5) ;1/m if zeroed[0] ne -1 then unspher[zeroed]=0 res=[[[res]],[[1000*unspher]]] bnames=[bnames,'unsphericity'] ;1/m (*1000) difcurv=(r*q^2-2*p*q*s+t*p^2)/((p^2+q^2)*sqrt(1+p^2+q^2))-((1+q^2)*r-2*p*q*s+(1+p^2)*t)/(2*(1+p^2+q^2)^1.5) if zeroed[0] ne -1 then difcurv[zeroed]=0 res=[[[res]],[[1000*difcurv]]] bnames=[bnames,'difcurv'] ;1/m (*1000) rotor=((p^2-q^2)*s-p*q*(r-t))/((p^2+q^2)^1.5) if zeroed[0] ne -1 then rotor[zeroed]=0 res=[[[res]],[[1000*temporary(rotor)]]] bnames=[bnames,'rotor'] ;1/m (*1000) khe=unspher-difcurv if zeroed[0] ne -1 then khe[zeroed]=0 res=[[[res]],[[1000*temporary(khe)]]] bnames=[bnames,'khe'] ;1/m (*1000) kve=unspher+difcurv if zeroed[0] ne -1 then kve[zeroed]=0 res=[[[res]],[[1000*temporary(kve)]]] bnames=[bnames,'kve'] ;1/m (*1000) kmin=meancurv-unspher if zeroed[0] ne -1 then kmin[zeroed]=0 res=[[[res]],[[1000*temporary(kmin)]]] bnames=[bnames,'kmin'] ;1/m (*1000) kmax=meancurv+temporary(unspher) if zeroed[0] ne -1 then kmax[zeroed]=0 res=[[[res]],[[1000*temporary(kmax)]]] bnames=[bnames,'kmax'] ;1/m (*1000) totGausscurv=(r*t-s^2)/((1+p^2+q^2)^2) if zeroed[0] ne -1 then totGausscurv[zeroed]=0 res=[[[res]],[[(1e6)*temporary(totGausscurv)]]] bnames=[bnames,'totGausscurv'] ;1/m2 (*1000000) totringcurv=((p^2+q^2)*s-p*q*(r-t)^2)/(((p^2+q^2)*(1+p^2+q^2))^2) if zeroed[0] ne -1 then totringcurv[zeroed]=0 res=[[[res]],[[1000*temporary(totringcurv)]]] bnames=[bnames,'totringcurv'] ;1/m2 (*1000) totaccumcurv=temporary(meancurv)^2-temporary(difcurv)^2 if zeroed[0] ne -1 then totaccumcurv[zeroed]=0 res=[[[res]],[[(1e6)*temporary(totaccumcurv)]]] bnames=[bnames,'totaccumcurv'] ;1/m2 (*1000000) longconv=-(r*p^2+2*p*q*s+t*q^2)/(p^2+q^2) if zeroed[0] ne -1 then longconv[zeroed]=0 res=[[[res]],[[1000*temporary(longconv)]]] bnames=[bnames,'longcurv'] ;1/m (*1000) crosssecconv=-(t*p^2-2*p*q*s+r*q^2)/(p^2+q^2) if zeroed[0] ne -1 then crosssecconv[zeroed]=0 res=[[[res]],[[1000*temporary(crosssecconv)]]] bnames=[bnames,'crossseccurv'] ;1/m (*1000) ;calculating contributing area areacontr=xps*yps viz=lindgen(ns,nl)+([0,ns+1,1,-ns+1,ns,-ns,ns-1,-1,-ns-1])[Envi_Get_Data(fid=fidldd,dims=dimsldd,pos=posldd)] n=histogram(viz,min=0l,max=n_elements(viz)-1) for i=0l,ns*nl-1 do begin if n[i] ne 0 then continue areacontr[viz[i]]+=areacontr[i] ii=viz[i] while (n[ii] eq 1) and (ii ne viz[ii]) do begin areacontr[viz[ii]]+=areacontr[ii] ii=viz[ii] endwhile if n[ii] gt 1 then n[ii]-- endfor res=[[[res]],[[areacontr]]] bnames=[bnames,'areacontr'] ;m2 ;calculating level difference from the top dZ=Envi_Get_Data(fid=fid,dims=dims,pos=pos) temp=sort(temporary(areacontr)) for i=0l,n_elements(temp)-1 do dZ[viz[temp[i]]]>=dZ[temp[i]] dZ=(temporary(dZ)-Envi_Get_Data(fid=fid,dims=dims,pos=pos)) > 0 res=[[[res]],[[temporary(dZ)]]] bnames=[bnames,'desvTopo'] ;m ;calculating level difference in flow direction temp=Envi_Get_Data(fid=fid,dims=dims,pos=pos) temp=temp[viz]-((temp+temp[viz[viz]])/2) dver=fltarr(ns,nl) for i=0l,n_elements(temp)-1 do if abs(dver[viz[i]]) lt abs(temp[i]) then dver[viz[i]]=temp[i] temp=0b res=[[[res]],[[temporary(dver)]]] bnames=[bnames,'desvVert'] ;m ;calculating slope in flow direction declD8=float(Envi_Get_Data(fid=fid,dims=dims,pos=pos)) declD8=100*((declD8-declD8[viz])>0)/(sqrt((xps+xps[viz])^2+(yps+yps[viz])^2)/2) temp=where(finite(declD8) eq 0) if temp[0] ne -1 then declD8[temporary(temp)]=0 res=[[[res]],[[temporary(declD8)]]] bnames=[bnames,'declD8'] ;% ;calculating maximum Strahler order ordStrMax=replicate(1b,ns,nl) idren=lindgen(ns,nl) k=2l while 1 do begin idren=where(histogram(viz[idren],min=0l) gt 1) if idren[0] eq -1 then break ordStrMax[idren]=k while 1 do begin temp=where((ordStrMax[idren] eq k) and (ordStrMax[viz[idren]] eq (k-1))) if temp[0] eq -1 then break idren=viz[idren[temp]] ordStrMax[idren]++ endwhile idren=where(ordStrMax eq k) k++ endwhile ildd0=where(Envi_Get_Data(fid=fidldd,dims=dimsldd,pos=posldd) eq 0) if ildd0[0] ne -1 then begin for i=0l,n_elements(ildd0)-1 do begin temp=where(viz eq ildd0[i]) temp=temp[where(temp ne ildd0[i])] temp=temp[where(ordStrMax[temp] eq max(ordStrMax[temp]))] ordStrMax[ildd0[i]]=ordStrMax[temp]+((n_elements(temp) gt 1) and (ordStrMax[temp] ne 0)) endfor endif res=[[[res]],[[temporary(ordStrMax)]]] bnames=[bnames,'ordStrMax'] ;adimensional res[*,0,*]=0 res[0,*,*]=0 res[*,nl-1,*]=0 res[ns-1,*,*]=0 if result.outf.in_memory eq 1 then begin Envi_Enter_Data,res,r_fid=fid,bnames=bnames,Map_Info=Envi_Get_Map_info(fid=fid) endif else begin Openw,fidc,result.outf.name,/Get_Lun Writeu,fidc,res Free_Lun,fidc Envi_Setup_Head,Fname=result.outf.name,ns=ns,nl=nl,nb=n_elements(bnames),$ Data_Type=size(res,/Type),bnames=bnames,/Write,/Open,$ interleave=0,Map_Info=Envi_Get_Map_info(fid=fid) endelse print,'end' END