How to find Summary Statistics of Image in Matlab

I've MODIS satellite raster (tiff image) and this raster is of evapotrainspiration (ET) dataset collected by this satellite. I want to find mean ET over my study area, however, this mean value of the raster is different from the mean value of the ET calculated in ESRI software e.g., ArcGIS. I am using this simple code

im=imread('MOD16A2.A2002001.h24v05.006.2017075173934_ET.tif');
imm=im*0.1 %multipling scale factor according to MODIS documentation
imm(imm == 65535) = nan; %datavalues outside the study area as nan
nemean=mean2(imm)

The mean value comes:

nemean =

4.9621

But the actual mean value of ET is

337.422222

i'm sharing ET image with this question. Why mean value of raster comes different in matlab? I think, I'm doing some mistake please correct?

Answer 1

[Image Analyst]

Perhaps you have the wrong image. The images has values only up to a gray level of 71, if you don't count all those values that are saturated at 65535. And the mean is 20.621096. How are you defining your "Study area"? There are no values around 337 or 10 times that in the image.

originalImage = imread('MOD16A2.A2002001.h24v05.006.2017075173934_ET.tif');
% imm=originalImage*0.1 %multipling scale factor according to MODIS documentation
subplot(2, 2, 1);
imshow(originalImage, [])
title('Original 16 bit image');
subplot(2, 2, 2);
[counts, edges] = histcounts(originalImage, 65536);
% Suppress huge spikes at 0 and 65535
counts(1) = 0;
counts(end) = 0;
% Find last non-zero bin
lastBin = find(counts > 0, 1, 'last')
bar(edges(1:lastBin), counts(1:lastBin))
grid on;
title('Histogram')
% Change 65535 to 0.
mask = originalImage < 65535;
maskedImage = originalImage .* uint16(mask);
subplot(2, 2, 3);
imshow(maskedImage, []);
meanGrayLevel = mean(originalImage(mask), 'omitnan')
caption = sprintf('Mean = %f', meanGrayLevel)
title(caption);

[Muhammad Usman Saleem]

Many thanks for your kind reply sir @Image Analyst

I'm using same image in both softwares e.g., Matlab and Arcmap. The image is cliped according to my shaefile (attached) and its pixel values contains ET multiple by 0.1 scale factor. I want to calculate mean ET within my polygon shapefile.

When I calculate summary statistics of ET under field Id 'Pixel_ET' in arcmap its giving me mean value of ET as 337.422222 and when I read the same image in Matlab software its means value of ET is different!.

Arcmap generate some metedata files with tiff image, you can check summary statistics of 'Pixel_ET' field of raster (Esri raster of same image attached). Field 'Pixel_ET' has been calculated by multiplying value*count*0.1

[Image Analyst]

Which of the files in the zip file is the "shaefile"? Is that a binary image or some ROI? What values do you want to compute the mean over? Probably everything that is not 65,535 right? I'll see if tomorrow I can load it into Photoshop and sees if it agrees with MATLAB. I don't know what Arcmap does. Maybe it multiplies the gray level by 10 to convert from gray levels back to ET value. Or maybe it does something different from that. I don't know.

[Muhammad Usman Saleem]

I've shared two zip files ' MOD16A2.A2002001.h24v05.006.2017075173934_ET.zip' file contain Esri tiff file along with its metadata dear and 'Shapefile.zip' which contains ROI (polygon of my area of interest).

Yes, you are right sir! I want to find means of all pixels except 65,535. I want to calculate mean value of pixels fall within my shapefile. I perafer Matlab as it is user friendly. The ultimate purpose of this question, is calculate accurate mean value of pixel values fall within my polygon! Thank you @Image Analyst

[Image Analyst]

My code does that. The mean gray level is 20.62. Since the gray level is 0.1 times the ET value, the mean ET value is 206.2.

[Image Analyst]

@Muhammad Usman Saleem So what else do we need to do here such that you feel like your problem is solved and you can "Accept this answer"?

[Muhammad Usman Saleem]

@Image Analyst Sir, Many thanks for your kind reply. I definately accept your answer as it satisfy me.

One more task I needed, I want to skip pixel values ranging 32761 to 32767 and 65536 while calculating mean of pixel values falling within my shapefile domain. Will you please edit your kind answer so that I skip pixel values ranging from 32761 to 32767 ( by MODIS documentation this range contains buildup landcover) and 65536 (Background), I not need these pixel values in my Mean value. Please modify your kind code.

[Image Analyst]

@Muhammad Usman Saleem, Thanks for your (future) acceptance. Try this:

clc;    % Clear the command window.
fprintf('Beginning to run %s.m ...\n', mfilename);
close all;  % Close all figures (except those of imtool.)
clearvars;
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 18;
originalImage = imread('MOD16A2.A2002001.h24v05.006.2017075173934_ET.tif');
% imm=originalImage*0.1 %multipling scale factor according to MODIS documentation
subplot(2, 3, 1);
imshow(originalImage, [])
impixelinfo;
title('Original 16 bit image');
subplot(2, 3, 2:3);
[counts, edges] = histcounts(originalImage, 65536);
% Suppress huge spikes at 0 and 65535
counts(1) = 0;
counts(end) = 0;
% Find last non-zero bin
lastBin = find(counts > 0, 1, 'last')
bar(edges(1:lastBin), counts(1:lastBin))
grid on;
title('Histogram')
% "I want to skip pixel values ranging 32761 to 32767 and 65536"
% There is no 65536 so I'll assume you mean 65535.
pixelsToExclude = (originalImage >= 32761 & originalImage <= 32767) | (originalImage >= 65535);
subplot(2, 3, 4);
imshow(pixelsToExclude, [])
title('Pixels to Exclude');
% The pixels to consider is just the inverse of the pixels we want to exclude.
mask = ~pixelsToExclude;
subplot(2, 3, 5);
imshow(mask, [])
title('Pixels to Include');
maskedImage = originalImage .* uint16(mask);
subplot(2, 3, 6);
imshow(maskedImage, []);
meanGrayLevel = mean(originalImage(mask), 'omitnan')
caption = sprintf('Mean = %f', meanGrayLevel)
title(caption);