以前就想寫,但是因為envi可以就一直沒弄。今天正好有個機會,就做了這個事情。ENVI中在主窗口中pixel locator可以實現,但是當我們需要讀入很多的數據的時候,也就是批量處理的時候,顯然編程來的更快。這里只是寫了單獨輸入參數的pro,批處理的是偶可以再寫一個pro,讀入坐標到數組,利用for循環調用就行了。
來說一下思路:
首先,我們很容易通過envi提供的一些函數獲取影像的基本信息,包括dn值的二維數組,坐標信息,像元大小,以及左上角坐標。(envi_open_file,envi_file_query,envi_get_data,envi_get_projection,這四個函數會經常使用。)
其次,根據輸入的經緯度,利用坐標轉換函數將經緯度轉換為圖像對應的坐標。
接下來,利用和左上角的坐標差值,除上對應的xsize或者ysize就得到了行列號,sample和line。
最后,用sample和line作為索引,從獲取的二維數組中讀取dn值即可。
注意:經緯度中,緯度對應y,經度對應x,x的坐標差除上xsize得到的是列,y的坐標差除以ysize得到的是行。
而且,idl中,envi也是,數組的訪問是[column,row]的形式取的,所以最后是data[sample,line]得到的就是正確的。如果不確定的話,可以和envi中的pixel locator進行對比。
附上代碼:
;+
; :Author:Zhigang
; :Copyright:Niglas
; Email:zhigang_niglas@163.com
;-
PRO locateDN_HJ,latitude,longtude
;
;LOAD FUNCTIONS' MODULES OF ENVI
COMPILE_OPT IDL2
ENVI,/RESTORE_BASE_SAVE_FILES
ENVI_BATCH_INIT
;define the path
imagePath = 'E:\temp\HJ1B-CCD1-451-76-20130926-L20001058628-1.TIF'
;open HJ image
ENVI_OPEN_FILE,imagePath,r_fid = fid
ENVI_FILE_QUERY, fid, dims=dims;some parameters will be used to get data
;get the projection information
image_proj = ENVI_GET_PROJECTION(fid = fid)
;create a geographic projection, can express the latitude and longtitude
geo_proj = ENVI_PROJ_CREATE(/geo)
;convert input lat and long to coordinate under image projection
;NOTED:In the WGS-84, X is longtude, Y is latitude.
ENVI_CONVERT_PROJECTION_COORDINATES,longtude,latitude,geo_proj,image_x,image_y,image_proj
;read metadata from image
mapinfo=ENVI_GET_MAP_INFO(fid=fid)
;help,mapinfo;query the mapinfo structure, understand the MC is corner coordinate,PS is pixel Size
; print,mapinfo.MC
; print,mapinfo.PS
;
;Geolocation of UpperLeft
;
ULlat=mapinfo.MC(3);Y is latitude
ULlon=mapinfo.MC(2);X is longtude
;2. Pixel Size
Xsize=mapinfo.PS(0)
Ysize=mapinfo.PS(1)
;calculate the row and column according to x,y
sample = FIX(ABS((image_x- ULlon)/Xsize));abs is determin the positive value, fix is get integer number
line = FIX(ABS((image_y - ULlat)/Ysize))
;print,thisRow,thisColumn
;get data via row and column
DN_data= ENVI_GET_DATA(fid = fid,dims = dims,pos = 0)
;help,DN_data
;get_data
dn = DN_data[sample,line]
;write to file
PRINT,dn
;Exit ENVI
ENVI_BATCH_EXIT
END