IDL通过经纬度定位获取DN值

以前就想写,但是因为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: empHJ1B-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

原文地址:https://www.cnblogs.com/ahnucao/p/4938776.html