Utilisateur:Ikonact/Cartographie/wmt addtopo

function mapdata = wmt_addtopo(x, y, data, cmap, levels, varargin)
%#WMT_ADDTOPO plot a relief layer on a map.
%#
%# WMT_ADDTOPO(x, y, z, clr, [lvl1 lvl2 lvl3 ..  lvln], 'name', layername);
%# x and y - vectors defining the lat and lon grids
%# z - topographic data matrix
%# clr - color map to be used
%# [lvl1 lvl2 lvl3 ..  lvln] - topographic levels to plot
%# layername - name of the layer (usually 'topo').

%# The WMT tool is a collection of Octave/Matlab scripts able to generate geographic maps images.
%#
%# Copyright 2012 ikonact
%# http://commons.wikimedia.org/wiki/User:Ikonact
%#
%# This file is part of the WMT tool.
%#
%# WMT is free software: you can redistribute it and/or modify
%# it under the terms of the GNU General Public License as published by
%# the Free Software Foundation, either version 3 of the License, or
%# (at your option) any later version.
%#
%# WMT is distributed in the hope that it will be useful,
%# but WITHOUT ANY WARRANTY; without even the implied warranty of
%# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%# GNU General Public License for more details.
%#
%# You should have received a copy of the GNU General Public License
%# along with WMT.  If not, see <http://www.gnu.org/licenses/>.

global wmt_param;
name = 'topo';
shadowlims = [];

for i = 2:2:(nargin-5)
    switch(varargin{i-1})
        case 'name'
            name = varargin{i};
        case 'shadow'
            shadowlims = varargin{i};
        otherwise
            disp(varargin{i});
    end
end

if strcmp(wmt_param.type, 'svg')
    %#%# Topo data in svg format
    [CC hC] = wmtp_drawcontourf(x, y, data, levels);
    colormap(cmap);
    Cld = get(hC, 'Children');

    %#print svg file
    fprintf(wmt_param.fid, '<g id="%s" inkscape:groupmode="layer" inkscape:label="%s" clip-path="url(#IDmain)">\n', name, name);
    fprintf(wmt_param.fid, '<!-- %d %d -->\n', wmt_param.pictsize(1), wmt_param.pictsize(2));
    for j=length(Cld):-1:1
        if strcmp(get(Cld(j), 'Type'), 'patch')
            Iso = get(Cld(j), 'CData');
            if ~isnan(Iso)
                wmtp_plotpath(wmt_param.fid, Cld(j), wmt_param.pictsize, ['class="' name '' int2str(find(levels==Iso)-1) '"']);
            end
        end
    end
    fprintf(wmt_param.fid, '</g>\n');

    %#add shadow
    if ~isempty(shadowlims)
        dd = data(2:end, 1:end-1) - data(1:end-1, 2:end);
        [CC hC] = wmtp_drawcontourf(x(2:end), y(2:end), dd, shadowlims);
        Cld = get(hC, 'Children');
        fprintf(wmt_param.fid, '<g id="topo-shadow" inkscape:groupmode="layer" inkscape:label="topo-shadow" clip-path="url(#IDmain)">\n');
        for j=length(Cld):-1:1
            if strcmp(get(Cld(j), 'Type'), 'patch')
                Iso = get(Cld(j), 'CData');
                if ~isnan(Iso)
                    wmtp_plotpath(wmt_param.fid, Cld(j), wmt_param.pictsize, ['class="shadow' int2str(find(shadowlims==Iso)) '"']);
                end
            end
        end
        fprintf(wmt_param.fid, '</g>\n');
    end

    %# hold on;

    %#%# raster format
else

    %#first cut data according to limits
    y1 = length(find (y<=wmt_param.limits(1)))+1;
    y2 = length(find (y<=wmt_param.limits(2)));
    x1 = length(find (x<=wmt_param.limits(3)))+1;
    x2 = length(find (x<=wmt_param.limits(4)));
    data = data(y1:y2, x1:x2);
    data = flipud(data);

    %# now plot the relief with shadows
    mapdata = ones(size(data));
    mapdatarel = zeros(size(mapdata, 1)-1, size(mapdata, 2)-1, 3);

    %#define topo levels
    lngth = size(levels, 2);
    for lv = 2:lngth;
        mapdata(find((data>=levels(lv-1)) & (data<levels(lv)))) = lv;
    end
    mapdata(data <= levels(1)) = 1;
    mapdata(data >= levels(lngth)) = lngth;

    %#mapdata(find(mapdata==0)) = nan;

    %# define relief
    [r, g, b] = ind2rgb(mapdata, cmap);
    r1 = r(2:end, 2:end);
    g1 = g(2:end, 2:end);
    b1 = b(2:end, 2:end);

    if ~isempty(shadowlims)
        dff = data(2:end, 2:end)-data(1:end-1, 1:end-1);
        dff(dff>=0) = 0;
        dff(dff<=-350) = -350;
        dff = ones(size(dff)) - (dff / min(min(dff)));
        r1 = dff.*r1;
        g1 = dff.*g1;
        b1 = dff.*b1;
    end

    mapdatarel(:,:,1) = r1;
    mapdatarel(:,:,2) = g1;
    mapdatarel(:,:,3) = b1;

    filename = [tempname() '.jpg'];
    imwrite(mapdatarel, filename, 'Quality', 100);

    pictsize = size(mapdatarel);

    %# print image
    fprintf(wmt_param.fid,'<g id="topo" inkscape:groupmode="layer" inkscape:label="topo" clip-path="url(#IDmain)">\n');
    fprintf(wmt_param.fid,'<image x="0" y="0" width="%d" height="%d"  transform="scale(%0.5f %0.5f)" ', pictsize(2), pictsize(1), wmt_param.pictsize(1)/pictsize(2), wmt_param.pictsize(2)/pictsize(1));
    %#fprintf(wmt_param.fid,'<image x="0" y="0" width="%d" height="%d" ', wmt_param.pictsize(2), wmt_param.pictsize(1));
    fprintf(wmt_param.fid,'xlink:href="data:image/jpeg;base64,');
    ifid = fopen(filename);
    while ~feof(ifid)
        %#     fprintf(wmt_param.fid, '%s\n', wmtp_base64encode(fread(ifid,57), '')); %fread(ifid,57)does not work with wikmedia rendering
        fprintf(wmt_param.fid, '%s\n', wmtp_base64encode(fread(ifid), ''));
    end
    fclose(ifid);

    fprintf(wmt_param.fid,'"/>\n');
    fprintf(wmt_param.fid,'</g>\n');
end