用Matlab提取三维二进制文件所需站点数据

lq_nuist

已知投影后的gpp和npp头文件信息如下： gpp_geog.hdr (582 Bytes, 下载次数: 0)

2018-1-5 12:31 上传

点击文件名下载附件
阅读权限: 10 ； npp_geog.hdr (378 Bytes, 下载次数: 0)

2018-1-5 12:31 上传

点击文件名下载附件
阅读权限: 10 ENVI
description = {
  Create New File Result [Thu Dec 07 22:22:43 2017]}
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

需要提取的16个站点经纬度信息： site.xls (28 KB, 下载次数: 0)

2018-1-5 12:36 上传

点击文件名下载附件
阅读权限: 10
其中（i,j）根据头文件最左上角坐标（54.67476805，67.80228509）和步长（1.8023360000e-002，2.2305700000e-002）计算得到。
按我的理解是，gpp数据一共有13个波段，分别是13年的数据，每一层有2186*3080的数据，现在需要提取每一层固定（i,j）的数值。

我自己编写的matlab程序如下： gpp&npp_extract_point_script.m (1.92 KB, 下载次数: 0)

2018-1-5 12:34 上传

点击文件名下载附件
阅读权限: 10
%Extract gpp&npp datas at the fixed location(9) based on latitude and longitude.
%This program is written by Liu Qian,2017-12-22.
%Email:liuqiannuist@163.com.

clear;
clc;
row=1;  % extract from row 1
fid=fopen('gpp_geog.dat','r'); % open binary data and only read. If extract npp data, read npp_geog.dat.
conclusion=cell(16*13,14); % write conclusion into (16 citys * 13 years, 14 columns)
[num,city]=xlsread('site.xls','Sheet1','A:A'); % city{2,1} is the first city.
[city_number,txt]=xlsread('site.xls','Sheet1','H:H'); % city_number from 1 to 16
[lat,txt]=xlsread('site.xls','sheet1','J:J'); % city latitude
[lon,txt]=xlsread('site.xls','sheet1','I:I'); % city longitude
[city_i,int]=xlsread('site.xls','Sheet1','L:L'); % (i,j) is where the city
[city_j,int]=xlsread('site.xls','Sheet1','K:K');

for year=2000:2013
    data=fread(fid,[3080,2186]); % read binary data to [m*n], rank as columns, so (j,i) is where the city
    for number=1:16
         if number==city_number(number,1);
            i=round(city_i(city_number(number,1),1)); % i is int
            j=round(city_j(city_number(number,1),1)); % j is int

            conclusion{row,1}=city{city_number(number,1)+1,1};
            conclusion{row,2}=city_number(number,1);
            conclusion{row,3}=year;
            conclusion{row,4}=lat(number,1);
            conclusion{row,5}=lon(number,1);
            conclusion{row,6}=data(j-1,i-1); %1
            conclusion{row,7}=data(j-1,i); %2
            conclusion{row,8}=data(j-1,i+1); %3
            conclusion{row,9}=data(j,i-1); %4
            conclusion{row,10}=data(j,i); %5
            conclusion{row,11}=data(j,i+1); %6
            conclusion{row,12}=data(j+1,i-1); %7
            conclusion{row,13}=data(j+1,i); %8
            conclusion{row,14}=data(j+1,i+1); %9
            row=row+1;
         end
    end
end
fclose(fid);
xlswrite('conclusion.xlsx',conclusion,'Sheet1','A2');


我得到的结果是： conclusion_gpp.xlsx (25.38 KB, 下载次数: 0)

2018-1-5 12:41 上传

点击文件名下载附件
阅读权限: 10 conclusion_npp.xlsx (25.43 KB, 下载次数: 0)

2018-1-5 12:41 上传

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

由于gpp.dat数据太大，此处无法上传。

发现结果数据不对，可能是我三维数据没有读对，也可能是指针没走对，按照我之前的思想，2186*3080*13不等于175054880个字节，检查了很久，不知道怎么改，希望懂得人能给予解答，多谢。

幽居寒舍

我就是来整点金币的

lq_nuist

自我回复
编程解决不了的问题，那就......不用编程了
经过两个月的随性思考，选择曲线救国的方法，直接用ENVI提取所需站点信息。
————————————正经说事专用分割线————————————————————
1、首先要提取数据，要认识数据结构。
    全国GPP和NPP的ENVI数据。数据类型：二进制后缀为.dat格式文件。
    数据存储格式：BSQ格式，从左到右，从上到下存储，共13个波段，每个波段代表处理过的一年的数据，即一共13年的GPP和NPP数据。地理坐标为经纬度。
2、我要做什么：提取已知经纬度的16个站点及其九宫格上共9个点的GPP和NPP数值。
3、怎么做：利用ENVI中ROI工具获取指定点的像素值。
step1：整理要提取站点的数据，以“序号 纬度 经度” 的格式保存为.txt文件。
step2：在ENVI经典界面中打开GPP.dat图像，选择overlay->Region Of Interest 打开ROI Tool。
           在ROI Tool中，选择ROI_Type->Input Points from ASCII，选择文本格式的实地调查数据反演点.txt。
           注意参数选择：
          x：选择经度；y：选择纬度；
          These point comprise：Individual Points；
          投影坐标（Select Map Based Projection）：Geographic Lat/Lon；
          注：投影坐标与实测数据中坐标值的投影参数保持一致。
          设置好投影信息后，单击OK这样将实地调查的点位置信息加载到图像中。
step3： 在ROI TOOL中，选择File->Output ROIs to ASCII。选择GPP.dat图像。
            在Output ROIs to ASCII Parameters面板中，选择ROI点，单击Edit Output ASCII Form，在输出内容设置面板中，选择ID、经纬度（Geo Location）、和波段像元值（Band Values）。
    这样，就将16个站点的GPP数值导出来了。同样可以计算好九宫格上其他站点的经纬度，重复上面步骤，导出数值。
    我们发现，导出来的经纬度与输入实测点的经纬度不完全一致，这是因为一般影像中像元的坐标是取中心点的经纬度，而我们实测的点位不一定恰好对应着影像像元的中心点，所以当输入的经纬度与影像上单个像元的经纬度不一致时，就会采用就近原则，与最临近的像素点匹配，输出该点的经纬度坐标。如环境小卫星的空间分辨率是30米，换算成经纬度在几秒的范围。
    将相应的像素值导入Excel表中，与实测值一一对应。
注：ENVI中利用ROI Tool导出的像元信息是按照从上到下，从左到右的“Z”字形顺序，因此我们导入的实测点最好事先进行排序，这样导入Excel表中就方便多了，不用一个个根据经纬度手动调整。

lq_nuist

ENVI中从图像中获取指定点的像素值：
http://blog.sina.com.cn/s/blog_764b1e9d0101dqte.html