- 积分
- 6918
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-1-1
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
这个程序是IDL编写的,我完全不懂这门语言。但其目的大致如下:对于一个降水场,设定一个临界值。该区域内大于该临界值的格点共同组成了一个降水主体,而在该降水主体内显然存在着不同的一个个不连续的降水个体。该程序的主要目的就是区分开这些降水个体。依次在mask矩阵中用1.2.3……标记,而非降水主体成员用0表示。希望各位大神帮帮忙,告诉小弟作者的编程思想,要是能够转成fortran,必感激不尽!!!!!!!!!代码如下::
;pro zero_in_on_maxes,grid1,grid2,thresh,critmass,ss1,ss2,ss3,ss4,entity_mask
pro isolate,grid,thresh,ie,area,iloc,jloc,maxval,mask
;IDL version of fortran subroutine ISOLAT (from McidAS NAWT command)
;Originally written by Andy Negri
;This may not be optimized at all for IDL!
;--- Searches array grid and defines entities by a threshold value.
;
;--- Input parameters are 'grid' (the array to be "contoured") and
;--- 'thresh' (the threshold value for an entity).
;
;--- Output is:
;--- ie (# of entities),
;--- area (# pixels in each entity),
;--- iloc,jloc (coord's of centroid)
;--- maxval (maximum value in each entity)
;---
;--- Additional output is the array 'mask', a mask of the original array
;--- with points less than thresh in grid set to zero in mask and
;--- points greater than or equal to thresh in grid set to the entity
;--- identifier (ie) in mask.
s=size(grid)
nx=s(1)
ny=s(2)
maxentities=800
maxentitysize=8000 ;maximum size of an entity
itogo=make_array(maxentitysize,/int,value=0)
jtogo=itogo
area=make_array(maxentities,/long,value=0)
iloc=area & jloc=area & maxval=area
isq=[-1,0,1] ;square search area
mask=make_array(nx,ny,/int,value=0)
ie=0
;--- Search for entities - loop
for jc=1,ny-2 do begin & for ic=1,nx-2 do begin ;100 loop
if mask(ic,jc) gt 0 or grid(ic,jc) lt thresh then goto, label100
;--- Brightness greater than threshold, increment cloud counter ie
ie=ie+1 ;Found one!
if ie ge maxentities then begin
print,'Too many entities - do not save any more'
ie=ie-1
goto, done
endif
;--- Starting index for fragment of cloud ie
ntogo=0 ;was originally 1
ndone=-1 ;was originally 0
;--- Place point's location in circular buffer
itogo(0)=ic ;was originally itogo(1)
jtogo(0)=jc ;same for this one
;--- Update histogram
area(ie)=area(ie)+1
iloc(ie)=iloc(ie)+ic
jloc(ie)=jloc(ie)+jc
maxval(ie)=grid(ic,jc)
;--- Update mask, marking point as found
mask(ic,jc)=ie
addonpoints:
if ndone ge ntogo then begin
iloc(ie)=round(float(iloc(ie))/float(area(ie)))
jloc(ie)=round(float(jloc(ie))/float(area(ie)))
goto, label100
endif
ndone=ndone+1
if ndone gt maxentitysize then ndone=0 ;was originally ndone=1
;--- Get center of next 8-point search
i=itogo(ndone)
j=jtogo(ndone)
;--- Search the 8 neighbors
for l=0,2 do begin & for k=0,2 do begin ;loop 50
if k eq 1 and l eq 1 then goto, label50
id=i+isq(k)
jd=j+isq(l)
if id lt 0 or id gt nx-1 then goto, label50
if jd lt 0 or jd gt ny-1 then goto, label50
if mask(id,jd) gt 0 or grid(id,jd) lt thresh then goto, label50
;--- Brightness greater than threshold so it's another point of cloud ie
ntogo=ntogo+1
if ntogo gt maxentitysize then ntogo=0 ;was originally ntogo=1
;--- Save location in circular buffer
itogo(ntogo)=id
jtogo(ntogo)=jd
;--- Update # of pixels
area(ie)=area(ie)+1
iloc(ie)=iloc(ie)+id
jloc(ie)=jloc(ie)+jd
maxval(ie)=max([maxval(ie),grid(id,jd)])
;--- Update mask, marking point as found
mask(id,jd)= ie
label50:
endfor & endfor
;--- Check that the area is large enough (> 2 pixels)
;--- If not, scratch out this entity
if area(ie) lt 3 then begin
area(ie)=0
iloc(ie)=0
jloc(ie)=0
maxval(ie)=0
for l=0,2 do begin & for k=0,2 do begin ;loop 60
id=i+isq(k)
jd=j+isq(l)
if id lt 0 or id gt nx-1 then goto, label60
if jd lt 0 or jd gt ny-1 then goto, label60
mask(id,jd)=0
label60:
endfor & endfor
ie=ie-1
goto, label100
endif
goto, addonpoints
label100:
endfor & endfor
;Clean up - shrink array sizes down to where there is data
area=area(0:ie)
iloc=iloc(0:ie)
jloc=jloc(0:ie)
maxval=maxval(0:ie)
done:
end
pro getentity,mask,entID,ss
;Get the grid box including entity 'entID'
s=size(mask) & nx=s(1) & ny=s(2)
xx=indgen(nx)#replicate(1,ny)
yy=replicate(1,nx)#indgen(ny)
ssm=where(mask eq entID,nmask)
if nmask ge 1 then begin
xstart=max([ 0,min(xx(ssm))-1])
ystart=max([ 0,min(yy(ssm))-1])
xend =min([nx-1,max(xx(ssm))+1])
yend =min([ny-1,max(yy(ssm))+1])
index=indgen(nx,ny)
ss=index(xstart:xend,ystart:yend)
endif else ss=-1
end
pro zero_in_on_maxes,grid1,grid2,thresh,critmass,ss1,ss2,ss3,ss4,entity_mask
;Given two grids of data, this procedure zeros in on all of the local maxes
;and returns the subscripts, 'ss', of the boxes containing data surrounding
;the four largest local maxes. Each box is defined by the area enclosing
;the isopleth given by 'thresh', where the value in the box can originate in
;either grid1 or grid2.
;Take the maximum value of grid1 and grid2
grid=grid1>grid2
;Isolate entities with local maxima
;nentities = number of entities isolated
;area = number of grid points in each entity
;iloc,jloc = indeces of centroid of each entity
;maxval = maximum value inside each entity
;mask = mask where each entity is marked with its identifying number,
; all other gridpoints are 0
isolate,grid,thresh,nentities,area,iloc,jloc,maxval,mask
if nentities gt 0 then begin
;Get the integrated value inside each entity (i.e., area * average value)
intvalue=fltarr(nentities)
for i=0,nentities-1 do intvalue(i)=total(grid(where(mask eq i+1)))
;Rank entities in order from large to small
order=reverse(sort(intvalue))
endif
;Get the boxes around the four largest entities
;Mimimum integrated size is 'critmass'
nentities=min([nentities,4])
ss=-1 & ss1=-1 & ss2=-1 & ss3=-1 & ss4=-1
n=0
for i=0,nentities-1 do begin
entID=order(i)+1
if intvalue(order(i)) ge critmass then begin
getentity,mask,entID,ss
if n eq 0 then ss1=ss
if n eq 1 then ss2=ss
if n eq 2 then ss3=ss
if n eq 3 then ss4=ss
n=n+1
endif
endfor
nentities=n
;New code to write entity mask to file
entity_mask=mask & entity_mask[*,*]=0
for i=0,nentities-1 do begin
entID=order(i)+1
newID=i+1
ss=where(mask eq entID)
entity_mask(ss)=newID
endfor
sz=size(entity_mask)
;Do a final check for repeated entities
if nentities eq 4 and ss4(0) eq ss3(0) then ss4=-1
if nentities ge 3 and ss3(0) eq ss2(0) then ss3=-1
if nentities ge 2 and ss2(0) eq ss1(0) then ss2=-1
done:
end
|
|