Matlab scripts to generate hill-shaded (3D effect) orienteering maps

These scripts allow you to generate hill-shaded orienteering maps using Matlab. The main script (generateHillShadeMap.m) takes an input image which can be downloaded using Fetch_map python script. It also uses a depth information in xyz ascii format, which can be downloaded from the open data service of National Land Survey of Finland (Elevation model 10 m, height system N2000, XYZ). The second script (loadXYZtoMat.m) is used to join and pre-process all xyz-ascii data files for the main script. The main script requires 'wgs2utm' function by Alexandre Schimel and additional 'loadMAPFile.m' function to parse calibration data from *.map file generated by Fetch_map.

Fig. 1. An example of 5x5km (1:16000) hill-shaded map using different diffuse shading and aspect shading combinations

Usage example:

An example of script usage to make hill-shaded map for a part of Urho Kekkonen National Park using Python and Matlab.

python fetch_map.py --corners -s rka ukk2 7566000 534500 7596000 564500 1:16000

By using the above example, Fetchmap generates 'ukk2.png' and 'ukk2.map' files that are both needed in the following Matlab scripts. The scripts are supposed to be run at the same directory where map files and xyz data files are. The script then produces hill-shaded map ('ukk2-hillshade.png').

loadXYZtoMat('./','./');
generateHillShadeMap('ukk2.png', 'all_xyz.mat', 'ukk2-hillshade.png');

How to set maps to android phone or tablet

First download OruxMaps program to your device. Then load OruxMapsDesktop program to your computer. Set the map file (generated by Fetch_map) as calibration file and just generated hillshaded map as image file. Select a suitable destination directory and create map. After a while, the just created map (a whole folder) can be moved into your mobile devices 'mapfiles' folder which is inside oruxmaps folder. If there is too little space left on your device, you can configure oruxmaps to use some folder at external memory card as map folder. After the map folder has been moved to your phone or tablet, you can use OruxMaps to browse for the just generated map and select it. OruxMaps can be used to navigate on the hillshaded map and record travelled path.

generateHillShadeMap.m

Version 2.0, Thanks to Jakke Kulovesi for implementing aspect shading and more exact hillshading algorithms.

%============================================================================
% Copyright (C) 2014, Heikki Hyyti
%
% Modified (06/2014) by Jakke Kulovesi to implement the method from:
% http://www.reliefshading.com/analytical/shading_methods.html
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the "Software"),
% to deal in the Software without restriction, including without limitation
% the rights to use, copy, modify, merge, publish, distribute, sublicense,
% and/or sell copies of the Software, and to permit persons to whom the
% Software is furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
% FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
% DEALINGS IN THE SOFTWARE.
%============================================================================
%
% generateHillShadeMap(inputMapFilename, inputXYZMatFilename, outputFilename)
%
% Description:
%    Converts map and depthmap and generates hillshaded output map image.
%    Needs wgs2utm by Alexandre Schimel, which can be downloaded from there:
%    http://www.mathworks.com/matlabcentral/fileexchange/14804-wgs2utm--version-2-
%
% Input:
%    inputMapFilename: Filename of input map image (got using fetch_map python script).
%    inputXYZMatFilename: Filename of XYZ.mat file generated using
%       loadXYZtoMat.m function.
%    outputFilename (optional): Filename to save resulting map image.
%    slopeMinMax (optional): To enable aspect shading of deep slopes. The value 
%       is a vector of slope minimum value in radians (default: 0) and Slope 
%       maximum value in radians (default: pi/4). If only aspect shading 
%       (without diffuse shading is needed, use minimum and maximum value of 0.
%
function generateHillShadeMap(inputMapFilename, inputXYZMatFilename, outputFilename, slopeMinMax)
    % compute name for map file from image filename (change ending)
    [pieces, delims] = regexp(inputMapFilename, '\.', 'split', 'match');
    j = [pieces(1:end-1); [delims(1:end-1) {''}]];
    calibrationFilename = [[j{:}] '.map'];
    
    % test input filenames
    if (~exist(inputMapFilename,'file'))
        disp(['Input Map image (' inputMapFilename ') not found!']);
        return;
    elseif (~exist(calibrationFilename,'file'))
        disp(['Input Map calibration file (' calibrationFilename ') not found!']);
        return;
    elseif (~exist(inputXYZMatFilename,'file'))
        disp(['Input depth data file (' inputXYZMatFilename ') not found!']);
        return;
    end
    
    % Load map from image file (got using fetch_map python script)
    map = imread(inputMapFilename);
    N_N = size(map,1); 
    N_E = size(map,2); 
    aspectShading = false;
    
    if (exist('slopeMinMax','var') && length(slopeMinMax) > 1)
        slopeMin = slopeMinMax(1);
        slopeMax = slopeMinMax(2);
        aspectShading = true;
    end
    
    % values from calibration map-file (got using fetch_map python script)
    calibration = loadMAPFile(calibrationFilename);

    Lon = calibration.MMPLL(:,1);
    Lat = calibration.MMPLL(:,2);
    [MMP_E,MMP_N] = wgs2utm(Lat,Lon,calibration.utmzone,'N'); %todo: N / S

    maxE = max(MMP_E);
    minE = min(MMP_E);
    maxN = max(MMP_N);
    minN = min(MMP_N);

    % compute depth map from mat-file generated using loadXYZtoMat.m
    data = load(inputXYZMatFilename);
    dxy = 10; %native resolution of xyz data

    N_x = ceil( (maxE-minE) / dxy );
    N_y = ceil( (maxN-minN) / dxy );

    depthMap = nan(N_y,N_x);

    %filter out z=-999 and x and y indexes which are out of area
    accepted = data.z > -999 & data.x >= minE & data.x <= maxE & data.y >= minN & data.y <= maxN;

    x = data.x(accepted);
    y = (minN+maxN) - data.y(accepted); %reverse y (same as map image)
    z = data.z(accepted);

    x_idx = 1 + floor((x - minE) / dxy);
    y_idx = 1 + floor((y - minN) / dxy);

    for i = 1:size(x,1)
        depthMap(y_idx(i), x_idx(i)) = z(i);
    end

    %remove some variables to reduce memory usage
    clear x y z x_idx y_idx accepted;

    % compute hillshade (exact diffuseShade implemented by Jakke Kulovesi)
    delta_x = [depthMap(:,2:end) - depthMap(:,1:end-1) zeros(size(depthMap,1),1)] ./ dxy;
    delta_y = [depthMap(2:end,:) - depthMap(1:end-1,:); zeros(1, size(depthMap,2))] ./ dxy;
    %hillshade = round(128 + 128*(delta_y+delta_x)); %(dummy approximation using gradients)
    diffuseShade = (1/sqrt(3)).*(delta_x + delta_y + 1)./sqrt(delta_x.^2 + delta_y.^2 + 1);
    diffuseShade = diffuseShade - (1/sqrt(3)); % reduce mean from diffuseShade
    
    if (aspectShading)
        % compute hillshade (aspectShade implemented by Jakke Kulovesi)
        lightingAzimuth = 1/4*pi;	
        steepness = abs(atan(delta_x + delta_y));    
        aspectShade = 0.5 .* (cos(atan2(delta_x,delta_y) - lightingAzimuth) + 1);
        aspectShade = 0.5*(aspectShade - 0.5); %reduce mean from aspectShade and reduce gain

        mask = (steepness-slopeMin) / slopeMax;
        mask(mask < 0) = 0;
        mask(mask > 1) = 1;

        % combine the two different shading methods
        diffuseShade = diffuseShade .* (1-mask) + aspectShade .* mask;
    end
    
    %compute hillshade image from shading(s)
    hillshade = round(128 + 256*diffuseShade); 
     
    nanidx = isnan(hillshade);
    hillshade(nanidx) = 128;
    
    %remove some variables to reduce memory usage
    clear diffuseShade aspectShade steepness mask nanidx;
    
    % resize hillshade to same resolution as map
    hillshade = imresize(uint8(hillshade), [N_N N_E]);

    % integrate hillshade to map image
    gain = 1.2;
    for i = 1:3
        map(:,:,i) = uint8(1.5*((double(map(:,:,i)).*(0.3342+double(hillshade)/383)).^gain)./(255^(gain-1)));
    end

    % save image to file
    if (exist('outputFilename','var') && ~isempty(outputFilename))
        imwrite(map, outputFilename);
    end

    % plot result
    figure(1); clf;
    imshow(map);
    axis equal;
end
    

loadXYZtoMat.m

%============================================================================
% Copyright (C) 2014, Heikki Hyyti
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the "Software"),
% to deal in the Software without restriction, including without limitation
% the rights to use, copy, modify, merge, publish, distribute, sublicense,
% and/or sell copies of the Software, and to permit persons to whom the
% Software is furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
% FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
% DEALINGS IN THE SOFTWARE.
%============================================================================
%
% loadXYZtoMat(sourceDir, outputDir)
%
% Description:
%    Loads all *.xyz files in the source directory and generates one mat
%    file that has the same xyz data.
%
% Input:
%    sourceDir: Directory to search *.xyz files.
%    outputDir: Directory to write 'all_xyz.mat' file.
%
function loadXYZtoMat(sourceDir, outputDir)
    
    if (sourceDir(end) ~= '/')
        sourceDir = [sourceDir '/'];
    end

    if (outputDir(end) ~= '/')
        outputDir = [outputDir '/'];
    end


    inputFiles = dir([sourceDir '*.xyz']);
    
    x = [];
    y = [];
    z = [];
    
    for i = 1:size(inputFiles,1)
        data = load([sourceDir inputFiles(i).name]);
        
        x = [x; data(:,1)];
        y = [y; data(:,2)];
        z = [z; data(:,3)];
    end
    
    save([outputDir 'all_xyz.mat'], 'x', 'y', 'z');
end

loadMAPFile.m

%============================================================================
% Copyright (C) 2014, Heikki Hyyti
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the "Software"),
% to deal in the Software without restriction, including without limitation
% the rights to use, copy, modify, merge, publish, distribute, sublicense,
% and/or sell copies of the Software, and to permit persons to whom the
% Software is furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
% FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
% DEALINGS IN THE SOFTWARE.
%============================================================================
%
% calibration = loadMAPFile(filename)
%
% Description:
%    Loads UTM zone and Map corner coordinates in WGS84 from calibration
%    file (map) generated by fetch_map python script.
%
% Input:
%    filename: Filename of calibration file (got using fetch_map python script).
%
% Output:
%    calibration: Calibration data (currently only utmzone and MMPLL params).
%
function calibration = loadMAPFile(filename)
    fileID = fopen(filename);
    
    calibration = struct('utmzone', [], 'MMPLL', []);
    
    while (true)
        tline = fgetl(fileID);
        if (tline == -1) 
            break;
        end
        strings = textscan(tline, '%s', 'Delimiter', ', ', 'MultipleDelimsAsOne', 1);
        if (strcmp(strings{1}(1),'Projection'))
            calibration(1).utmzone = floor(str2double(strings{1}(4))/6)+31;
        elseif (strcmp(strings{1}(1),'MMPLL'))
            calibration(1).MMPLL = [calibration(1).MMPLL; ...
                str2double(strings{1}(3)) str2double(strings{1}(4))];
        end
    end
    fclose(fileID); 
end