如何算MODIS积雪对不同海拔区域?

Question

我致力于MOD10C2积雪(8天),当我夹MODIS积雪对不同区域的海拔高度对.i.e总趋势的结果。高海拔低显示积雪和低海拔显示更高的积雪。我用matlab inpolygon函数(inpolygon (yq xq,十五,青年志愿))我无法找出犯错误。

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

hdfvars = {‘Eight_Day_CMG_Snow_Cover’,‘Eight_Day_CMG_Clear_Index’,……“Eight_Day_CMG_Cloud_Obscured”、“Snow_Spatial_QA '};测试管理扫描= ' E: \ \ hdf测试”;dinfo = dir (fullfile(管理扫描' * .hdf '));num_files =长度(dinfo);文件名= fullfile(管理扫描,{dinfo.name});Eight_Day_CMG_Snow_Cover =细胞(num_files, 1);Eight_Day_CMG_Clear_Index =细胞(num_files, 1);Eight_Day_CMG_Cloud_Obscured =细胞(num_files, 1);Snow_Spatial_QA =细胞(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 (3600、7200、24);因为我= 1:3600 j = 1:7200 B2 (i, j,一24)=重塑(Eight_Day_CMG_Snow_Cover {i, j}, [1 3 2]);结束结束

然后生成纬度经度和子集感兴趣的区域

朗= 180:0.05:180;lat = 90:0.05:90;subsetqa = Eight_Day_CMG_Snow_Cover (2001:2817、4801:5741:);

然后试图提取高程多边形,我试着(月度)和(8天)部分积雪数据,海拔从SRTM DEM提取多边形。我试图以这种方式使用matlab

shapefile = ' shapefile。轴马力”;S = shaperead (shapefile);N =长度(S);i = 1: N的阴谋(S(我)方式,(我).Y)抓住结束% %朗=负载(“testlon.mat”);朗=经度。testlon;lat =负载(“testlat.mat”);lat = lat。testlat;(X, Y) = meshgrid(经度、纬度);data =负载(“testarray.mat”);data =数据。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

然后我希望以这种方式转化为矩阵输出

测试= cell2mat (cellfun (@transpose,我想,“统一”,0));

注意:这张照片显示出只有一个海拔的多边形区域,我们有六个不同海拔区域。

Answer 1

(如果只有一个)的问题是在你导入数据到.mat文件。如果提供的文件我做

data =负载(“sc_8days.mat”);data =数据。sc8days;[nx、纽约、d] =大小(数据);显示亮度图像(数据(:,:1));

我得到这个废话:

和d = 23, MOD10C2数据只有四个乐队:Eight_Day_CMG_Snow_Cover, Eight_Day_CMG_Cloud_Obscured 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,你正试图访问元素,j。

不管怎么说,这样太麻烦,避免重塑和使用细胞阵列通过导入:

hdfvars = {‘Eight_Day_CMG_Snow_Cover’,‘Eight_Day_CMG_Clear_Index’,……“Eight_Day_CMG_Cloud_Obscured”、“Snow_Spatial_QA '};测试管理扫描= ' E: \ \ hdf测试”;dinfo = dir (fullfile(管理扫描' * .hdf '));num_files =长度(dinfo);文件名= fullfile(管理扫描,{dinfo.name});Eight_Day_CMG_Snow_Cover = 0 (3600、7200、num_files);Eight_Day_CMG_Clear_Index = 0 (3600、7200、num_files);Eight_Day_CMG_Cloud_Obscured = 0 (3600、7200、num_files);Snow_Spatial_QA = 0 (3600、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。