Uncategorized posts

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

Lets suppose that two shapes describe the same geographic area each digitized by different people at different coordinate system e.g. different transform, rotation and scaling and you want to match the two shapes.

This may seem unusual task. Unfortunately I had to deal with it. In particular I was given a Modflow model which was rotated and translated so that the left bottom corner of the domain was the (0,0) coordinate and I had to convert it to the coordinate system I was using. Also the two models described the same area, yet the digitization of the outline was made by different people therefore there were significantly different in the details.

The approach I took might not be the best and I welcome suggestions.

As I have some background in optimization I decided to run an optimization problem where the decision variables would be the translate, rotate and scale factors and the objective function would be the minimization of the area that remains after a Boolean subtraction between the two shapes.

Although this may sound a bit complicated task, it turned out to be relatively trivial using Matlab as it provides all the functionality one may need.

Lets name the outline of the modflow domain, modflow shape and the one with actual coordinates the target shape.

First we need to have the two shapes into Matlab workspace. For the modflow shape I rasterized the modflow grid and convert it to polygon shapefile in QGIS. The target shape was already a shapefile. The two shape files can be easily loaded with the shaperead matlab function:

modflow = shaperead('CVHM_outline');
target = shaperead('B118_outline_simple');
subplot(1,2,1); plot( modflow(1,1).X, modflow(1,1).Y )
subplot(1,2,1); title('Modflow')
axis equal
subplot(1,2,2); plot( target(1,1).X, target(1,1).Y )
subplot(1,2,2); title('Target')
axis equal

modflow = shaperead('CVHM_outline');

target = shaperead('B118_outline_simple');

subplot(1,2,1); plot( modflow(1,1).X, modflow(1,1).Y )

subplot(1,2,1); title('Modflow')

axis equal

subplot(1,2,2); plot( target(1,1).X, target(1,1).Y )

subplot(1,2,2); title('Target')

axis equal

It is quite important to have the two shapefiles somehow simplified otherwise the boolean operations might slow down the process.

Next we have to write the objective function. The arguments will be a vector 4×1 (X, Y, R S) i.e. x,y translation, rotation and scaling.

Create a new file, e.g. obj_fun.m and paste the code from this link:

Now we are ready to run the optimization problem. I have found that the pattern search is very reliable optimization algorithm for this problem.

First I run the optimization without scaling because small changes can make big differences. This will give an idea about the scaling limits.

 % Set the optimization options first
pa_opt = psoptimset('CompletePoll', 'on', 'CompleteSearch','on', ...
                    'MaxIter', 1000, 'PlotFcns',{@psplotbestf, @psplotbestx},...
                    'PollMethod', 'GPSPositiveBasis2N', 'SearchMethod',@GPSPositiveBasis2N);
 % run the optimization
[x, fval] = patternsearch(@obj_fun,[0 0 0], [], [], [], [], [],[], [], pa_opt);

% Set the optimization options first

pa_opt = psoptimset('CompletePoll', 'on', 'CompleteSearch','on', ...

'MaxIter', 1000, 'PlotFcns',{@psplotbestf, @psplotbestx},...

'PollMethod', 'GPSPositiveBasis2N', 'SearchMethod',@GPSPositiveBasis2N);

% run the optimization

[x, fval] = patternsearch(@obj_fun,[0 0 0], [], [], [], [], [],[], [], pa_opt);

The figure above shows the result of the first optimization. The red shape is the translated and rotated geometry. The translation and rotation looks good. But definitely we need scaling to make it better. We can also estimate that the scaling factor has to be greater than 1 and less than 1.5. Therefore we will set those values as constrains. In addition we will set the starting point of the next optimization run the solution of the previous run.

[xf, fval] = patternsearch(@obj_fun,[xf(1) xf(2)  xf(3) 1.1], [], [], [], [], ...
                           [xf(1)-1e5 xf(2)-1e5 30 1.0], ...
                           [xf(1)+1e5 xf(2)+1e5 32 1.5], [], pa_opt);

[xf, fval] = patternsearch(@obj_fun,[xf(1) xf(2) xf(3) 1.1], [], [], [], [], ...

[xf(1)-1e5 xf(2)-1e5 30 1.0], ...

[xf(1)+1e5 xf(2)+1e5 32 1.5], [], pa_opt);

The figure below shows the final results of the optimization run

I should note that I had to perform quite a few optimization runs to reach this solution. In addition I performed few runs with fewer desicion variables. For example keep the rotation constant and optimize for scaling etc.

As you can see the result is not perfect and there are discrepancies especially in the north area. Therefore it may need to optimize for rotation.

Subsurface

George Kourakos

Equal_Areas.m

Uncategorized Posts

How to align 2D shapes