위성 영상의 지정 위치 값 추출

위성 영상에서 특정 위치의 값을 읽어오고자 하는 상황입니다. ENVI Interface를 이용하고 있고, 값을 눈으로 확인하고자 하는 상황이라면, ENVI Menu의 Display > Cursor Value 를 실행하면 될 것입니다. 영상 위에 마우스가 이동하는 대로 좌표와 해당 위치의 값을 표출합니다.

Cursor Value Window

Cursor Value window

이를 ENVI Batch 프로그래밍으로 처리할 때는 어떻게 할까요? 다음 예제와 설명을 적당히 활용하면 될 것 같습니다.

1) ENVI Batch 모드를 시작합니다.

ENVI> envi, /restore_base_save_files
ENVI> envi_batch_init

2) File을 지정하고 이를 ENVI에서 Open합니다. 프로그래밍으로는 아래와 같이 하면 됩니다.

ENVI> file=’C:\mydata\myimage.img’
ENVI> envi_open_file, file, R_FID=fid        ;File 변수에 있는 파일을 읽고 관리 번호를 fid 변수에 넣습니다.

이제 ENVI 파일 포맷 대로 헤더 파일인  myimage.hdr의 정보 역시 ENVI가 관리하기 시작합니다. 위경도 좌표를 이용하여 원하는 위치를 지정한다면, a) 위경도좌표계 b) 영상파일의 좌표계 c) 영상의 픽셀 좌표계 의 단계로 원하는 위치를 찾아가야 합니다. 이 경우 위경도좌표계 투영법과 영상파일의 좌표계(TM, UTM 등) 투영법 정의가 필요합니다. 다음 단계의 예제를 보시는 게 더 쉽습니다.

3) 위경도 좌표계(GEOGRAPHIC)와 영상에 사용된 좌표계를 정의합니다.

ENVI> iProj=ENVI_PROJ_CREATE(/geographic)         ;위경도 좌표계 정의
ENVI> oProj=ENVI_GET_PROJECTION(FID=fid)

ENVI_PROJ_CREATE는 투영법과 지구타원체 등을 지정해서 좌표계를 정의합니다.  관련 글을 참고하세요(ENVI_PROJ_CREATE 함수에서 TM 좌표계 설정). 위경도 체계(GEOGRAPHIC)과 UTM 체계는 간단히 키워드로 정의할 수 있습니다. 영상의 좌표계는 영상의 헤더 파일로부터 ENVI_GET_PROJECTION 함수를 이용해 가져올 수 있습니다.

4) ENVI의 좌표변환 루틴은 단일값을 넣어 주어도 되고, 배열을 넣어주면 모든 배열 요소에 대한 좌표 변환을 수행합니다. 예를 들어, 3)의 정의를 이용하여 위도 36.6 경도 127.0 지점과 위도 35.7, 경도 128.6 지점 두개를 변환하고자 한다면,

ENVI> lat=[36.6,   35.7]
ENVI> lon=[127.0,   128.6]

ENVI> ENVI_CONVERT_PROJECTION_COORDINATES, lon, lat, iProj, XMAP, YMAP, oProj

ENVI> print, XMAP, YMAP
199928.56 344740.36
344318.39 245627.43

가장 중요한 한 줄은 ENVI_CONVERT_PROJECTIN_COORDINATES 인데요, 문법 정의가 이렇습니다.

ENVI_CONVERT_PROJECTION_COORDINATES, iXmap, iYmap, iProj, oXmap, oYmap, oProj

iProj 변수에 정의된 좌표계에서 iXmap, iYmap 좌표를 oProj 변수에 정의된 좌표계의 좌표로 변환하여 oXmap 변수와 oYmap 변수에 넣는 프로시저입니다. iProj 나 oProj에 해당하는 변수의 값, 즉 ENVI의 좌표계 정의 구조체는 3)과 같이 정의할 수 있는 거지요.

36.6N, 127.0E -> 199928.56, 344318.39 로 변환되고
35.7N, 128.6E -> 344740.36, 245627.43 으로 변환된 것을 보면 영상의 좌표계는 Korea 중부 원점 TM 좌표계입니다.

5) 영상의 Map 좌표계 좌표로부터 영상 배열의 픽셀 좌표값을 구하는 것도 비슷합니다.

ENVI> ENVI_CONVERT_FILE_COORDINATES, fid, X, Y, XMAP, YMAP
ENVI> print, X, Y
1115.2936      1260.1054
511.72030      610.41126

XMAP, YMAP 변수에 있는 영상 지도 좌표를 픽셀 좌표로 변환하여 X, Y 변수에 대입한 것입니다. 이 프로시저는 반대로 픽셀좌표로부터 영상 지도 좌표(TM, UTM 등)로 변환하는 것도 가능한데, 이 경우 /TO_MAP 키워드를 붙여주면 됩니다.

6) 이제 픽셀 좌표를 알았으므로, 영상을 읽어와 해당 픽셀좌표 값을 뽑아내면 됩니다.  

ENVI> ENVI_FILE_QUERY, fid, DIMS=dims
ENVI> data=envi_get_data(fid=fid, DIMS=dims, pos=[0])
ENVI> print, data[x, y]
0.444358 0.465551

ENVI_FILE_QUERY 프로시저를 이용하여 영상의 다양한 정보를 알아올 수 있습니다. 밴드의 개수는 몇개인지, 영상의 크기는 얼마인지(DIMS). 그리고 이러한 정보를 이용하여 영상에서 원하는 영역과 원하는 밴드를 배열로 읽어올 수 있습니다(envi_get_data). DIMS를 위와 같이 지정함으로서 공간적으로는 영상 전역을 모두 읽어오며, 밴드는 첫번째 밴드만(POS=[0]) 읽어오게 됩니다. 이제 data 변수는 2차원 배열이 됩니다(pos=[0, 1] 과 같이 2개 이상의 밴드를  지정한다면 3차원 배열이 되겠지요.

두 위치의 값을 읽어낸 결과는 0.444358과 0.46551 로 나왔습니다. 예제 영상을 NDVI 밴드로 사용했거든요.

LAT, LON 변수에 수 십개, 수 백개의 위치를 지정해도 방법은 같습니다.