게시판 답변이 생성됐습니다
-
글
-
img = read_tiff(filename, geotiff=geoinfo)
위와 같이 읽으셨다면 img라는 내용말고도 변수감시창에 GEOTAG라는 내용이 확인되실 것입니다. 여기에 지리정보 및 기타 정보가 포함되어 있습니다. 이를 이용해주시면 됩니다.
GEOTIFF 포맷의 Landsat8 데이터를 읽어보니 아래와 같은 Geotag가 확인이 되는데요, img 배열을 원하시는 수식에 맞게 처리하시어 이를 new_img라는 배열로 저장하셨다면
g_tags = { $
ModelPixelScaleTag: [30, 30, 0d], $
ModelTiepointTag: [0., 0., 0., 175500., 4269300., 0], $
GTModelTypeGeoKey: 1, $ ; (ModelTypeProjected)
GTRasterTypeGeoKey: 2, $ ; (RasterPixelIsArea)
ProjectedCSTypeGeoKey: 32652, $ ; (user-defined)
Geogangularunitsgeokey : 9102, $
Geoglinearunitsgeokey : 9001, $
Gtcitationgeokey : ‘UTM zone 52 N with WGS84′ $
}
WRITE_TIFF,’calculated_data.tif’, new_img, geotiff=g_tags
와 같이 처리하시면 됩니다.
물론 geotag 내용은 파일에 따라 다를테니 확인하시어 변경해주시면 됩니다.
도움말의 WRITE_TIFF 내용을 확인해보시면 정확할 것 같습니다.ENVI에서 roi 설정 및 이에 대한 파일을 저장해놓으셨다면 거의 다 하신 것입니다 ㅎㅎ
다만 예측하셨던 envi_get_roi 계열이 아니라 envi_restore_rois를 이용하시면 됩니다.
이후 도움말의 envi_roi_to_image_doit 내용을 보시면 example에 잘 나와 있습니다. 역시 example 내부에서도 envi_restore_rois 프로시저를 사용한 것이 보이네요 ㅎㅎ한장 강의 속 링크로 넘어가지 않는 이유는 아래의 설명을 확인하시면 도움이 될 것 같습니다.
http://idl.selab.re.kr/?p=1766
그리고 현재 GRIB에 대한 라이브러리는 아래 링크에서 다운로드 하실 수 있습니다.
http://www.exelisvis.com/docs/GRIB_Routines.html혹은 첫번째 링크의 안내처럼 ExelisVIS.com 에 회원가입 하신 후 extensions 쪽에서 다운로드 하셔도 됩니다.
usersym은 plot 프로시저 구동 시 psym 키워드에 해당하는 심볼에 대해 사용자가 직접 만드는 프로시저입니다.
단순히 원하는 지점을 포인트로 표현하고 싶으시다면 plot 프로시저 사용 시 psym=3 이라는 키워드를 추가해주시면 됩니다.for k=0, lon_num*lat_num-1 do begin
value = bias_matrix_2D[0,k]
하신 것과 정반대로 하시면 될 것 같습니다.
실제데이터 내부의 특정 행 혹은 열을 value, x, y 등과 같이 선언하신 것 같은데요, 비어있는 배열을 미리 만들어놓은 후
temp = fltarr(XX, YY)
temp[*, i] = value[1, *]
와 같이 쓰시면 됩니다.IDL 설치 시의 기본 shp 파일로는 정확한 마스킹이 안된다는 의미신지요?
그렇다면 해상도가 좋은 shp 파일을 다운로드하셔서 사용하시면 될 것 같습니다.
NOAA에서 제공하는 GSHHS 데이터를 추천합니다(http://www.ngdc.noaa.gov/mgg/shorelines/data/gshhg/latest/)
원하시는 정밀도가 어느정도인지는 모르겠습니다만, GSHHS의 i(Intermediate)가 IDL에서 hires(high-resolution)키워드를 사용했을 때의 정밀도와 유사합니다. 이보다 정밀도가 높은 h(high), f(full) 데이터를 이용하시면 해결이 되지 않을까 싶습니다.작성한지 조금 오래됐지만 아래의 한장 강의 내용을 확인해주십시오. 내용은 유효합니다.
http://idl-envi.co.kr/callcenter/tips/file/3_1241073553.pdf우선 1번째 질문은 해결이 된 것으로 판단합니다.
2번째 질문을 보고 정말 없나? 하고 소스코드를 보니 정말 없네요 ㅎㅎ
하지만 다행히도 mpfit에는 yfit이라는 키워드가 있습니다. yfit 키워드를 통해 획득되는 값은 y = ax+b에서 x에 정확히 1:1 매칭이 되는 y 값이 산출됩니다(정확히 표현하자면 1:1 매칭보다는 error값이 최소가 되는 방향으로의). 이렇게 획득된 값을 실제 y값과의 variance 비로 표현해주면 됩니다.
R2 = variance(yfit)/variance(y)
또는 regress 함수의 correlation 키워드를 사용하여 R값(상관계수)을 산출 후 제곱을 해주셔도 됩니다. R은 상관계수, R2는 결정계수를 나타냅니다.A[wh] = !values.f_nan
에서 바로 out of range가 나오나요? 이상하네요, 문법적으로는 문제가 없고, 저는 잘 작동하는데요 ㅎㅎ
다만, 원본처럼 A가 정수형일 때는 !values.f_nan처리를 해도 실제 NAN이 아닌 0으로 바뀌네요, 이전에 A를 미리 실수형으로 변경해야 정상적으로 NAN으로 변경됩니다.-999의 개수를 아시려면 where 함수 내 키워드 하나만 추가해주시면 됩니다.
wh = where(A eq -999, ct)
라고 하면 ct라는 변수는 A eq -999를 만족하는 개수값을 갖게 됩니다.우선 where를 제외하셔야 할 것 같습니다.
위 내용대로 실행한다면 where를 사용 시 결과는 0, 1, 2로만 표현이 됩니다.
즉, where(A[0:2] ne -999] 이나, where(A[3:5] ne -999] 나, where(A[6:8] ne -999]의 결과는 각각의 위치인 0, 1, 2, 3, 4, 5 가 아닌 0, 1, 2, 0, 1, 2, 1, 2와 같이 나옵니다. 매 for 문 구동시 위치가 다시 0으로 초기화되기 때문입니다.-999를 제외하기 위해 where와 ne을 설정하신거라 생각하는데요, for문 이전에 다음과 같이 where를 사용하시는게 더 좋을 것 같습니다.
wh = where(A eq -999)
A[wh] = !values.f_nan
입니다. 위와 같이 처리하게 되면 A 내부의 -999는 NAN(not a number)가 되고, 이에 대한 평균값 산출 시 nan 키워드를 사용해주시면 됩니다.
mean(A[3*i:3*i+2], /nan)단, 최초 배열이 정수형이기 때문에 정확한 값을 얻으시려면 실수형으로 바꿔주는 과정도 필요해보입니다.
2번째 축 생성 시 최초 축과의 비교? 값이 들어가게 됩니다. 키워드는 COORD_TRANSFORM 입니다.
IDL 도움말에서 확인하실 수 있습니다.
위 경우 yaxis=axis(‘Y’,location=’right’,title=’Accumulated rainfall [mm]‘, axis_range=[0,100], COORD_TRANSFORM=[0, 100./60]) 정도가 될 것 같습니다.IDL 버전 얼마를 사용 중이신가요?
8.3을 사용하신다면 가장 간단한 방법은 ORDERDHASH를 사용하는 것입니다.
header = ORDEREDHASH(‘a1′, findgen(5), ‘a2′, ‘This is header’, ‘a3′, 0.999, ‘a4′, 999L)
입니다.1. 기본적으로 title의 위치 조절 기능은 없는 것으로 알고 있습니다.
그렇기 때문에 window 생성 시 사이즈 조절과 map_set 실행 시 position 키워드를 이용하는 방법이 필요해보입니다. 그게 아니라면 (꼼수로) 제목을 넣지 않고, 나중에 xyouts 프로시저를 이용하는 방법도 괜찮을 것 같습니다.2. 질문을 정확히 이해하지 못했습니다. limit 키워드를 설정하셨다면 맵 안에 표현이 될텐데요, 밖에 생성되는 점은 어떤 함수나 프로시저를 이용했을 때인가요? 혹시 괜찮으시다면 해당 부분의 코드를 보내주시면 확인해보겠습니다.
IDL 8.3에서 확인을 해보니, mapcontinents함수의 /hires 키워드는 잘 작동합니다.
new graphics라고 하셨는데요, 그러면 map_set 프로시저 대신 map 함수를 통한 설정을 하셨는지 궁금합니다.
아래 예문처럼 실행했을 때 두 결과가 다르게 나타나는지 확인해주세요. 키워드 사용 후에도 결과가 같다면 다시 말씀해주셔요.IDL> maptest1 = map(‘geographic’, limit = [30,120,45,135])
IDL> mc1 = mapcontinents()
IDL> maptest1 = map(‘geographic’, limit = [30,120,45,135])
IDL> mc2 = mapcontinents(/hires)
