I am working on MOD10C2 (8days snow cover), when I clip the MODIS snow cover against different zones of elevation the results was against general trend .i.e. the higher elevation was showing lower snow cover and lower elevation was showing higher snow cover. I did the with matlab inpolygon function, (inpolygon(xq,yq,xv,yv)) I am could not figure out where I am making mistakes.
First of all I imported the data in this way
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);
filenames = fullfile( projectdir, {dinfo.name} );
Eight_Day_CMG_Snow_Cover = cell(num_files, 1);
Eight_Day_CMG_Clear_Index = cell(num_files, 1);
Eight_Day_CMG_Cloud_Obscured = 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
Then reshaped data in this way
B2 = zeros(3600,7200,24);
for i = 1:3600
for j = 1:7200
B2(i,j,1:24) = reshape(Eight_Day_CMG_Snow_Cover{i,j},[1 3 2]);
end
end
Then generate lat lon and subset the area of interest
lon = -180:0.05:180;
lat = -90:0.05:90;
subsetqa = Eight_Day_CMG_Snow_Cover(2001:2817,4801:5741,:);
then attempted to extract with elevation polygon, I tried both [monthly] and [8days] fractional snow cover data, elevation polygon are extracted from SRTM DEM. I attempted to use matlab in this way
shapefile = 'shapefile.shp' ;
S = shaperead(shapefile) ;
N = length(S) ;
for i = 1:N
plot(S(i).X,S(i).Y)
hold on
end
%%
lon = load('testlon.mat') ; lon = lon.testlon ;
lat = load('testlat.mat') ; lat = lat.testlat ;
[X,Y] = meshgrid(lon,lat) ;
data = load('testarray.mat') ; data = 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
Then output iwant converted into matrix in this way
test = cell2mat(cellfun(@transpose,iwant,'uniform',0));
NOTE: this picture is showing polygon of only one elevation zone and we have six different elevation zones.