- 积分
- 3064
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-6-20
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2013-1-3 11:42:56
|
显示全部楼层
发现网上有inside函数,http://www.astro.washington.edu/ ... brary28.html?INSIDE
Viewing contents of file '../idllib/contrib/esrg_ucsb/inside.pro'
function inside,xx,yy,xpoly,ypoly,psym=psym,nbin=nbin
;+
; ROUTINE: inside
;
; USEAGE: RESULT = inside(xx,yy)
; RESULT = inside(xx,yy,xpoly,ypoly,nbin=nbin)
;
; PURPOSE: find the indices of the data quantities xx and yy which
; lie inside a given polygonal region
;
; INPUT:
; xx vector of quantity x, ploted on x axis of scatter plot
; yy vector of quantity y, ploted on y axis of scatter plot
;
; OPTIONAL INPUT:
;
; xpoly vector of x coordinates vertices of polygon
; ypoly vector of y coordinates vertices of polygon
;
; OPTIONAL KEYWORD INPUT:
;
; psym if set, selected points are marked with symbol psym.
; if psym=3, the screen is erased and a new plot is drawn.
; this only works if the boundaries are specified
; interactively.
;
; nbin number of bins with which to resolve the polygon.
; A two element vector may be used to specify number
; in x and y directions seperately. (DEFAULT = 100)
;
; OUTPUT:
; returned value = vector of array indecies inside polygon
;
; PROCEDURE: If xpoly and ypoly are not supplied, the operator
; specifies the region using the interactive graphics
; procedure, TRACE.
;
;
; EXAMPLE
; x=randomn(seed,20,20)
; y=sin(5*x)+.8*randomn(seed,20,20)
; plot,x,y,psym=3
; i1=inside(x,y,psym=4) ; the indecies of the data points
; i2=inside(x,y,psym=5) ; which are inside polygonal area
;
;
; AUTHOR: Paul Ricchiazzi oct92
; Earth Space Research Group, UCSB
;-
if n_params() eq 0 then begin
xhelp,'inside'
return,0
endif
if keyword_set(xpoly) eq 0 then begin
trace,xpoly,ypoly,/restore
frmx= '("xpoly= [",g8.3,'+string(n_elements(xpoly)-1)+'(",",g8.3),"]")'
frmy= '("ypoly= [",g8.3,'+string(n_elements(xpoly)-1)+'(",",g8.3),"]")'
print,strcompress(string(form=frmx,xpoly),/remove_all)
print,strcompress(string(form=frmy,ypoly),/remove_all)
endif
xmin=min(xpoly)
xmax=max(xpoly)
ymin=min(ypoly)
ymax=max(ypoly)
xdif=xmax-xmin
ydif=ymax-ymin
xmin=xmin-.1*xdif
xmax=xmax+.1*xdif
ymin=ymin-.1*ydif
ymax=ymax+.1*ydif
xdif=xmax-xmin
ydif=ymax-ymin
;
if keyword_set(nbin) eq 0 then begin
nx=100
ny=100
endif else begin
if n_elements(nbin) eq 1 then begin
nx=nbin
ny=nbin
endif else begin
nx=nbin(0)
ny=nbin(1)
endelse
endelse
;
xxx=fix((nx-1)*((xx-xmin)/xdif > 0. < 1.))
yyy=fix((ny-1)*((yy-ymin)/ydif > 0. < 1.))
xb=(nx-1)*(xpoly-xmin)/xdif > 0. < (nx-1)
yb=(nx-1)*(ypoly-ymin)/ydif > 0. < (ny-1)
ii=polyfillv(xb,yb,nx,ny)
mm=bytarr(nx,ny)
mm(ii)=1
jj=where(mm(xxx,yyy) eq 1)
if keyword_set(psym) ne 0 then oplot,xx(jj),yy(jj),psym=psym
return,jj
end
|
|