Importing, Processing and Displaying EnMAP Hyperspectral Images with MATLAB - MATLAB and Python Recipes for Earth Sciences

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.

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 (369.9 MB) that we took as part of our project in Chew Bahir in the southern Ethiopian Rift. We first define the names of the files with the metadata and the spectral data using

clear, clc, close all

filename_metadata = ...
   ['ENMAP01-____L2A-DT0000038475_',...
   '20230831T083041Z_',...
   '003_V010400_202311',...
   '07T163005Z-METADATA.XML'];

filename_spectrum = ...
   ['ENMAP01-____L2A-DT0000038475_',...
   '20230831T083041Z_',...
   '003_V010400_202311',...
   '07T163005Z-SPECTRAL_IMAGE.TIF'];

We then read and store the wavelengths of the spectral bands in wv using

opts = xmlImportOptions(...
   "VariableSelectors",...
   "//wavelengthCenterOfBand");

wavelengths = readtable(filename_metadata,opts);

wv = table2array(wavelengths);
wv = str2double(wv);

clear opts wavelengths

We then read the spectral data using

hcube = geohypercube(filename_spectrum);
datacube = gather(hcube);

We then extract the geospatial information using

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

We then extract bands 38, 23, 12 to be concatenated to a RGB composite I1 using

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

We then convert the 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 by scaling image data to the range [0,1] using

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

We then further enhance the image for display by using an adaptive histogram equalization. Since the method replaces NaNs by zeros, it causes a long list of warnings. We therefore disable this specific type of warning before using adapthisteq.

warning('off','images:imhistc:inputHasNaNs')

I4(:,:,1) = adapthisteq(I4(:,:,1));
I4(:,:,2) = adapthisteq(I4(:,:,2));
I4(:,:,3) = adapthisteq(I4(:,:,3));

We are then ready to display the image using the geospatial information XData and YData:

imshow(I4,...
   'XData',xlimits,...
   'YData',ylimits)
axis on

Download example data and script

Download data and script (369.9 MB)

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