<p>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</p>

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

<p>I get this nonsense:</p>

<p><a href="https://i.sstatic.net/8M5kw.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/8M5kw.png" alt="enter image description here"></a></p>

<p>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 <a href="https://nsidc.org/data/mod10c2" rel="nofollow noreferrer">product specifications</a>).</p>

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

<pre><code>data = hdfread('MOD10C2.A2000185.006.2016068190452.hdf','Eight_Day_CMG_Snow_Cover');
imagesc(data);
</code></pre>

<p>I get:</p>

<p><a href="https://i.sstatic.net/ck4xl.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/ck4xl.png" alt="enter image description here"></a></p>

<p>Which does makes sense and show the actual data.</p>

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

<p>The error is in this line</p>

<pre><code>B2(i,j,1:24) = reshape(Eight_Day_CMG_Snow_Cover{i,j},[1 3 2]);
</code></pre>

<p>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.</p>

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

<pre><code>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};
    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
</code></pre>

<p>So in this case, the variable Eight_Day_CMG_Snow_Cover contains what you wanted to put in B2.</p>

CC BY-SA 4.0
2019-02-05T02:07:10.640Z
After importing all four bands of MOD10C2, I Concatenate 23 files (23 files of 8days snow cover) of only 'Eight_Day_CMG_Snow_Cover' on  3rd dimension in order to speed up the work.
irfan
CC BY-SA 4.0
2019-02-05T02:10:27.953Z
Makes sense, but you messed up something on the concatenation, the first figure above is the proof. Add your code for the import and concatenation to the question.
Camilo Rada
CC BY-SA 4.0
2019-02-05T15:33:47.387Z
@irfan See the change in the answer.
Camilo Rada
CC BY-SA 4.0
2019-02-06T08:29:18.680Z
@Dear Camilo Rada, it worked well for me. Thank you for your time. One last thing, how I can extract the lat and lon from this data. let say, I want a subset of MOD10C2 from 10 to 30 N and 40 to 60 E so I should have three .mat file .i.e (1) snow cover, (2) lat and (3) lon for further analysis.
irfan
CC BY-SA 4.0
2019-02-04T17:46:51.500Z
2019-02-05T19:17:21.467Z
2019-02-05T19:17:21.467Z
Camilo Rada
<p>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.</p>

<p>First of all I imported the data in this way</p>

<pre><code>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
</code></pre>

<p>Then reshaped data in this way</p>

<pre><code>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
</code></pre>

<p>Then generate lat lon and subset the area of interest</p>

<pre><code>lon = -180:0.05:180;
lat = -90:0.05:90;
subsetqa = Eight_Day_CMG_Snow_Cover(2001:2817,4801:5741,:);
</code></pre>

<p>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</p>

<pre><code>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
</code></pre>

<p>Then output iwant converted into matrix in this way</p>

<pre><code>test = cell2mat(cellfun(@transpose,iwant,'uniform',0));
</code></pre>

<p>NOTE: this picture is showing polygon of only one elevation zone and we have six different elevation zones.</p>

<p><a href="https://i.sstatic.net/UOCLB.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/UOCLB.jpg" alt="enter image description here"></a></p>

CC BY-SA 4.0
2019-02-02T15:43:39.673Z
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.
Camilo Rada
CC BY-SA 4.0
2019-02-04T03:12:02.803Z
Dear Camilo Raza, the question is edited and added all asked things. Looking for your valuable suggestions and help
irfan
CC BY-SA 4.0
2019-02-04T03:24:48.493Z
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');
Camilo Rada
CC BY-SA 4.0
2019-02-04T03:31:04.720Z
Updated my question. kindly have a look
irfan
CC BY-SA 4.0
2019-02-04T03:37:28.780Z
I don't see the base image that should be in the picture due to the line imagesc(data(:,:,1)); that should plot the Eight_Day_CMG_Snow_Cover on the image. Maybe try imagesc(lon, lat, data(:,:,1));
Camilo Rada
CC BY-SA 4.0
2019-02-04T03:44:53.563Z
First I loaded the snow cover data (data = testarray) then tried to extract the data of specific elevation zone one by one rather plotting the image of now cover. The snow cover data is in HDF and NC format.
irfan
CC BY-SA 4.0
2019-02-04T03:54:21.180Z
What I want you to do is to plot the polygons on top of the snow cover image, so we can see if the geolocation of both datasets match. Or share all the files used in the script trough dropbox, google drive or something similar.
Camilo Rada
CC BY-SA 4.0
2019-02-04T04:30:25.683Z
http://www.mediafire.com/folder/qo8b5eqd0d3nu/Snowcover    here I have three files 8days snow cover(HDF), monthly snow cover (Netcdf) and SRTM DEM. Kindly find it
irfan
CC BY-SA 4.0
2019-02-04T05:10:29.837Z
Those are not the files in the script. I mean shapefile.shp, testlon.mat, testlat.mat and testarray.mat
Camilo Rada
CC BY-SA 4.0
2019-02-04T09:47:20.367Z
updated, kindly have a look.   http://www.mediafire.com/folder/688zf2kkmu5mu/snowcover
irfan
CC BY-SA 4.0
2019-02-04T09:48:13.803Z
Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/89227/discussion-between-irfan-and-camilo-rada).
irfan
<mountains><remote-sensing><snow><gis>
How to figure the MODIS snow cover against different elevation zones?
irfan
# How to figure the MODIS snow cover against different elevation zones?

(Asked by: username_0 on 2019-02-02)

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.

[![enter image description here](https://i.sstatic.net/UOCLB.jpg)](https://i.sstatic.net/UOCLB.jpg)

**Tags:** mountains, remote-sensing, snow, gis

> **Comments:**
> **username_1**: 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.  
> **username_0**: Dear username_1 Raza, the question is edited and added all asked things. Looking for your valuable suggestions and help  
> **username_1**: 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');  
> **username_0**: Updated my question. kindly have a look  
> **username_1**: I don't see the base image that should be in the picture due to the line imagesc(data(:,:,1)); that should plot the Eight_Day_CMG_Snow_Cover on the image. Maybe try imagesc(lon, lat, data(:,:,1));  
> **username_0**: First I loaded the snow cover data (data = testarray) then tried to extract the data of specific elevation zone one by one rather plotting the image of now cover. The snow cover data is in HDF and NC format.  
> **username_1**: What I want you to do is to plot the polygons on top of the snow cover image, so we can see if the geolocation of both datasets match. Or share all the files used in the script trough dropbox, google drive or something similar.  
> **username_0**: http://www.mediafire.com/folder/qo8b5eqd0d3nu/Snowcover    here I have three files 8days snow cover(HDF), monthly snow cover (Netcdf) and SRTM DEM. Kindly find it  
> **username_1**: Those are not the files in the script. I mean shapefile.shp, testlon.mat, testlat.mat and testarray.mat  
> **username_0**: updated, kindly have a look.   http://www.mediafire.com/folder/688zf2kkmu5mu/snowcover  
> **username_0**: Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/89227/discussion-between-username_0-and-username_1-username_1).  

---

## 1 answers

### Answer by username_1 on 2019-02-04. Score: 2 | Accepted answer

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:

[![enter image description here](https://i.sstatic.net/8M5kw.png)](https://i.sstatic.net/8M5kw.png)

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](https://nsidc.org/data/mod10c2)).

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:

[![enter image description here](https://i.sstatic.net/ck4xl.png)](https://i.sstatic.net/ck4xl.png)

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

```

So in this case, the variable Eight\_Day\_CMG\_Snow\_Cover contains what you wanted to put in B2.

> **Comments:**
> **username_0**: After importing all four bands of MOD10C2, I Concatenate 23 files (23 files of 8days snow cover) of only 'Eight_Day_CMG_Snow_Cover' on  3rd dimension in order to speed up the work.  
> **username_1**: Makes sense, but you messed up something on the concatenation, the first figure above is the proof. Add your code for the import and concatenation to the question.  
> **username_1**: @username_0 See the change in the answer.  
> **username_0**: @Dear username_1, it worked well for me. Thank you for your time. One last thing, how I can extract the lat and lon from this data. let say, I want a subset of MOD10C2 from 10 to 30 N and 40 to 60 E so I should have three .mat file .i.e (1) snow cover, (2) lat and (3) lon for further analysis.  

---

