r_interp_opc_UTM.m 8.19 KB
function r_interp_opc_UTM(expNumber,pretreat,zone,resol,iFirst,nSkip,iEnd)
% R_INTERP_OPC_UTM - Uses the rectification data to interpolate RGB values on a
% regular grid (UTM) and save it to an image. It is specifically adapted to
% images taken from the Observatoire du Pic Champlain (OPC)
% 
% Syntax:  r_interp_opc_UTM(expNumber,pretreat,zone,resol,iFirst,nSkip,iEnd))
% 
% Inputs:
%    expNumber - The date number (yyyymmdd). Points to a directory 
%    pretreat  - 0: distorted images, 1: undistorted images, 2: stabilized,
%                3: undist + stable
%    zone      - The part of the scene to analyse. See the code for details
%    resol     - image resolution of the georectification (in meter)
%    iFirst    - The first image to be treated
%    nSkip     - The number of images skipped
%    iEnd      - The last image to be treated
%                (if iEnd is set to zero, the loop is done from iFirst to
%                 nFiles)
% 
% Example:
%    r_interp_opc_UTM(20140203,1,3,1.5,1,1,0)
% 
% Other m-files required: g_rect toolbox
%   (http://demeter.uqar.ca/g_rect)
% Subfunctions: none
% MAT-files required: the rectification parameter file
% 
% See also: G_RECT, R_UNDISTORT, R_ROTATE

% Author: Dany Dumont/Paul Nicot
% UQAR/ISMER
% email: dany_dumont@uqar.ca
% Website: http://www.ismer.ca/dumont-dany
% October 2017
% ______________________________________________________________________

estr = num2str(expNumber);
zstr = num2str(zone);

rootdir = '/share/archives/dumoda01/data/photogrammetry/opc';
expdir  = [rootdir,'/',estr(1:4),'/',estr];
if pretreat == 0;
    origdir = [expdir,'/original'];
    intdir  = [expdir,'/interp'];
    calfile = [rootdir,'/',estr(1:4),'/calib/opc_rect_' estr '_UTM.mat'];
elseif pretreat == 1;
    origdir = [expdir,'/original_undist'];
    intdir  = [expdir,'/interp_undist'];
    calfile = [rootdir,'/',estr(1:4),'/calib/opc_rect_' estr '_dist_UTM.mat'];
elseif pretreat == 2;
    origdir = [expdir,'/original_stable'];
    intdir  = [expdir,'/interp_stable'];
    calfile = [rootdir,'/',estr(1:4),'/calib/opc_rect_' estr '_stable_UTM.mat'];
elseif pretreat == 3;
    origdir = [expdir,'/undist_stable'];
    intdir  = [expdir,'/interp_undist_stable'];
    calfile = [rootdir,'/',estr(1:4),'/calib/opc_rect_' estr '_dist_stable_UTM.mat'];
end
zonedir = [intdir,'/zone',zstr];
tifdir  = [zonedir,'/tif'];



%% Load the calibration file
load(calfile);
% files for context image
S = shaperead('/home/nicp0003/LND_LSLE_line_UTM.shp');
[res,Rr] = geotiffread([rootdir '/' estr(1:4) '/calib/resolution_' estr(1:4) '_UTM.tif']);
res = double(res);

%% Creation of the grid

if zone == 1      % chenal du bic
    llon = 508152;  rlon = 514063;
    llat = 5353671; ulat = 5360778;

elseif zone == 2  % entree BHH + au large
    llon = 510208;  rlon = 512933;
    llat = 5353902; ulat = 5358298;

elseif zone == 3  % entree BHH
    llon = 510208;  rlon = 512939;
    llat = 5353902; ulat = 5355976; 

% elseif zone == 4  % mouillage AWAC
%     llon = 508522;  rlon = 509631;
%     llat = 5354097; ulat = 5355766;

elseif zone == 5
    llon = 512050;  rlon = 513254;
    llat = 5353830; ulat = 5355470;
        
elseif zone == 6 % entree + fond BHH
    llon = 512050;  rlon = 514500;
    llat = 5353830; ulat = 5355800;
    
elseif zone == 7 % BHH complet
    llon = 510000;  rlon = 514400;
    llat = 5353600; ulat = 5355800;
end


% update grid with resolution
%str_resol = num2str(resol*100);
%lon1 = llon ; lon2 = rlon;
x1vec =  llon:resol:rlon;
y1vec =  llat:resol:ulat;
y1vec = y1vec';

[X1,Y1] = meshgrid(x1vec,y1vec);

cd(origdir);

imfiles = dir('*.jpg');
nFiles  = size(imfiles,1);

%% Masking
mask = ones(5472,3648);
mask(:,   1:1300) = NaN;  
mask(:,3000:end ) = NaN;

LON = mask.*LON; %#ok<NODEF>
LAT = mask.*LAT; %#ok<NODEF>

x1 = LON(:); x2 = LON(:);
y1 = LAT(:);

if ~exist(intdir,'dir')
    mkdir(intdir)
end
if ~exist(zonedir,'dir')
    mkdir(zonedir)
end
if ~exist(tifdir,'dir')
    mkdir(tifdir)
end

save([zonedir,'/zone',zstr,'.mat'],'X1','Y1');

if iEnd == 0
    iEnd = nFiles;
end

%% plot context image
h = figure; set(h,'visible','off');

lev1 = [2:2:10];
mapshow(res, Rr,'DisplayType', 'contour','LineStyle','--', ...
    'LevelList',lev1,'edgecolor', 'none','LineColor', 'black',...
    'ShowText', 'on');

hold on;
lev2 = [15,25,50,75,100];
mapshow(res, Rr,'DisplayType', 'contour','LineStyle','--', ...
    'LevelList',lev2,'edgecolor', 'none','LineColor', 'black', ...
    'ShowText', 'on');

mapshow(S,'color',[0.3 0.3 0.3],'LineWidth',1);

ylim([5353000, 5363000]);
xlim([506000, 515000]);

line([rlon llon],[llat llat],'Color',[1 0 0]);
line([llon rlon],[ulat ulat],'Color',[1 0 0]);
line([rlon rlon],[llat ulat],'Color',[1 0 0]);
line([llon llon],[llat ulat],'Color',[1 0 0]);

iname = [zonedir,'/georect_context_',estr,'.png'];
saveas(h,iname);


%% georect loop
for n = iFirst:nSkip:iEnd
    
    imfilename = imfiles(n).name;
    imNumber   = imfilename(1:13);
    
    I0 = imread(imfilename);
    I0 = double(I0);
 
    if pretreat == 2
            K0 = mask.*I0';
        else
            K0 = mask.*mean(I0,3)';
    end

 
    k0 = K0(:);
    
    k0(isnan(x2)) = [];
    x1(isnan(x1)) = [];
    y1(isnan(y1)) = [];
    
    Fk = TriScatteredInterp(x1,y1,k0);
    
    K = flipud(Fk(X1,Y1));
    
    I = uint8(K);
    
    if pretreat == 0;
            filename = [zonedir,'/',imNumber,'_',zstr,'_UTM.jpg'];
        elseif pretreat == 1;
            filename = [zonedir,'/',imNumber,'_',zstr,'_undist_UTM.jpg'];
        elseif pretreat == 2;
            filename = [zonedir,'/',imNumber,'_',zstr,'_stable_UTM.jpg'];
        elseif pretreat == 3;
            filename = [zonedir,'/',imNumber,'_',zstr,'_undist_stable_UTM.jpg'];
    end
 
    
    imwrite(I,filename,'jpg');

    
    % create geotiff file UTM
    if pretreat == 0;
            filename = [zonedir,'/tif/',imNumber,'_',zstr,'_UTM.tif'];
        elseif pretreat == 1;
            filename = [zonedir,'/tif/',imNumber,'_',zstr,'_undist_UTM.tif'];
        elseif pretreat == 2;
            filename = [zonedir,'/tif/',imNumber,'_',zstr,'_stable_UTM.tif'];
        elseif pretreat == 3;
            filename = [zonedir,'/tif/',imNumber,'_',zstr,'_undist_stable_UTM.tif'];
    end

    mstruct = defaultm('utm');
    mstruct.geoid = referenceEllipsoid('WGS84','meters');
    mstruct.zone = num2str(19);
    mstruct = defaultm(mstruct);

    R = maprasterref( ...
         'RasterSize', size(I), ...
         'XWorldLimits', [llon rlon], ...
         'YWorldLimits', [llat ulat], ...
         'ColumnsStartFrom', 'north', ...
         'RasterInterpretation', 'postings');

    key = struct( ...
        'GTModelTypeGeoKey',[], ...
        'GTRasterTypeGeoKey',[], ...
        'ProjectedCSTypeGeoKey',[]);

    key.GTModelTypeGeoKey = 1;    % for Projected Coordinate System
    key.GTRasterTypeGeoKey = 2;   % for PixelIsPoint (postings)
    key.ProjectedCSTypeGeoKey = 32619;  % for UTM 19N
    
    geotiffwrite(filename,I,R,'GeoKeyDirectoryTag',key);

    disp([' ',imNumber]); disp(datetime);

end

%% write a README.txt with the georectification informations
fidname = [zonedir,'/README.txt'];
fid = fopen(fidname,'w');

     fprintf(fid,'%s\n','README.txt');
     fprintf(fid,'%s\n','This file contains informations about the georectification performed');
     fprintf(fid,'%s\n','with the matlab script r_interp_opc_UTM.m');
     fprintf(fid,'%s\n','----------------------');
     fprintf(fid,'%s\n',['Date of processing : ' char(datetime) ]);
     fprintf(fid,'%s\n',['Lower Left corner  : ' num2str(llat) ' / ' num2str(llon) ]);
     fprintf(fid,'%s\n',['Upper right corner  : ' num2str(ulat) ' / ' num2str(rlon) ]);
     fprintf(fid,'%s\n',['Resolution in x and y : ' num2str(resol) ' m/pixel' ]);
     fprintf(fid,'%s\n',['Georectified directory : ' char(intdir) ]); 
     fprintf(fid,'%s\n',['First georectified image : ' imfiles(1).name ]);
     fprintf(fid,'%s\n',['Last georectified image : ' imfiles(end).name ]);
     fprintf(fid,'%s\n',['georectified matrix file : ' calfile ]);
     fprintf(fid,'%s\n','----------------------');

quit
%end