如何计算不同高程区域的MODIS积雪?

Question

我正在研究MOD10C2(8天积雪)，当我针对不同海拔区域剪辑MODIS积雪时，结果与总体趋势相反。海拔越高，积雪越少，海拔越低，积雪越高。我用matlab做了inpolygon函数，(inpolygon(xq,yq,xv,yv))我不知道我哪里出错了。

首先，我以这种方式导入数据

hdfvars = {'Eight_Day_CMG_Snow_Cover'， 'Eight_Day_CMG_Clear_Index'，…“Eight_Day_CMG_Cloud_Obscured”、“Snow_Spatial_QA '};projectdir = 'E:\test\hdf test';Dinfo = dir(fullfile(projectdir， '*.hdf'));Num_files = length(dinfo);文件名= fullfile(projectdir， {dinfo.name});8_day_cmg_snow_cover = cell(num_files, 1);第八条day_cmg_clear_index = cell(num_files, 1);eight_day_cmg_cloud_hazy = cell(num_files, 1);Snow_Spatial_QA = cell(num_files, 1); for K = 1 : num_files this_file = filenames{K}; Eight_Day_CMG_Snow_Cover{K} = hdfread(this_file, hdfvars{1}); Eight_Day_CMG_Clear_Index{K} = hdfread(this_file, hdfvars{2}); Eight_Day_CMG_Cloud_Obscured{K} = hdfread(this_file, hdfvars{3}); Snow_Spatial_QA{K} = hdfread(this_file, hdfvars{4}); end

然后以这种方式重塑数据

B2 = 0 (3,600,7200,24);B2(i,j,1:24) =重塑(Eight_Day_CMG_Snow_Cover{i,j}，[1 3 2]);结束结束

然后生成lat lon并对感兴趣的区域进行子集

Lon = -180:0.05:180;Lat = -90:0.05:90;subsetqa = 8_day_cmg_snow_cover (2001:2817,4801:5741，:);

然后尝试用高程多边形进行提取，我尝试了[月]和[8天]分数积雪数据，高程多边形提取自SRTM DEM。我尝试用matlab这样做

Shapefile = ' Shapefile . '轴马力”;S = shaperead(shapefile);N =长度(S);for i = 1:N plot(S(i).X,S(i).Y) hold on end %% lon = load('testlon.mat');Lon = Lon。testlon;Lat = load(' testat .mat');Lat = Lat。testlat;[X,Y] =网格(lon,lat);Data = load('testarray.mat');数据=数据。testarray; [nx,ny,d] = size(data) ; %%Extract data iwant = cell(d,N) ; for i =1:d A = data(:,:,i) ; for j = 1:N idx = inpolygon(X(:),Y(:),S(i).X,S(i).Y) ; iwant{i,j} = A(idx) ; end

然后用这种方法将输出转换成矩阵

Test = cell2mat(cellfun(@转置，iwant，'uniform'，0));

注意:这张照片只显示了一个高程带的多边形，我们有六个不同的高程带。

Answer 1

问题(如果只有一个)在于您将数据导入到.mat文件的方式。如果有你提供的文件，我知道

Data = load('sc_8days.mat');数据=数据。sc8days;[nx,ny,d] = size(data);显示亮度图像(数据(:,:1));

我知道这是无稽之谈:

， d=23，而MOD10C2数据只有四个波段:Eight_Day_CMG_Snow_Cover、eight_day_cmg_cloud_hazy、Eight_Day_CMG_Clear_Index和Snow_Spatial_QA(参见产品规格）.

但是，如果我直接从您共享的.hdf文件加载数据并绘制

data = hdfread('MOD10C2.A2000185.006.2016068190452.hdf'，'Eight_Day_CMG_Snow_Cover');显示亮度图像(数据);

我得到:

这很有意义，并显示了实际的数据。

因此，您将数据导入.mat文件的方式会把事情搞砸。

错误在这一行

B2 (i, j,一24)=重塑(Eight_Day_CMG_Snow_Cover {i, j}, [1 3 2]);

因为您只使用来自一个文件的数据为每个文件分配像素i,j，并且应该生成和错误，我不太确定这段代码如何运行，因为Eight_Day_CMG_Snow_Cover的大小是Kx1，并且您正在尝试访问元素i,j。

不管怎样，这种方式太麻烦了，通过像这样导入避免重塑和使用单元格数组:

hdfvars = {'Eight_Day_CMG_Snow_Cover'， 'Eight_Day_CMG_Clear_Index'，…“Eight_Day_CMG_Cloud_Obscured”、“Snow_Spatial_QA '};projectdir = 'E:\test\hdf test';Dinfo = dir(fullfile(projectdir， '*.hdf'));Num_files = length(dinfo);文件名= fullfile(projectdir，{dinfo.name});Eight_Day_CMG_Snow_Cover = 0 (3,600,7200,num_files);Eight_Day_CMG_Clear_Index = 0 (3,600,7200,num_files);eight_day_cmg_cloud_hazy = 0 (3,600,7200,num_files);Snow_Spatial_QA = 0 (3,600,7200,num_files); for K = 1 : num_files this_file = filenames{K}; Eight_Day_CMG_Snow_Cover(:,:,K) = hdfread(this_file, hdfvars{1}); Eight_Day_CMG_Clear_Index(:,:,K) = hdfread(this_file, hdfvars{2}); Eight_Day_CMG_Cloud_Obscured(:,:,K) = hdfread(this_file, hdfvars{3}); Snow_Spatial_QA(:,:,K) = hdfread(this_file, hdfvars{4}); end

因此在本例中，变量Eight_Day_CMG_Snow_Cover包含您想要放入B2的内容。