Greek Bathing Water Profiles

The quality of bathing waters in Greece is systematically monitored since 1988 according to Directive 76/160/EEC “on the quality of bathing waters” under a Programme organized and coordinated by the MEECC (Ministry of Environment and Energy). Since 2010, the quality of bathing water is monitored in accordance with the new Directive 2006/7/EC “on the management of bathing water quality” as transposed to Greek legislation with JMD 8600/416/E103/2009 (GG 356V/2009) under the “Country’s Monitoring Programme for bathing water quality” (hereinafter called the “Programme”).

This is the google fusion version

This is the Google Map version

Equal_Areas.m

function out_raster = Equal_Areas(in_raster, Nunits, id, excludeid, method, opt)
    % out_raster = Equal_Areas(in_raster, Nunits, id, excludeid, method, opt)
    % This function divides a raster into equal areas
    %
    % INPUTS:
    % in_raster     : raster to be divided. Each cell id in the raster is
    %                 treated as a different region
    % Nunits        : The number of cells that each each equal area cell will have (Units are cells).
    % id            : Divides only the area associated with the particular id.
    %                 If empty then all areas are divided. 
    %
    % excludeid     : any area associated with ids in this list will not be divided
    %                   This is usefull to exclude no data cell.
    % method        : switches between the scott and metis method
    % opt           : options for the metis case. Consists of the following fields  
    %          name : name of the simulation.
    %    metis_path : Path of the metis program. It can be empty if
    %                 Default_metis_path variable is defined
    %
    % OUTPUTS
    % outraster     : The raster with the equal area cells.


switch method
    case 'scott'
        out_raster = Equal_Area_Scott(in_raster, Nunits, id, excludeid);
        
    case 'metis'
        All_ids = unique(in_raster);
        Xcoords = 1:size(in_raster,2);
        Ycoords = size(in_raster,1):-1:1;
        out_raster = nan(size(in_raster));
        for i = 1:length(All_ids)
            if ~isempty(find(excludeid == All_ids(i),1))
                continue
            end
            N_pixel_per_region = sum(sum(in_raster == All_ids(i)));
            N_eq_areas = round(N_pixel_per_region / Nunits);
            clear p TR;
            [ii, jj] = find(in_raster == All_ids(i));
            p(:,1) = Xcoords(jj);
            p(:,2) = Ycoords(ii);
            TR = delaunay( p(:, 1), p(:, 2) );
            %run metis
            if N_eq_areas < 2
                N_eq_areas = 2;
            end
            opt = check_metis_options(opt);
            [ElPart, NDPart] = decompose_2D_mesh(TR, N_eq_areas, opt);
            
            out_raster(sub2ind(size(out_raster), ii, jj)) = All_ids(i)*1000+NDPart+1;
        end
        
    otherwise
        warning('Uknown method. Use scott or metis');
        out_raster = nan;
        
end

end


function [ElPart NDPart]=decompose_2D_mesh(MSH_2d,Nsub,opt)
    % This function decompose a 2D mesh into Nsub subdomains
    %
    %Input:
    % p_2d :[np x 2] coordinates of mesh points
    %
    % MSH_2d : Mesh structure
    %
    % Nsub : either a number or a string. In case of string specifies the file
    % that decribes the subdomain division. Otherwise runs the Metis and create
    % the file.
    %
    % Nlay is the number of layers. The elevation has to be adjusted later

    % opt structure variable with fields 
    %    name : name of the simulation. In case number Nsub, this is the name 
    %    metis_path : Path of the metis program. It can be empty if
    %                 Default_metis_path variable is defined
    %

    %create inputfile for mpmetis
frmt='%g';
for ii=2:size(MSH_2d,2)-1
    frmt=[frmt ' %g'];
end
frmt=[frmt ' %g\n'];
fid=fopen([opt.name '.mesh'],'w');
fprintf(fid,'%g\n',size(MSH_2d,1));
fprintf(fid,frmt,MSH_2d');
fclose(fid);

    %run metis
system([opt.metis_path ' ' opt.name '.mesh ' num2str(Nsub)])

    %read mesh partition
fid=fopen([opt.name '.mesh' '.epart.' num2str(Nsub)]);
A=textscan(fid,'%f');
fclose(fid);
ElPart=A{1,1};

fid=fopen([opt.name '.mesh' '.npart.' num2str(Nsub)]);
A=textscan(fid,'%f');
fclose(fid);
NDPart=A{1,1};

end

function opt = check_metis_options(opt)
    if ispc
        Default_metis_path='c:\Users\gkourakos\Downloads\metis-5.0.2\build\windows\programs\Release\mpmetis.exe';
    else
        Default_metis_path='/home/giorgos/Downloads/metis-5.1.0/build/Linux-x86_64/programs/mpmetis';
    end
    if isempty(opt)
        opt.metis_path=Default_metis_path;
        opt.name='tempMetis';
    end
    if isempty(opt.metis_path); 
        opt.metis_path=Default_metis_path; 
    end

end


function outraster = Equal_Area_Scott(raster, min_area, id, excludeid)
    % outraster = Equal_Area(inraster, min_area, Neq_areas, id)
    % This function divides a raster into equal area cells
    %
    % INPUTS:
    % raster      : raster to be devided. Each cell id in the raster is
    %                 treated as a different region
    % min_area      : The area of each equal area cell (Units are cells).
    % id            : Divides only the area associated with the particular id.
    %                 If empty then all areas are divided. 
    %
    % excludeid     : any area associated with ids in this list will not be divided
    %                   This is usefull to exclude no data cell.
    %
    % OUTPUTS
    % outraster     : The raster with the equal area cells.
    % 
    % Usage:
    % If the raster comes from ARCGIS layer then read it into the Matlab workspace using
    % [raster metadata]=readArcGisASCIIfile(asciiraster.asc);
    %
    % Next this function to create the equal area cells 
    %
    % and finaly convert them to ARCGIS raster ascii.asc file using 
    % createAsciiforRaster(filename,TAB,xl,yl,csz,nodata)
    % (The additional functions are part of the mSim software http://groundwater.ucdavis.edu/msim/)

    if isempty(id)
        id = unique(raster);
    end
    id = setdiff(id, excludeid);
    [nr nc] = size(raster);
    outraster = nan(nr,nc);

    for k = 1:length(id)
        inraster=raster==id(k);
        NumCells=sum(sum(inraster));
        NumCell_per_eq_area=round(min_area);
        Ncat=round(NumCells/NumCell_per_eq_area);
        tempraster=zeros(size(inraster));
        tempraster(~inraster)=nan;

        Strips=zeros(size(inraster));
        Nstrips=round(sqrt(Ncat));

        NumCell_strip=round(NumCells/Nstrips);
        cumsum_strips=cumsum(sum(inraster));
        istart=1;
        for i=1:Nstrips-1
            [~,iend]=min(abs(cumsum_strips-i*NumCell_strip));
            Strips(:,istart:iend)=i;
            istart=iend+1;
        end
        Strips(:,istart:end)=Nstrips;

        iCat=1;cntCells=0;
        up=true;
        for i=1:Nstrips
            fprintf('%g strip out of %g\n', [i Nstrips])
            %find the lower left unused pixel
            [Is Js]=ind2sub(size(Strips),find(Strips==i,1));
            if i~=Nstrips
                [Ie Je]=ind2sub(size(Strips),find(Strips==i+1,1));
            else
                Je=nc;
            end
            if up
                Is=1;Ie=nr;step=1;
            else
                Is=nr;Ie=1;step=-1;
            end
            for ii=Is:step:Ie
                for jj=Js:Je-1
                    if ~isnan(tempraster(ii,jj)) && tempraster(ii,jj)==0
                        tempraster(ii,jj)=iCat;
                        cntCells=cntCells+1;
                        if cntCells==NumCell_per_eq_area
                            iCat=iCat+1;
                            if iCat>Ncat
                                iCat=Ncat;
                            end
                            cntCells=0;
                        end
                    end
                end
            end
            up=~up;
        end
        outraster(inraster) = id(k)*1000+tempraster(inraster);
    end
end

 

obj_fun.m

function f = obj_fun(X)
 % Although matlab provides better ways to overload the objective function
 % with additional data, I'll keep things simple here and use the following
 % trick to bring the two shapes into the function workspace
target = evalin('base', 'target');
c2vsim = evalin('base', 'c2vsim');

 % Because the shapefiles contain islands I split them into polygons and use
 % only the first which is the outline.
[X1, Y1] = polysplit(c2vsim.X, c2vsim.Y);
[X2, Y2] = polysplit(target.X, target.Y);


DX = X(1); % This is the x translation
DY = X(2); % This is the y translation
 % first I run the optimzation without rotation and scaling
 % So I give here some default values
if length(X) > 2
    TH = deg2rad(X(3)); % This is the rotation
else
    TH = 0;
end
if length(X) > 3
    S = X(4);
else
    S = 1;
end

 % translate shape to origin
xc = mean(X1{1,1});
yc = mean(Y1{1,1});
X1{1,1} = X1{1,1} - xc;
Y1{1,1} = Y1{1,1} - yc;

 % first scale
X1{1,1} = S*X1{1,1};
Y1{1,1} = S*Y1{1,1};

 % then rotate around the origin (0,0)
 % in case the shape is away from the origin it may be necessary to
 % translate first so that the centroid of the shape is on the origin,
 % then rotate and translate back again
XR = X1{1,1}*cos(TH) -  Y1{1,1}*sin(TH);
YR = X1{1,1}*sin(TH) +  Y1{1,1}*cos(TH);

XR = XR + xc;
YR = YR + yc;

 % and then translate according to the desicion variables
XR = XR + DX;
YR = YR + DY;


 % find their intersection
[x, y] = polybool('intersection', XR, YR, X2{1,1}, Y2{1,1});

if isempty(x)
    % if the two shapes do not intersect they are still far away therefore
    % return the distance between their barycenters plus a constant that
    % will make the objective function value greater than any possible
    % value if they are intersected.
    f = sqrt((mean(XR) - mean(X2{1,1}))^2 + (mean(YR) - mean(Y2{1,1}))^2) + 1000;
else
    % If they intersect compute the area that remains after we subtract the
    % two shapes. We want this area minimum. Note that we have divided the
    % area by a 10000 sq km. to make this value small.
    % We perform two subtraction operators, to avoid the case where one is
    % completely inside the other and the subtracted area is zero.
    [x1, y1] = polybool('subtraction', X2{1,1}, Y2{1,1}, XR, YR);
    [x2, y2] = polybool('subtraction', XR, YR, X2{1,1}, Y2{1,1});

    [xs, ys] = polysplit(x1, y1);
    f1 = 0;
    for i = 1:size(xs,1)
        f1 = f1 + polyarea(xs{i,1}, ys{i,1})/(10000*10000);
    end

    [xs, ys] = polysplit(x2, y2);
    f2 = 0;
    for i = 1:size(xs,1)
        f2 = f2 + polyarea(xs{i,1}, ys{i,1})/(10000*10000);
    end
    % we return the maximum subtracted area.
    f = max(f1, f2);
end