function asc_cleanup_and_plot_difference(polygon, fname1, fname2, varargin)

OPT = struct(...
    'fieldfun', @lower);
OPT = setproperty(OPT, varargin);

if fname2 == false
    fname2 = fname1;
end

% read both files
[z1, x1, y1] = arcascread(fname1);
[z2, x2, y2] = arcascread(fname2);
% flip in y direction (because y seems to be positive south)
z1 = flipud(z1);
z2 = flipud(z2);

% find the corresponding part of the grids
[X, ix1, ix2] = intersect(x1, x2);
[Y, iy1, iy2] = intersect(y1, y2);

if any(cellfun(@isempty, {X, Y}))
    disp('no overlap found')
    return
end

% calculate the z difference
zdiff = z1(iy1,ix1) - z2(iy2,ix2);


[x, y] = cellfun(@(pg) deal(pg.(feval(OPT.fieldfun, 'x')), pg.(feval(OPT.fieldfun, 'y'))), polygon,...
    'uniformoutput', false);

figure
subplot(1,2,1)
title('original')
pcolor(X, Y, zdiff)
shading interp
hold on
for ipol = 1:length(polygon)
    % plot all polygons
    plot(x{ipol}, y{ipol})
end
axis image

[X,Y] = meshgrid(X,Y);
IN = false(size(X));
for ipol = 1:length(x)
    % derive the grid points that are in any of the polygons
    IN = IN | inpolygon(X,Y, x{ipol}, y{ipol});
end
% get the overlapping part of the duplicate grid
z1part = z1(iy1,ix1);
% derive which of the points must become nan
nanidx = IN & zdiff==0;
% set those points to nan
z1part(nanidx) = nan;
% allocate the new z
z1new = z1;
% modify the z values of concern
z1new(iy1,ix1) = z1part;

% plot the modified z difference
subplot(1,2,2)
title('post-cleanup')
pcolor(X, Y, z1new(iy1,ix1) - z2(iy2,ix2))
shading interp
hold on
for ipol = 1:length(x)
    plot(x{ipol}, y{ipol})
end
axis image


fprintf('%i values set to NaN\n', sum(nanidx(:)));
if all(isnan(z1new))
    % delete the file if all values are nan (after modication)
    fprintf('"%s" now contains nans only\n', fname1)
    delete(fname1);
    fprintf(2, '"%s" deleted\n', fname1)
elseif sum(nanidx(:)) == 0
    % do nothing if no change has to be made
    fprintf(1, '"%s" kept as is\n', fname1)
else
    % overwrite the file if part of it has changed
    [X1,Y1] = meshgrid(x1,y1);
    arcgridwrite(fname1,X1,Y1,z1new)
    fprintf(2, '"%s" modified\n', fname1)
end

function [z, x, y] = arcascread(fname)
    OPT.block_size          = 1e6;
    OPT.z_scalefactor       = 1; %scale factor of z values to metres altitude

fid = fopen(fname);

s = textscan(fid,'%s %f',6);

ncols        = s{2}(strcmpi(s{1},'ncols'       ));
nrows        = s{2}(strcmpi(s{1},'nrows'       ));
xllcorner    = s{2}(strcmpi(s{1},'xllcorner'   ));
yllcorner    = s{2}(strcmpi(s{1},'yllcorner'   ));
cellsize     = s{2}(strcmpi(s{1},'cellsize'    ));
nodata_value = s{2}(strcmpi(s{1},'nodata_value'));
if isempty(ncols)||isempty(nrows)||isempty(xllcorner)||isempty(yllcorner)||isempty(cellsize)||isempty(nodata_value)
    error('reading asc file')
end

%% read file chunkwise
WB.tic = tic;
small_number  = 1e-16;

kk = 0;
while ~feof(fid)
    % read the file
    kk       = kk+1;
    D{kk}    = textscan(fid,'%f64',floor(OPT.block_size/ncols)*ncols,'CollectOutput',true); %#ok<AGROW>
    D{kk}{1} = reshape(D{kk}{1},ncols,[])'; %#ok<AGROW>
    if all(abs(D{kk}{1}(:) - nodata_value) < small_number)
        D{kk}{1} = nan; %#ok<AGROW>
    else
        D{kk}{1}(abs(D{kk}{1}-nodata_value) < small_number) = nan; %#ok<AGROW>
    end
end
fclose(fid);
if length(D) == 1
    z = D{1}{1};
else
    D2 = cellfun(@(x) x{1}, D, 'UniformOutput', false);
    z = cat(1, D2{:});
end
x = (xllcorner + 10) + (0:(ncols-1))*cellsize;
y = (yllcorner + 10) + (0:(nrows-1))*cellsize;

