MATLAB提取BSQ格式文件经纬度

lq_nuist

哈~大前天为了解决一个月前发的提问帖（至今无人问津。。。）用Matlab提取三维二进制文件所需站点数据http://bbs.06climate.com/forum.p ... 58234&fromuid=72509(出处: 气象家园)
困扰我两个月的问题，从MATLAB小白到自己完整的编一个程序并且调试没有错误！（数据量不大，虽然用EXCEl两个小时也能解决，我还是强迫自己花了整整一天的时间用MALTAB算出来了，嘻）
小白上手很痛苦，基础知识都不懂，还好网上资源广，左拼右凑终成功，怕我过了一阵就忘了，资源再找起来又比较费事，所以把有用的经验资料啥的发在帖子里，以便今后查找用。两个月的烦恼终于撂了，激动的我关了电脑浪了整整两天没有学习，今日上线一看帖子咋地屏蔽了，唉，家园新人啥都不懂。。。怕我过一阵子又忘了当时咋编出来的，趁现在记得，总结一下过程（自己写的）。。。

———————————————————正经说事分割线———————————————————
问题：现有全国13年的GPP.dat和NPP.dat的二进制数据及其头文件，需要提取出指定16个站点及其周围九宫格上共9个点的数据。
          ENVI头文件如下：
          samples = 3080
          lines   = 2186
          bands   = 13
          header offset = 0
          file type = ENVI Standard
          data type = 12
          interleave = bsq
          sensor type = Unknown
          byte order = 0
          map info = {Geographic Lat/Lon, 1.0000, 1.0000, 67.80228509, 54.67476805, 2.2305700000e-002, 1.8023360000e-002, WGS-84, units=Degrees}
          coordinate system string =            {GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]}
wavelength units = Unknown
从头文件中可以获知：该数据的排列方式为BSQ格式（PS：遥感影像的三种存储格式：BSQ是按波段保存，也就是一个波段保存后接着保存第二个波段；BIL是按行保存，就是保存第一个波段的第一行后接着保存第二个波段的第一行，依次类推；BIP是按像元保存，即先保存第一个波段的第一个像元，之后保存第二波段的第一个像元，依次保存。）地理坐标系为WGS-84，投影坐标系为经纬度，一共有13个波段，每个波段是处理后的年GPP或NPP值，13个波段即13年的数值，BSQ格式说明数据是按年排列的，每一个波段的图像上共有2186行，3080列，最左上角点的经纬度为 67.80228509°E，54.67476805°N，与经度平行方向的步长为2.2305700000e-002°，与纬度方向平行的步长为1.8023360000e-002°。
为了更清晰地描述数据结构，示意图如下：

图1

方案一：直接用MATLAB读取二进制文件并将数据提取出来 gpp&npp_extract_point_script.m (1.92 KB, 下载次数: 0, 售价: 2 贡献)

2018-1-26 18:25 上传

点击文件名下载附件
阅读权限: 20 。

提取结果值都在0到255之间，实际值在0~100000之间，不知程序错在哪里，按照书本上的打开二进制文件写的，程序也没有错误，但是结果不对。思考无果，暂时搁置。

方案二：利用ENVI的ROI工具直接提取图像的像素值。为了避免抄袭，此处贴出我参考的资料的原链接http://blog.sina.com.cn/s/blog_764b1e9d0101dqte.html。

根据教程，要提取指定点的像素值，首先需要将点的经纬度整理成以下格式：lat_lon.txt

图2

现已知16个站点的经纬度，还需要计算出九宫格上其他站点的经纬度，虽然只有16个站点，每个站点加减X或Y Pixel Size即可算出经纬度，但在......诶，我突然发现在Excel中计算很快也很方便哦，九个站点的经纬度分别为（lon - Y Pixel Size, lat - X Pixel Size),（lon - Y Pixel Size, lat )（lon - Y Pixel Size, lat + X Pixel Size), （lon , lat - X Pixel Size), （lon , lat), （lon , lat + X Pixel Size), （lon + Y Pixel Size, lat - X Pixel Size), （lon + Y Pixel Size, lat), （lon + Y Pixel Size, lat + X Pixel Size)，X和Y Pixel Size 都是已知的，已知站点经纬度直接Excel拉一下就可以求出九宫格的经纬度啦，所以说我为啥还花了一天的时间用MATLAB，用行列号来间接地求经纬度呢......好吧，只能自我安慰在只知道行列号的情况下，用图1的公式可以得到任意点的经纬度并保存为图2的格式。

利用MATLAB计算任一点的经纬度程序： lat_lon.m (2.11 KB, 下载次数: 0, 售价: 2 贡献)

2018-1-26 18:56 上传

点击文件名下载附件
阅读权限: 20 .

计算好经纬度后，根据http://blog.sina.com.cn/s/blog_764b1e9d0101dqte.html可以导出16*9个站点的GPP和NPP数值。

我们发现导出的数值序号并不是按照MATLAB计算得到的lat_lon.txt中的序号排列的，这是因为ENVI中利用ROI Tool导出的像元信息是按照从上到下，从左到右的“Z”字形顺序，因此用EXCEL打开导出的ASCII文件时，还需要根据经纬度做相应的调整，标注每个站点对应的经纬度及其周边9个点的数值，再利用EXCEl或者SPASS等做数据处理，画图等。

CuiGuoSheng

佩服你的学习精神！赞！