一个实用的由MODIS正弦投影反推经纬度的小工具

好久没发日志了。。。
最近在没日没夜地做我的小论文，其间遇到很多问题，其中一个就是：如何得到MODIS数据网格的经纬度，而且要准确、可靠。一开始我在谷歌上搜了很久，但是很遗憾没有找到任何一个让我满意的方案，所以我就编写了如下的一个小工具：

一、该程序包含两个“工作模式”：点模式 和 文件模式
（1）点模式：
适用于单点推算经纬度，在Col、Row、X、Y处输入目标MODIS网格的列号、行号、X坐标或Y坐标。
这里的判断规则是：
如果输入了有效的行列号，则程序会无视XY坐标，以行列号作为目标网格的位置
如果想要以XY坐标作为目标网格的位置，则将行列号置空

如上图，首先在File处选择你的MODIS数据文件，设置好col和row，点击运行，则会在Output处显示推算结果，如果一个输入文件里有多个维度的子数据集，那么本程序会逐一计算其在各维度情况中的经纬度。比如，本例中，MOD09产品内有两种维度的数据，分别是1200×1200和2400×2400，对应于1公里分辨率和500米分辨率

（2）文件模式：
单点模式只能算是“题中应有之义”，很多情况下，我们需要的是大范围的网格经纬度信息，这时候“文件模式”就派上用场了。

如上图，在Input中选择一个有效的MODIS hdf文件，并设置好输出目录，点击运行。
则会在输出目录下得到 不同维度情况下的经度和纬度文件，如下：

如图，文件的命名规则是：[条带号].[X维度长度]x[Y维度长度].[经度或纬度].tif

二、验证
对于任何一个程序，验证都是最重要的。
首先，既然我们要验证，那么我们以什么为标准呢。MRT（MODIS Reprojection Tool)是官方推荐的投影转换工具，所以我觉得以MRT为标准应当是合理的。
（1）首先验证单点模式
我选取任意一个MOD11产品的文件，使用本程序推算它的左上角的中心经纬度，即行列号均为1，如下

得到的经纬度是69.284和29.996，然后使用MRT来验证
MRT在载入文件时会显示一个”Input File info”，如下：

注意这里的条带号是h24v06, 然后看下面的经纬度范围，左上角（Upper-left Corner）是69.28和 30，和本程序得到的结果做一下对比，差了0.004度，这是因为本程序算的是左上角网格的中心点的经纬度，而MRT给出的是左上角顶点的经纬度，考虑到0.01度差不多是1公里，0.004度差不多是半个网格。所以，认为本程序通过了验证（与后面的wgs84转换相比较，可以认为这个结果已经很精确了）。

（2）验证文件模式
这次使用MOD09GA产品，里面有两种维度的数据

使用Arcgis打开500米分辨率的文件（即2400×2400）：

上图展示的是经度文件
这次我们验证右下角，游标点击右下角，得到经纬度图层的栅格值：

然后对照MRT中右下角顶点的值：

同前所述，MRT算的是顶点的坐标，本程序算的是右下角网格的中心点坐标，两相比较，认为是精确的。
综上，本程序通过了验证。

三、 关于sin投影和wgs84的讨论
之前，我一贯的认知是对于经纬度，就直接使用wgs84坐标系就OK了，所以一开始是想简单地做投影转换。而且有一种侥幸的想法是：做wgs84转换后的经纬度应该和“原始”（或者叫正确）的经纬度差不多，but，，，这世上有很多事物是经不起考验的-_-!，以本程序为例：
使用MRT的算法根据sin投影反推得到的经纬度是：

69.284 和 29.995

然后，转换到wgs84的话：

69.183 和 29.962

经度差了0.101度，纬度差了0.033度，这是什么概念呢，对于MOD11产品，相当于经向偏了10几个栅格，纬向偏了4、5个栅格，我认为影响还是比较大的（毕竟做了基准面变换）。

所以，我在本程序中提供了一个可选项，勾选WGS84，则程序会计算wgs84转换后的经纬度。
到底使用哪种方式我觉得是值得斟酌的，应该是跟具体问题相关的，在我的小论文中，经纬度需要有一个可靠的来源，所以我选择使用MRT中的函数直接基于sin投影反推。我可以保证的是，如果不勾选WGS84，本程序是完全准确地按照MRT推算经纬度的算法来计算的
OK，程序就介绍到这里，它的下载地址是：

http://pan.baidu.com/s/1pJkMvgv

最近做了很多工作，这个小工具是其中一个“娱乐项目”，为了完成我的一篇小论文，我还写了一个更加实用的小工具，会在月底论文投出之前整理发布~