转载请注明原文地址
本文纯属交流学习，请勿用作其他用途。匿名登录ModisFTP站点可能导致侵权，所造成的一切法律后果，本人概不负责。
一 介绍
Modis免费分发，光谱通道丰富，产品体系成熟，在多个行业和领域有广泛成功的应用。已成为重要的遥感数据源之一。一般若需获取modis数据，要注册wist账号，查询订购（免费）并等待回复mail，整个流程一般约需数小时。为了避免等待，本文用IDL语言实现了modis产品的地理范围查询，返回的url直接添加到迅雷下载任务列表。
图 迅雷添加modis 下载任务
运行环境：IDL7.0以上版本，迅雷5.0版本以上
二 原理背景
- Modis官方ftp站点实时记录了5分钟产品的经纬度范围，并保存在txt文件中。Txt各列以逗号分隔，可以用记事本查看。
- Modis官方ftp站点还保存有全部常规产品hdf格式文件。
- IDL对象IDLnetUrl可以获取ftp站点目录列表，运行ftp命令，获得http协议url指向的文件。
- IDL对象IDLcomIDispatch对象可以调用com /ole对象；迅雷ThunderAgent.Agent是一种ole组件，可用于新建添加迅雷下载任务。
三 流程步骤
- 1) 指定日期、经纬度多边形polygon
- 2) 获得指定日期5分钟产品的地理范围txt文件，获得指定日期5分钟产品列表（简洁模式）
- 3) 解析txt中各5分钟产品经纬度范围、日夜模式，判断与指定的经纬度polygon是否相交；有交集则记录5分钟产品的url
- 4) 将所有符合要求的hdf文件url添加为迅雷下载任务
四 源码
IDL 源码;-----------------------------------------------------------------
pro ModisTest
QUERY_DOWNLOAD,2010,7,18,[110,110,117,117],[34,30,30,34],PRODUCTNAME='021KM'
end
;-----------------------------------------------------------------
PRO QUERY_DOWNLOAD,YEAR,MONTH,DAY,XPOLY,YPOLY,PRODUCTNAME=PRODUCTNAME
GET_MODIS_METAFILE,YEAR,MONTH,DAY,geometa1,MODLIST,/TERRA,PRODUCTNAME=PRODUCTNAME
QUERY_MODIS,geometa1,MODLIST,XPOLY,YPOLY
GET_MODIS_METAFILE,YEAR,MONTH,DAY,geometa2,MYDLIST,/AQUA,PRODUCTNAME=PRODUCTNAME
QUERY_MODIS,geometa2,MYDLIST,XPOLY,YPOLY
IF SIZE(MODLIST,/TYPE) EQ 7 THEN L = [MODLIST]
IF SIZE(MYDLIST,/TYPE) EQ 7 THEN L = [L,MYDLIST]
IF N_ELEMENTS(L) GT 0 THEN ADD_DOWNLOAD,L
END
;-----------------------------------------------------------------
PRO QUERY_MODIS,DATA,FILELIST,XPOLY,YPOLY
FLIST =FILELIST
FILELIST = -1
IF SIZE(DATA,/TYPE) NE 7 OR SIZE(FLIST,/TYPE) NE 7 OR N_ELEMENTS(XPOLY) NE N_ELEMENTS(YPOLY) THEN RETURN
G = {GEOMETA,G_ID:'',STARTTIME:'',SET:0B,ORBIT_NUMBER:1L,DAYNIGHT:0B,BOX:DINDGEN(4),X:DINDGEN(4),Y:DINDGEN(4)}
META = [G]
FOR I = 3,N_ELEMENTS(DATA)-1 DO BEGIN
IF N_ELEMENTS(DATA) LE 3 THEN RETURN;
TMP = STRSPLIT(DATA[I],",",/EXTRACT)
IF N_ELEMENTS(TMP) EQ 17 THEN BEGIN
G.G_ID = (STRSPLIT(TMP[0],'.',/EXTRACT))[2]
G.STARTTIME = TMP[1]
G.SET = TMP[2]
G.ORBIT_NUMBER = TMP[3]
G.DAYNIGHT = BYTE(TMP[4])
G.BOX = TMP[5:8]
G.X = TMP[9:12]
G.Y = TMP[13:16]
META = [META,G]
ENDIF
ENDFOR
IF N_ELEMENTS(META) GT 1 THEN META = META[1:*] ELSE RETURN
IDX = INTARR(480);
I = FIX(STRMID(FILE_BASENAME(FLIST),18,4))/5
IDX[I] = INDGEN(N_ELEMENTS(I))
ret = ['']
FOR I = 0,N_ELEMENTS(META)-1 DO BEGIN
id = IDX[FIX(META[I].G_ID)/5]
;
;换日线附近一般不满足要求
IF (MAX(META[I].X) - MIN(META[I].X)) GT 300 THEN CONTINUE
;
IF TOTAL(INSIDE(XPOLY,YPOLY,META[I].X,META[I].Y)) GT 0 $
AND META[I].DAYNIGHT LT BYTE('N') $
AND ID GT 0 THEN RET = [RET,FLIST[ID]]
ENDFOR
FILELIST = N_ELEMENTS(RET) GT 1 ? RET[1:*] : -1
END
;-----------------------------------------------------------------
PRO GET_MODIS_METAFILE,YEAR,MONTH,DAY,GEOMETA,FLIST,TERRA=TERRA,AQUA=AQUA,PRODUCTNAME=PRODUCTNAME
GEOMETA = -1
FLIST = -1;
CATCH, errorStatus
IF (errorStatus NE 0) THEN BEGIN
CATCH, /CANCEL
r = DIALOG_MESSAGE(!ERROR_STATE.msg, TITLE='URL Error', $
/ERROR)
PRINT, !ERROR_STATE.msg
RETURN
ENDIF
CASE 1 OF
KEYWORD_SET(TERRA):BEGIN
prefix = 'MOD'
SATELLITE = 'TERRA'
BREAK
END
KEYWORD_SET(AQUA):BEGIN
prefix = 'MYD'
SATELLITE = 'AQUA'
BREAK
END
ELSE:RETURN
ENDCASE
PRODUCTNAME = KEYWORD_SET(PRODUCTNAME) ? PRODUCTNAME : '021KM'
site = 'ftp://ladsweb.nascom.nasa.gov'
JDAY = JULDAY(MONTH,DAY,YEAR,0,0,0) - JULDAY(1,1,YEAR,0,0,0)
FORMAT = '("geoMeta/6/'+SATELLITE+'/",i4,"/","'+PREFIX+'03_",I4,"-",I02,"-",I02,".txt")
DIR_FMT = '("/allData/5/","'+prefix+PRODUCTNAME+'/",I4,"/",I03,"/")'
url = STRING([year,year,month,day],format=FORMAT)
DIR = site+ STRING([YEAR,JDAY],FORMAT = DIR_FMT)
; create a new url object
oUrl = OBJ_NEW('IDLnetUrl',$
CALLBACK_FUNCTION ='ddcall',$
URL_SCHEME = 'ftp',$
URL_HOST = 'ladsweb.nascom.nasa.gov',$
URL_PATH = url,$
URL_USERNAME = 'Anonymous',$
URL_PASSWORD = '',$
FTP_CONNECTION_MODE = 0)
GEOMETA = oUrl->Get( /STRING_ARRAY )
FLIST = dir + oUrl->GetFtpDirList(url=DIR,/short)
oUrl->CloseConnections
OBJ_DESTROY, oUrl
END
;-----------------------------------------------------------------
FUNCTION ddcall, status, progress, data
PRINT, status
RETURN, 1
END
;-----------------------------------------------------------------
PRO ADD_DOWNLOAD,URL
oThunder = OBJ_NEW('IDLcomIDispatch$ProgId$ThunderAgent_Agent');
IF OBJ_VALID(oThunder) THEN BEGIN
FOR I = 0,N_ELEMENTS(URL)-1 DO oThunder->ADDTASK,url[I],"","","","",1,0,5
oThunder->COMMITTASKS,1
OBJ_DESTROY ,oThunder
ENDIF
END
;-----------------------------------------------------------------
FUNCTION Inside, x, y, xpts, ypts, INDEX=index
On_Error, 1
sx = Size(xpts)
sy = Size(ypts)
IF (sx[0] EQ 1) THEN NX=sx[1] ELSE Message, 'X coordinates of polygon not a vector.'
IF (sy[0] EQ 1) THEN NY=sy[1] ELSE Message, 'Y coordinates of polygon not a vector.'
IF (NX EQ NY) THEN N = NX ELSE Message, 'Incompatable vector dimensions.'
; Close the polygon.
tmp_xpts = [xpts, xpts[0]]
tmp_ypts = [ypts, ypts[0]]
; Set up counters.
i = indgen(N)
ip = indgen(N)+1
nn = N_Elements(x)
X1 = tmp_xpts(i) # Replicate(1,nn) - Replicate(1,n) # Reform([x],nn)
Y1 = tmp_ypts(i) # Replicate(1,nn) - Replicate(1,n) # Reform([y],nn)
X2 = tmp_xpts(ip) # Replicate(1,nn) - Replicate(1,n) # Reform([x],nn)
Y2 = tmp_ypts(ip) # Replicate(1,nn) - Replicate(1,n) # Reform([y],nn)
dp = X1*X2 + Y1*Y2 ; Dot-product
cp = X1*Y2 - Y1*X2 ; Cross-product
theta = Atan(cp,dp)
ret = Replicate(0L, N_Elements(x))
i = Where(Abs(Total(theta,1)) GT 0.01, count)
IF (count GT 0) THEN ret(i)=1
; Make this a scalar value if there is only one value.
IF (N_Elements(ret) EQ 1) THEN ret=ret[0]
; If the index keyword is set, then return indices.
IF (Arg_Present(index)) THEN ret=(Indgen(/Long, N_Elements(x)))(Where(ret eq 1))
RETURN, ret
END
五 感谢
感谢 qqz的【原】根据Modis 分块计算经纬度程序 给我的启示。