Displaying Single Pixel Spectra from EnMAP Hyperspectral Images with MATLAB

On April 1, 2022, a satellite was launched as part of the German Environmental Mapping and Analysis Program (EnMAP) to collect hyperspectral data, which can be obtained free of charge. Here is a MATLAB script to import and display the data and display the spectra of single pixels selected from an EnMAP image.

The EnMAP satellite acquires 242 spectral bands in the range between 420 nm and 2450 nm with a ground resolution of 30 m x 30 m. The data are free of charge, but you have to apply for permission to include new images. The free EnMAP-Box is available for use with QGIS for data evaluation.

However, the data can also be analyzed with MATLAB, for which there is a Hyperspectral Imaging Library available for the Image Processing Toolbox, released by MathWorks. As an example, we download an image that we took as part of our project in Chew Bahir in the southern Ethiopian Rift and display the spectra of individual pixels selected from the image.

We first define the number of pixels to be selected and the degree of the polynomial used for baseline correction.

ni = 6;
np = 0;

We then create the list of files contained in our directory using

delete('*SW.tiff')
filelist_meta = dir(fullfile('*.XML'));
filename_meta = {filelist_meta.name};
filelist_data = dir(fullfile('*.TIF'));
filename_data = {filelist_data.name};

We then choose the ID of the file to work with using

i = 1;

We then read the wavelengths using

opts = xmlImportOptions(...
   "VariableSelectors",...
   "//wavelengthCenterOfBand");
wavelengths = readtable(filename_meta{i},opts);
wavelengths = table2array(wavelengths);
wavelengths = str2double(wavelengths);

We then read the spectral data

hcube = geohypercube(filename_data{i});
datacube = gather(hcube);

We then extract the geospatial information using

R = hcube.Metadata.RasterReference;
xlimits = R.XWorldLimits;
ylimits = R.YWorldLimits;

We then extract spectral bands 38, 23, 12 using

I1(:,:,1) = datacube(:,:,38);
I1(:,:,2) = datacube(:,:,23);
I1(:,:,3) = datacube(:,:,12);

We then convert the image data to double to perserve the NaNs in the image. Converting the image to uint16 would replace NaNs by zeros.

I2 = double(I1);
I2(I2 == -32768) = NaN;

We then enhance the image by scaling image data to the range [0,1].

I3 = I2/max(I2(:));

Since the image values have a long tail towards the right, a simple way to make the histogram more uniform is to take the square root of the values.

I4 = sqrt(I3);
I4 = abs(I4);

The image is further enhanced using

I4 = imadjust(I4,[0.1 0.6],[]);

We then display the image using

figure('Position',[100 800 800 800])
imshow(I4)
axis on

We then mark and display single pixels using

clear position
for i = 1 : ni
    [position(i,1),position(i,2)] = ginput(1);
    line(position(i,1),position(i,2),...
        'LineStyle','None',...
        'Marker','o',...
        'MarkerFaceColor','b')
     text(position(i,1),position(i,2),num2str(i),...
       'FontSize',18,...
        'Color','r')
end
position = round(position);
position = double(position);

print -dpng -r300 singlepixelimage.png

We then extract the spectra of the individual pixels and perform the baseline correction using simply np degree  polynomial detrending

for i = 1 : ni
clear pixel
pixel = datacube(position(i,1),position(i,2),:);
medspectra(:,i) = double(pixel);
p = polyfit(wavelengths(wavelengths<1300|...
   wavelengths>1450),...
   medspectra(wavelengths<1300|...
   wavelengths>1450,i),np);
   medspectra(:,i) = ...
   medspectra(:,i)-polyval(p,wavelengths);
end

Water bands should be NaNs according to the EnMAP FAQs using

medspectra(wavelengths>1330 & ...
   wavelengths<1450,:) = NaN;
medspectra(wavelengths>1795 & ...
   wavelengths<1970,:) = NaN;

We display the spectra of ni pixels

for i = 1 : ni
legstr(i,:) = ['Pixel ', sprintf('%d',i)];
end

figure('Position',[100 800 1600 600])
axes('Position',[0.1 0.1 0.8 0.8],...
   'XLim',[400 2500],...
   'XGrid','on',...
   'YGrid','on',...
   'XMinorGrid','on',...
   'YMinorGrid','on')
for i = 1 : ni
   line(wavelengths,medspectra(:,i)+100*(i-1),...
      'LineWidth',1)
end
legend(legstr,...
   'Location','southeast',...
   'Box','off')
xlabel('Wavelength (microns)')

print -dpng -r300 singlepixelspectra.png

Download example data and script

Download the EnMAP data used in the example.

References

Guanter, L., Kaufmann, H., Segl, K., Foerster, S., Rogass, C., Chabrillat, S., Kuester, T., Hollstein, A., Rossner, G., Chlebek, C., Straif, C., Fischer, S., Schrader, S., Storch, T., Heiden, U., Mueller, A., Bachmann, M., Mühle, H., Müller, R., Habermeyer, M., Ohndorf, A., Hill, J., Buddenbaum, H., Hostert, P., Linden, S. van der, Leitão, P., Rabe, A., Doerffer, R., Krasemann, H., Xi, H., Mauser, W., Hank, T., Locherer, M., Rast, M., Staenz, K., Sang, B., 2015. The EnMAP Spaceborne Imaging Spectroscopy Mission for Earth Observation. Remote Sensing 7, 8830–8857. https://doi.org/10.3390/rs70708830