LS因子计算代码（IDL)

LS因子计算代码（IDL)

;*****************************************坡度、汇水面积、流向
function downslope,DEM,cellSize,size1
slopeang=FLTARR(size1[0]+4,size1[1]+4)
celllen=FLTARR(size1[0]+4,size1[1]+4)
outflow=FLTARR(size1[0]+4,size1[1]+4)
slope=FLTARR(3,size1[0]+4,size1[1]+4)
; inflow=outflow
cols =size1[0]+4
rows=size1[1]+4
diagCellSize = cellSizesqrt(2)
; cellSize_data=[[1/(sqrt(2)*cellSize),1/cellSize,1/(sqrt(2)*cellSize)],[1/cellSize,1/cellSize,1/cellSize],[1/(sqrt(2)*cellSize),1/cellSize,1/(sqrt(2)*cellSize)]]
for i=2,size1[1]+1 do begin
for j=2,size1[0]-1 do begin
flag=0

  bestCellLen = 0.0
  bestDirection =0.0
  bestAngle = 0.0
  angle = 0.0

  c = i * cols + j

  nw = c - cols - 1
  n = c - cols
  ne1 = c - cols + 1
  w = c - 1
  e = c + 1
  sw = c + cols - 1
  s = c + cols
  se = c + cols + 1

; window1=DEM[i-1:i+1,j-1:j+1]
; if (i eq 2 ) then begin
; window1[0,]=DEM[i,j]
; endif
; if (j eq 2 ) then begin
; window1[,0]=DEM[i,j]
; endif
; if (i eq size1[0]-1 ) then begin
; window1[2,]=DEM[i,j]
; endif
; if (j eq size1[0]-1) then begin
; window1[,2]=DEM[i,j]
; endif
;
; location=where(DEM eq DEM[i,j])
; ;print,location[0],c
; slope1=float(Atan(abs(-window1+DEM[i,j])cellSize_data)!RADEG)
; ;print,slope1
if (DEM[w] ne -9999 && DEM[w] lt DEM[c])then begin
if (outFlow[w] eq !E_value) then begin

    endif else begin
      angle =(Atan((abs(DEM[c] - DEM[w]) / cellSize)))*!RADEG
      ;print,angle,'W'
      if (angle gt bestAngle) then begin
        bestAngle = angle
        bestDirection = !W_value
        bestCellLen = cellSize
        flag =1
      endif
      if (angle eq 0.0 && ~flag) then begin

        bestAngle = !MIN_SLOPE_value
        bestDirection = !W_value
        bestCellLen = cellSize
        flag = 1
      endif
    endelse
  endif
  if (DEM[nw] ne -9999 && DEM[nw] lt DEM[c])then begin
    if (outFlow[nw] eq !SE_value) then begin

    endif else begin
        angle =(Atan((abs(DEM[c] - DEM[nw]) / diagCellSize)))*!RADEG

;print,angle,‘nW’
if (angle gt bestAngle) then begin
bestAngle = angle
bestDirection = !nw_value
bestCellLen = diagCellSize
flag =1
endif
if (angle eq 0.0 && ~flag) then begin

        bestAngle =  !MIN_SLOPE_value;
        bestDirection = !nw_value
        bestCellLen = diagCellSize;
        flag = 1
      endif
    endelse
  endif
  if (DEM[n] ne -9999 && DEM[n] lt DEM[c])then begin
    if (outFlow[n] eq !S_value) then begin;防对流，

    endif else begin
      angle =(Atan(abs(DEM[c] - DEM[n]) / cellSize))*!RADEG
       ;print,angle,'n'
      if (angle gt bestAngle) then begin
        bestAngle = angle
        bestDirection = !n_value
        bestCellLen = cellSize
        flag =1
      endif
      if (angle eq 0.0 && ~flag) then begin

        bestAngle = !MIN_SLOPE_value
        bestDirection = !n_value
        bestCellLen = cellSize
        flag = 1
      endif
    endelse
  endif
  if (DEM[ne1] ne -9999 && DEM[ne1] lt DEM[c])then begin
    if (outFlow[ne1] eq !SE_value) then begin

    endif else begin
       angle =(Atan(abs(DEM[c] - DEM[ne1]) /  diagCellSize ))*!RADEG
        ;print,angle,'nW'
      if (angle gt bestAngle) then begin
        bestAngle = angle
        bestDirection = !ne_value
        bestCellLen =  diagCellSize
        flag =1
      endif
      if (angle eq 0.0 && ~flag) then begin

        bestAngle =  !MIN_SLOPE_value;
        bestDirection = !ne_value
        bestCellLen = diagCellSize;
        flag = 1
      endif
    endelse
  endif
  if (DEM[e] ne -9999 && DEM[e] lt DEM[c])then begin
    ;  if (outFlow[e] eq !S_value) then begin
    ;
    ;  endif else begin
    angle =(Atan(abs(DEM[c] - DEM[e]) / cellSize)*!RADEG)
    if (angle gt bestAngle) then begin
      bestAngle = angle
      bestDirection = !e_value
      bestCellLen = cellSize
      flag =1
    endif
    if (angle eq 0.0 && ~flag) then begin

      bestAngle = !MIN_SLOPE_value
      bestDirection = !e_value
      bestCellLen = cellSize
      flag = 1
    endif
    ;endelse
  endif
  if (DEM[se] ne -9999 && DEM[se] lt DEM[c])then begin
    ;  if (outFlow[se] eq !SE_value) then begin
    ;
    ;  endif else begin
     angle =(Atan(abs(DEM[c] - DEM[se]) /  diagCellSize )*!RADEG)
    if (angle gt bestAngle) then begin
      bestAngle = angle
      bestDirection = !Se_value
      bestCellLen =  diagCellSize
      flag =1
    endif
    if (angle eq 0.0 && ~flag) then begin

      bestAngle =  !MIN_SLOPE_value;
      bestDirection = !Se_value
      bestCellLen = diagCellSize;
      flag = 1
    endif

; ;endelse
endif
if (DEM[s] ne -9999 && DEM[s] lt DEM[c])then begin
; if (outFlow[e] eq !S_value) then begin
;
; endif else begin
angle =(Atan(abs(