一直以来处理modis产品都是用的 modis conversion toolkit(MCTK),用 IDL 来跑感觉好像也没什么问题,好像并没有去验证结果到底对不对,因为用的气溶胶数据 MOD04_L2,没法验证到底对不对,不像反射率数据可以对照地图影像。
流程也都是重投影--重采样--拼接--裁剪,其实应该也没有问题,因为最主要的步骤重投影是通过 MCTK 做的,错不了。
那么问题来了,为什么重投影是最重要的呢?
MODIS产品都是正弦投影 sinusoidal 投影坐标系,而我们一般要用的都是 WGS84 地理坐标系,不转换根本没法用,中国地区形变巨大,而且别的数据一般都不会有这种投影的。
MOD04处理完,要处理MOD02和03了,还是用 MCTK,就有坑出现了。
坑1:拼接结果偏离行政边界
因为官方说数据都是1km的,所以为了图省事就把重采样这一步给去掉了,步骤是重投影--波段合成(02+03)--拼接。拼完的结果和行政边界一对比,咋偏了这么多。后来发现重投影到 WGS84 每景影像是不同分辨率的(度为单位都在0.01左右)。
解决方法:
重投影后把每个影像重采样到0.01度,再拼接就不会偏离了。不过现在感觉这个解决方法也有点问题,可能是拼接的时候参数没设置好导致的。(看了一遍 ENVI_DOIT, 'MOSAIC_DOIT',好像没有参数控制这个,应该还是重采样的问题 )
OK,上面都弄完现在要搞 MCD19A2 的 MAIAC 数据了,我的天,MCTK 不支持这玩意,那这重投影可咋做(MCTK就是用来转投影的)。行吧,硬撸 IDL 代码,不就是搞个重投影。
坑2:自定义 geotiff 坐标信息
MODIS 产品自身不带投影信息,在每个 hdf 文件的全局属性里才有,光读数据导出 tif 是没有坐标系的。找不到好的智能解决方法,只能用蠢方法。导出一个 tif 无非就是数组加坐标信息,坐标信息 modis 的 ENVI 似乎没有(我没找到啊!),只能自己来定义。自己定义又不知道每个参数都是干啥的咋设置,没在官方帮助文档里找到,只看到列了一排看不懂的参数。
解决方法:
先根据这个 http://blog.sina.com.cn/s/blog_764b1e9d01010v69.html#cmt_5691CA7B-7F000001-13E754CBE-868-8A0 给一个hdf加上头文件,导出一个有坐标属性的tif,然后用 read_tiff 函数读取,得到 geotiff 属性,就可以复制过来自定义了。
坑3:自定义的 geotiff 起始坐标小数不见了
解决方法:
以为是哪个属性没设置对,找了好久终于找到了每个参数的含义,参考这个 http://blog.sina.com.cn/s/blog_764b1e9d01010v69.html#cmt_5691CA7B-7F000001-13E754CBE-868-8A0。发现都没问题,应该是处理过程 envi 精度丢失了。原来 envi 在处理的时候会先把小数转到 float 型,再处理,在数字后面加上D就保证用double了,成功保留了小数
到这里感觉取得了突破性的胜利,剩下只要转投影就OK了,这还不简单,有 ENVI_CONVERT_FILE_MAP_PROJECTION 工具。
坑4:转完投影又套不上行政边界了
解决方法:
原来是函数里的一个参数 warp_method (解决形变利用的方法)没设置好,默认是 0: Rotation, scaling, and translation (RST),改成 2: Triangulation 就没问题了。
坑5:每个转完的拼起来好像中间有缺失条带
解决方法:
先拼接后转投影,完美解决条带问题,不用再去插值了