# How to figure the MODIS snow cover against different elevation zones?

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};
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' ;

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.

• Can you please add more details? Like how you generated the elevation zones poligons? Whats your source elevation data? You can also put pieces of code and links to the data files you are working on, so we can help better. In a first instance I would think that one the datasets are not geolocated properly, or the resolution of the elevation data is too coarse. Feb 2, 2019 at 15:43
• Dear Camilo Raza, the question is edited and added all asked things. Looking for your valuable suggestions and help Feb 4, 2019 at 3:12
• The code looks incomplete. What is "tp5k6k"? Can you add the figure that would result from: figure; imagesc(data(:,:,1)); hold on; plot(S(1).X,S(1).Y,'r'); Feb 4, 2019 at 3:24
• Updated my question. kindly have a look Feb 4, 2019 at 3:31
• Those are not the files in the script. I mean shapefile.shp, testlon.mat, testlat.mat and testarray.mat Feb 4, 2019 at 5:10

The problem (if there is only one) is in the way you are importing the data into .mat file. If with the files you provided I do

data = load('sc_8days.mat') ;
data = data.sc8days ;
[nx,ny,d] = size(data) ;
imagesc(data(:,:,1));


I get this nonsense:

and d=23, while the MOD10C2 data has only four bands: Eight_Day_CMG_Snow_Cover, Eight_Day_CMG_Cloud_Obscured, Eight_Day_CMG_Clear_Index, and Snow_Spatial_QA (see product specifications).

However, if I load the data and plot directly from the .hdf file you shared, with

data = hdfread('MOD10C2.A2000185.006.2016068190452.hdf','Eight_Day_CMG_Snow_Cover');
imagesc(data);


I get:

Which does makes sense and show the actual data.

Therefore, the way you are importing the data into your .mat file is messing things up.

The error is in this line

B2(i,j,1:24) = reshape(Eight_Day_CMG_Snow_Cover{i,j},[1 3 2]);


Because you are assigning the pixel i,j of each file using data from one file only, and also that should generate and error, I'm not quite sure how this code can run, as the size of Eight_Day_CMG_Snow_Cover is Kx1 and you are attempting to access elements i,j.

Anyway, that way is too cumbersome, avoid reshaping and using cell arrays by importing like this:

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 = zeros(3600,7200,num_files);
Eight_Day_CMG_Clear_Index = zeros(3600,7200,num_files);
Eight_Day_CMG_Cloud_Obscured = zeros(3600,7200,num_files);
Snow_Spatial_QA = zeros(3600,7200,num_files);

for K = 1 : num_files
this_file = filenames{K};