Giorgos Kourakos | Subsurface

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

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

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

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

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

up=~up;

end

outraster(inraster) = id(k)*1000+tempraster(inraster);

end

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

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

Subsurface

George Kourakos

Author: Giorgos Kourakos

Greek Bathing Water Profiles

Equal_Areas.m

obj_fun.m