1
\ begingroup美元

我正在研究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));

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

在这里输入图像描述

\ endgroup美元
11
  • 1
    \ begingroup美元 请介绍具体情况。比如你是如何生成海拔区域多边形的?你的源海拔数据是什么?您还可以将代码片段和链接放到您正在处理的数据文件中,这样我们可以更好地提供帮助。在第一个例子中,我会认为一个数据集没有正确地定位,或者高程数据的分辨率太粗糙。 \ endgroup美元
    - - - - - -卡米洛·Rada
    2019年2月2日15:43
  • \ begingroup美元 亲爱的卡米洛·拉扎,这个问题已经编辑并添加了所有问过的问题。期待您的宝贵建议和帮助 \ endgroup美元
    - - - - - -伊尔凡
    2019年2月4日3:12
  • \ begingroup美元 代码看起来不完整。什么是“tp5k6k”?你能不能加上由:figure;显示亮度图像(数据(:,:1));抓住;情节((1)方式,年代(1).Y, ' r '); \ endgroup美元
    - - - - - -卡米洛·Rada
    2019年2月4日3:24
  • \ begingroup美元 更新了我的问题。请您看一下 \ endgroup美元
    - - - - - -伊尔凡
    2019年2月4日3:31
  • 1
    \ begingroup美元 这些不是脚本中的文件。我指的是shapefile。轴马力,testlon。垫,testlat。Mat和testarray。Mat \ endgroup美元
    - - - - - -卡米洛·Rada
    2019年2月4日5:10

1回答1

2
\ begingroup美元

问题(如果只有一个)在于您将数据导入到.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的内容。

\ endgroup美元
4
  • \ begingroup美元 在导入MOD10C2的所有四个波段后,为了加快工作速度,我在三维上拼接了仅“Eight_Day_CMG_Snow_Cover”的23个文件(23个8days snow cover)。 \ endgroup美元
    - - - - - -伊尔凡
    2019年2月5日2:07
  • \ begingroup美元 有道理,但是你在连接上搞砸了一些东西,上面的第一个图是证明。将导入和连接的代码添加到问题中。 \ endgroup美元
    - - - - - -卡米洛·Rada
    2019年2月5日2:10
  • \ begingroup美元 @irfan看到答案的变化。 \ endgroup美元
    - - - - - -卡米洛·Rada
    2019年2月5日15:33
  • 1
    \ begingroup美元 @亲爱的Camilo Rada,这对我来说很有效。谢谢你的宝贵时间。最后一件事,如何从数据中提取经纬度和经度。比方说,我想要一个MOD10C2的子集,从10到30 N和40到60 E,所以我应该有三个.mat文件。E(1)积雪,(2)纬度,(3)lon作进一步分析。 \ endgroup美元
    - - - - - -伊尔凡
    2019年2月6日8:29

你的答案

点击“张贴您的答案”,即表示您同意我们的服务条款隐私政策而且饼干的政策

这不是你想要的答案?浏览带标签的其他问题问自己的问题