- This topic has 2개 답변, 2명 참여, and was last updated 7 years, 1 month 전에 by
.
-
글
-
-
안녕하세요,
바람장미(풍향에 따른 풍속의 빈도)를 idl로 표출하는 도중에 질문이 있어서 글 올립니다.
코드는 cgWindRose 코요테 프로그램을 이용하였는데요,
이 프로그램을 간단하게 설명하면, 풍속을 일정 간격으로(5 m/s, 10 m/s 등) 나누고 풍향을 여러 방위로(8방위, 16방위 등) 나누어서
입력된 풍향, 풍속 데이터가 각 방위 그리고 각 풍속별로 몇 퍼센트인가를 계산하여 그림으로 표출하는 프로그램입니다.
그런데 해당 프로그램에서 사용되는 펑션인 Hist_nd는 풍속 구간을 일정하게 나누도록 계산됩니다.
예를들어 5m/s씩 나눈다고 하면 0-5, 5-10 , 10-15, 15-20, 20 m/s 이상, 이런식으로 구간이 생성되는데
저는 0-0.3, 0.3-3.3, 3.3-7.9, 7.9-13.8, 13.8 이상
이렇게 일정하지 않은 5개의 구간을 설정하여 풍향별 빈도를 알아내야하는 상황입니다.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185;+; NAME:; HIST_ND;; PURPOSE:;; Perform an N-dimensional histogram, also known as the joint; density function of N variables, ala HIST_2D.;; CALLING SEQUENCE:; hist=HIST_ND(V,[BINSIZE,MIN=,MAX=,NBINS=,REVERSE_INDICES=]);; INPUTS:;; V: A NxP array representing P data points in N dimensions.;; BINSIZE: The size of the bin to use. Either an N point vector; specifying a separate size for each dimension, or a scalar,; which will be used for all dimensions. If BINSIZE is not; passed, NBINS must be.;; OPTIONAL INPUTS:;; MIN: The minimum value for the histogram. Either a P point; vector specifying a separate minimum for each dimension, or; a scalar, which will be used for all dimensions. If; omitted, the natural minimum within the dataset will be; used.;; MAX: The maximum value for the histogram. Either a P point; vector specifying a separate maximmum for each dimension, or; a scalar, which will be used for all dimensions. If omitted,; the natural maximum within the dataset will be used.;; NBINS: Rather than specifying the binsize, you can pass NBINS,; the number of bins in each dimension, which can be a P point; vector, or a scalar. If BINSIZE it also passed, NBINS will; be ignored, otherwise BINSIZE will then be calculated as; binsize=(max-min)/nbins. Note that *unlike* RSI's version; of histogram as of IDL 5.4, this keyword actually works as; advertised, giving you NBINS bins over the range min to max.;; KEYWORD PARAMETERS:;; MIN,MAX,NBINS: See above;; REVERSE_INDICES: Set to a named variable to receive the; reverse indices, for mapping which points occurred in a; given bin. Note that this is a 1-dimensional reverse index; vector (see HISTOGRAM). E.g., to find the indices of points; which fell in a histogram bin [i,j,k], look up:;; ind=[i+nx*(j+ny*k)]; ri[ri[ind]:ri[ind+1]-1];; See also ARRAY_INDICES for converting in the other; direction.;; OUTPUTS:;; hist: The N-Dimensional histogram, an array of size; N1xN2xN3x...xND where the Ni's are the number of bins; implied by the data, and/or the optional inputs min, max and; binsize.;; OPTIONAL OUTPUTS:;; The reverse indices.;; EXAMPLE:;; v=randomu(sd,3,100); h=hist_nd(v,.25,MIN=0,MAX=1,REVERSE_INDICES=ri);; SEE ALSO:;; HISTOGRAM, HIST_2D;; MODIFICATION HISTORY:;; Wed Nov 3 17:40:21 2010, J.D. Smith;; Handle 1D input with out of range elements.;;; Mon Mar 5 09:45:53 2007, J.D. Smith <jdsmith@as.arizona.edu>;; Correctly trim out of range elements from the; histogram, when MIN/MAX are specified. Requires IDL; v6.1 or later.;; Tue Aug 19 09:13:43 2003, J.D. Smith <jdsmith@as.arizona.edu>;; Slight update to BINSIZE logic to provide consistency; with HIST_2D.;; Fri Oct 11 10:10:01 2002, J.D. Smith <jdsmith@as.arizona.edu>;; Updated to use new DIMENSION keyword to MAX/MIN.;; Fri Apr 20 12:57:34 2001, JD Smith <jdsmith@astro.cornell.edu>;; Slight update to NBINS logic. More aggressive keyword; checking.;; Wed Mar 28 19:41:10 2001, JD Smith <jdsmith@astro.cornell.edu>;; Written, based on HIST_2D, and suggestions of CM.;;-;##############################################################################;; LICENSE;; Copyright (C) 2001-2010 J.D. Smith;; This file is free software; you can redistribute it and/or modify; it under the terms of the GNU General Public License as published; by the Free Software Foundation; either version 2, or (at your; option) any later version.;; This file is distributed in the hope that it will be useful, but; WITHOUT ANY WARRANTY; without even the implied warranty of; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU; General Public License for more details.;; You should have received a copy of the GNU General Public License; along with this file; see the file COPYING. If not, write to the; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,; Boston, MA 02110-1301, USA.;;##############################################################################function hist_nd,V,bs,MIN=mn,MAX=mx,NBINS=nbins,REVERSE_INDICES=ris=size(V,/DIMENSIONS)if n_elements(s) ne 2 then message,'Input must be N (dimensions) x P (points)'if s[0] gt 8 then message, 'Only up to 8 dimensions allowed'imx=max(V,DIMENSION=2,MIN=imn)if n_elements(mx) eq 0 then mx=imxif n_elements(mn) eq 0 then mn=imnif s[0] gt 1 then beginif n_elements(mn) eq 1 then mn=replicate(mn,s[0])if n_elements(mx) eq 1 then mx=replicate(mx,s[0])if n_elements(bs) eq 1 then bs=replicate(bs,s[0])if n_elements(nbins) eq 1 then nbins=replicate(nbins,s[0])endif else beginmn=[mn] & mx=[mx]endelseif ~array_equal(mn le mx,1b) then $message,'Min must be less than or equal to max.'if n_elements(bs) eq 0 then beginif n_elements(nbins) ne 0 then beginnbins=long(nbins) ;No fractional bins, pleasebs=float(mx-mn)/nbins ;a correct formulationendif else message,'Must pass either binsize or NBINS'endif else nbins=long((mx-mn)/bs+1)total_bins=product(nbins,/PRESERVE_TYPE) ;Total number of binsh=long((V[s[0]-1,*]-mn[s[0]-1])/bs[s[0]-1]);; The scaled indices, s[n]+N[n-1]*(s[n-1]+N[n-2]*(s[n-2]+...for i=s[0]-2,0,-1 do h=nbins[i]*temporary(h) + long((V[i,*]-mn[i])/bs[i])out_of_range=[~array_equal(mn le imn,1b),~array_equal(mx ge imx,1b)]if ~array_equal(out_of_range,0b) then beginin_range=1if out_of_range[0] then $ ;out of range lowin_range=total(V ge rebin(mn,s,/SAMP),1,/PRESERVE_TYPE) eq s[0]if out_of_range[1] then $ ;out of range highin_range AND= total(V le rebin(mx,s,/SAMP),1,/PRESERVE_TYPE) eq s[0]h=(temporary(h) + 1L)*temporary(in_range) - 1Lendifret=make_array(TYPE=3,DIMENSION=nbins,/NOZERO)if arg_present(ri) then $ret[0]=histogram(h,MIN=0L,MAX=total_bins-1L,REVERSE_INDICES=ri) $else $ret[0]=histogram(h,MIN=0L,MAX=total_bins-1L)return,retend위에 첨부한 Hist_nd 코드에서 h 부분을 수정해야 할 것 같은데요
for i=s[0]-2,0,-1 do h=nbins[i]*temporary(h) + long((V[i,*]-mn[i])/bs[i])
이 부분에 bs (binsize)가 스칼라로 고정이 되어있어서
어떤식으로 일정하지 않게 구간을 조정할 수 있는지 아이디어를 얻고싶습니다.
아이디어가 있으신분은 공유해주시면 정말 감사하겠습니다.
-
HIST_ND는 IDL이 자체적으로 갖고 있는 함수는 아니고 Coyote 라이브러리를 설치해야 쓸 수 있는 외부 기능입니다. 저는 개인적으로 HIST_ND를 사용해본 적이 없지만, 내용을 좀 확인해본 결과 IDL에 내장된 HIST_2D 함수의 N 차원 버전이라고 보면 될 것 같습니다. 그런데 IDL에서 사용되는 빈도분포 관련 함수들은 대개의 경우 그것이 내장 루틴이든 외부 루틴이든 간에 모두 등간격 bin들을 전제로 합니다. 따라서 만약 부등간격 bin에도 대응되도록 하려면 코드 자체의 세부 내용을 상당 부분 뜯어고치거나 아니면 유저가 직접 만드는 수 밖에 없는데, 사실 이것은 상당히 번거로운 작업일 것 같습니다. 그래서 아이디어만 생각해보면, 부등간격 bin들을 고려한다고 했을 때 가장 잘게 쪼갤 수 있는 최소 크기의 등간격을 가정하여 이 기준으로 결과를 얻은 후, 중간중간의 bin들을 서로 합치는 방식이 어떨까 합니다. 위에서 알려주신 간격들의 경우라면 0.1을 간격으로 하여 결과를 먼저 얻는 방식입니다. 구간이 0-0.3, 0.3-3.3, 3.3-7.9, 7.9-13.8, 13.8 이상인 경우 그냥 IDL 자체에 내장된 HISTOGRAM 함수를 사용하는 예를 들어 본다면 과정은 대략 다음과 같습니다. 데이터 값들로 구성된 배열을 data라고 가정합니다.
h = HISTOGRAM(data, MIN=0, MAX=13.8, BINSIZE=0.1)
이렇게 하면 h는 총 139개의 값들로 구성됩니다. 세부 구간들이 0~0.1, 0.1~0.2, …., 13.7~13.8, 13.8~13.9와 같은 방식으로 나눠지기 때문입니다. 그러면 이렇게 결과를 얻은 후 만약 0~0.3 구간의 결과가 필요하다면 다음과 같이 배열의 부분합을 구하면 됩니다.
h1 = TOTAL(h[0:2]) ; 0~0.3 구간
h2 = TOTAL(h[3:32]) ; 0.3~3.3 구간유사한 방식을 HIST_2D를 사용할 경우에 적용한다면 아마 다음과 같은 방식이 될 것 같습니다.
h1 = TOTAL(h[0:2, 0:2]) ; 0~0.3 구간
h2 = TOTAL(h[3:32, 3:32]) ; 0.3~3.3 구간그리고 HIST_ND의 경우라면 차원만 더 확장된다고 보면 되겠지요. 이 아이디어가 맞을지는 모르겠으나 어쨌든 이와 유사한 방식을 적용해 보시면 어떨까 합니다.
-
-
