The Use of SRTM Data to Visualize my Holiday Resort in the Rhine Rift

During our early summer vacation my family and I will spend a few days near Landau in der Pfalz, where I grew up many years ago. The place is a beautiful town located in the Rhine Rift, east of the Pfälzerwald (Palatinate Forest). Here is a MATLAB script to visualize part of the Rhine Rift, similar to the one used to display the digital elevation model of the Kenya Rift, whose exploration I have devoted most of my career to.

The free 1-arc minute Shuttle Radar Topography Mission (SRTM) data are very popular digital terrain models (DEM). Chapter 7.5 of the MRES book demonstrates how to import, process and display a single SRTM tile. We first load the data set N49E008.hgt as described in the MRES book:

clear, close all, clc
fid = fopen('N49E008.hgt','r');
SRTM = fread(fid,[1201,inf],'int16','b');
fclose(fid);

Next we transpose and flip the data sets up and down, before we replace the 16-bit signed integer no-data identifier -32768 by the MATLAB no-data identifier NaN (not-a-number).

SRTM = SRTM'; SRTM = flipud(SRTM);
SRTM(SRTM == -32768) = NaN;

A popular way to eliminate gaps in digital elevation models is by filling gaps with the arithmetic means of adjacent elements:

for i = 2 : 1200
    for j = 2 : 2400
        if isnan(SRTM(i,j)) == 1
            SRTM(i,j) = nanmean(nanmean(...
            SRTM(i-1:i+1,j-1:j+1)));
        end
    end
end
clear i j

Next we define a coordinate grid.

[LON,LAT] = meshgrid(8:1/1200:9,49:1/1200:50);

Here is the script of a beautiful display of the data. Please note the use of the function daspect to reduce the length of the z-axis.

figure1 = figure('Color',[1 1 1],...
    'Position',[50 50 1200 600]);
axes1 = axes('Units','centimeters',...
    'FontSize',8,...
    'YLim',[49 49.5]);
hold(axes1,'all');
surf1 = surf(LON,LAT,SRTM,...
    'SpecularExponent',20,...
    'FaceLighting','phong',...
    'FaceColor','interp',...
    'EdgeColor','none');
light1 = light('Parent',axes1,...
    'Style','local',...
    'Position',[145 70 900000]);
set(gca,'View',[25 20])
daspect([1 1 10000])
demcmap(SRTM)
plot3(8.118,49.2012,200,...
    'Marker','o',...
    'MarkerSize',6,...
    'MarkerEdgeColor','w',...
    'MarkerFaceColor','w')
text(8.118+0.01,49.2012,200,...
    'Landau in der Pfalz',...
    'Color','w',...
    'FontSize',12)
xlabel('Longitude')
ylabel('Latitude')
zlabel('Altitude (m asl)')