Tuesday, October 18, 2011

2a25 scan time read

path = ''
file = '2A25.080922.61845.6.HDF'
filename = path+file
;fileout = 'dim.txt'
;openw, lun, fileout, /get_lun
hdfid = hdf_open(file)
result = hdf_vd_vdatalist(hdfid)
help, result, /structure

index = hdf_vd_find(hdfid, 'scan_time')
vdataid = hdf_vd_attach(hdfid, index)
nread = hdf_vd_read(vdataid, data, fields=scan_time)

Wednesday, March 9, 2011

cylindrical Coordinates

CV_COORD

The CV_COORD function converts 2D and 3D coordinates between the rectangular, polar, cylindrical, and spherical coordinate systems.

This routine is written in the IDL language. Its source code can be found in the file cv_coord.pro in the lib subdirectory of the IDL distribution.

Syntax

Result = CV_COORD( [, /DEGREES] [, /DOUBLE] [, FROM_CYLIN=cyl_coords | , FROM_POLAR=pol_coords | , FROM_RECT=rect_coords | , FROM_SPHERE=sph_coords] [, /TO_CYLIN | , /TO_POLAR | , /TO_RECT | , /TO_SPHERE] )

Return Value

If the value specified in the "FROM_" keyword is double precision, or if the DOUBLE keyword is set, then all calculations are performed in double precision and the returned value is double precision. Otherwise, single precision is used. If none of the "FROM_" keyword are specified, 0 is returned. If none of the "TO_" keywords are specified, the input coordinates are returned.

Arguments

This function has no required arguments. All data is passed in via keywords.

Keywords

DEGREES

If set, then the input and output coordinates are in degrees (where applicable). Otherwise, the angles are in radians.

DOUBLE

Set this keyword to force the computation to be done in double-precision arithmetic.

FROM_CYLIN

A vector of the form [angle, radius, z], or a (3, n) array of cylindrical coordinates to convert.

FROM_POLAR

A vector of the form [angle, radius], or a (2, n) array of polar coordinates to convert.

FROM_RECT

A vector of the form [x, y] or [x, y, z], or a (2, n) or (3, n) array containing rectangular coordinates to convert.

FROM_SPHERE

A vector of the form [longitude, latitude, radius], or a (3, n) array of spherical coordinates to convert.

TO_CYLIN

If set, cylindrical coordinates are returned in a vector of the form [angle, radius, z], or a (3, n) array.

TO_POLAR

If set, polar coordinates are returned in a vector of the form [angle, radius], or a (2, n) array.

TO_RECT

If set, rectangular coordinates are returned in a vector of the form [x, y] or [x, y, z], or a (2, n) or (3, n) array.

TO_SPHERE

If set, spherical coordinates are returned in a vector of the form [longitude, latitude, radius], or a (3, n) array.

Examples

Convert from spherical to cylindrical coordinates:

sph_coord = [[45.0, -60.0, 10.0], [0.0, 0.0, 0.0]]    rect_coord = CV_COORD(FROM_SPHERE=sph_coord, /TO_CYLIN, /DEGREES)   

Convert from rectangular to polar coordinates:

rect_coord = [10.0, 10.0]    polar_coord = CV_COORD(FROM_RECT=rect_coord, /TO_POLAR)   

Version History

Pre 4.0
Introduced

Tuesday, February 8, 2011

Plotting wind vector using IDL

I'm a atmospheric scientist who want to take advantage of IDL tool for visualizing observed and simulated atmosphere and precipitation data set. Until now, the ability of IDL seems to be excellent in solving most of issues I have encountered everyday. I really love this tool. However, the IDL contains very basic issues that had made me frustrated. Here, I will show several problems and some of my efforts.

Thursday, February 3, 2011

plotting cross_section using IDL

pro cross_section, input_arr, lon, lat, lonp, latp, res, $
csection = cross, xaxis = xx
;+
;NAME:
; cal_cross
;PURPOSE:
; calculate 2d cross section (radial - altitude)
; from 3d field (nx, ny, nz) and 2d lon, lat (nx, ny)
;INPUT:
; 1) two long/lat points at A and B
; 2) res ( A-X )
;NOTE
;(horizontal plain)
; A
; |\
; | X
; | \
; | X dis
; | \
; dy | X
; | \
; | X
; | \
; C_X_X_X_X_B
; dx
;
; AB(dis) : AX(res) = CB(dx) : CX (res1)

;OUTPUT
; 2 dimensional array [num, nz]
; a---------a
; | |
; | |
; | |
; A_X_X_X_X_B
; x-axis in km
;
; cross_section.pro , v 1.0 by M.-S. Park (mpark@nps.edu)
;_________________________________________________________
check_dim = size(input_arr)
if check_dim(0) eq 3 then begin
nz = check_dim(1)
nz = 31
; 1. finding lon, lat array for cross section
; (1) finding AB (dis), dx, dy in km
dlon = lonp(1) - lonp(0)
dlat = latp(1) - latp(0)
if dlon lt 0. then dlon = -dlon
if dlat lt 0. then dlat = -dlat
dx = deg2km (dlon, lonp(0), latp(0))
dy = deg2km (dlat, lonp(0), latp(0))
print, dx(0), ' km along longitude.'
print, dy(1), ' km along latitude.'
dis = (dx(0)^2. + dy(1)^2.)^(0.5)
print, dis, ' km is dis between A and B'
; (2) finding res1 in km
xx = [ 0., dx(0)]
yy = [dy(1), 0.]
res1 = res * dx(0) / dis
; (3) finding num & making cross array
num = fix(dx(0)/res1)
cross = fltarr(num, nz)
in = where(min(latp) eq latp)
; (4) res1 in km to reg1_deg.
imsi = km2deg(res1,lonp(in), latp(in))
res1_deg=imsi(0)
; (5) calculating slope
reg = linfit(lonp, latp, /double)
; (6) making lon, lat array for cross section
xx = findgen(num)*res1
lon_arr = findgen(num)*res1_deg + lonp(0)
lat_arr = fltarr(num)

for inum = 0, num - 1 do begin
lat_arr (inum) = reg(1) * lon_arr(inum) + reg(0)
endfor
; 2. finding the nearest points
ind1 = intarr(2, num) ;the first nearestpoint
ind2 = intarr(2, num) ;the second nearestpoint

for inum = 0, num - 1 do begin
pnt = Nearestpoint (lon, lat, lon_arr(inum), lat_arr(inum))
print, 'nearest index', pnt
lon_imsi = lon
lat_imsi = lat
dx = lon_arr(inum) - lon(pnt(0), pnt(1))
dy = lat_arr(inum) - lat(pnt(0), pnt(1))
; print, pnt, pnt_2
for iz = 0, nz - 1 do begin
cross(inum, iz) = input_arr (pnt(0),pnt(1),iz)
endfor
endfor
;return
endif ; some conditions for checking input
end

Thursday, April 22, 2010

TRMM 3B42 자료 읽기

fn_in = 'hdfname'
print, fn_in

sd_id=HDF_SD_START(fn_in,/read)

sds_id=HDF_SD_SELECT(sd_id,0) ; rain
HDF_SD_GETINFO,sds_id,name=n,ndims=r,type=t,natts=nats,dims=dims
rain = fltarr(dims(0), dims(1))
HDF_SD_GETDATA, sds_id, rain

sds_id=HDF_SD_SELECT(sd_id,1) ; err
HDF_SD_GETINFO,sds_id,name=n,ndims=r,type=t,natts=nats,dims=dims
err = fltarr(dims(0), dims(1))
HDF_SD_GETDATA, sds_id, err
rain1 = transpose(rain)
err1 = transpose(err)
HDF_SD_END, sd_id

Tuesday, January 19, 2010

Making IDL colorbar



질문 :

1. 청색톤의 컬러테이블 1번을 이용하여 최종 컬러테이블의 0~127번을 생성한다.
2. 적색톤의 컬러테이블 7번을 이용하여 최종 컬러테이블의 128~255번을 생성한다.
이 때, 7번 컬러테이블은 반전되어야 한다.
3. 1번과 7번 컬러테이블 모두 0번 부근의 낮은 값들은 검은색 또는 검은 색이 많이 섞여 있으므로
이 색들은 제외한다.

1) 1번 컬러테이블을 로드하여 그 컬러테이블을 구성하는 Red, Green, Blue 값을 읽어 놓습니다.
여기서는 1번의 R, G, B 값이라는 의미로 R1, G1, B1 변수를 사용합니다.

IDL> loadct, 1, /SILENT
IDL> tvlct, r1, g1, b1, /GET


2) 위 값들 중에 컬러 테이블의 왼쪽 부분인 검은 영역을 잘라내기 위하여 대략 100번부터 255번
까지만 사용합니다. 그리고 이 값들을 128개로 재구성하기 위하여 Congrid 함수를 이용합니다.
100~255번까지 156개의 색조합 배열이 있었으나 이를 128개로 재구성하는 것입니다.

IDL> r_1st=congrid(r1[100:255], 128)
IDL> g_1st=congrid(g1[100:255], 128)
IDL> b_1st=congrid(b1[100:255], 128)

128개로 정리된 1번 테이블의 R, G, B 값은 R_1st, G_1st, B_1st 변수에 넣어 두었습니다.

3) 7번 컬러테이블을 로드하여 그 컬러테이블을 구성하는 Red, Green, Blue 값을 읽어 놓습니다.
여기서는 7번의 R, G, B 값이라는 의미로 R7, G7, B7 변수를 사용합니다.

IDL> loadct, 7, /SILENT
IDL> tvlct, r7, g7, b7, /GET

4) 마찬가지로 7번 테이블 중에서도 검은색 영역을 빼고 100~255번 까지의 팔레트만을 취할 것입니다.
그리고 이 값들을 역시 128개로 재구성할 것입니다. 붉은 톤을 담당하는 테이블은 반전시키기 위하여
REVERSE 함수를 사용합니다.

IDL> r_2nd=reverse(congrid(r7[100:255], 128))
IDL> g_2nd=reverse(congrid(g7[100:255], 128))
IDL> b_2nd=reverse(congrid(b7[100:255], 128))

5) 이제부터 최종 컬러테이블에 우리가 원하는 값을 넣습니다.
R_1st, G_1st, B_1st 가 128개의 컬러 조합이므로 0~127번을 설정하게 됩니다.

IDL> tvlct, r_1st, g_1st, b_1st
6) R_2nd, G_2nd, B_2nd 가 128개의 컬러 조합이므로 128번부터 시작하면 128~255번을 설정하게 됩니다.

IDL> tvlct, r_2nd, g_2nd, b_2nd, 128
7) 다음과 같이 원하는 컬러테이블이 만들어졌음을 확인할 수 있습니다(이미 원하는 컬러테이블을 사용할
수 있는 상태입니다).

IDL> tvscl, dist(300)

REFERENCE:
http://www.idluser.org/

Saturday, January 16, 2010

read netcdf data using IDL

Using cdf2idl

http://woce.nodc.noaa.gov/woce_v3/wocedata_1/cmdac/primer/idl.htm


Reading HDF/netCDF Data using IDL
http://www.atmos.umd.edu/~gcm/usefuldocs/hdf_netcdf/IDL_hdf-netcdf.html