get_tif.m 1.7 KB
function get_tif(filename,I,type)
% GET_TIF - genereates geotif files 
%
% Synthax: get_tif(filename,I,type)
% 
% Inputs:
%    filename - output filename (.tif)
%    I        - georectified image  
%    type     - 'dd' for decimal degre / 'utm' for UTM
%
%
% Example:
%    get_tif('201802029_2141_4m_UTM.tif',I,'utm')
%
% Author: Paul Nicot
% UQAR/ISMER
% March 2018
%______________________________________________________________________


rlon = evalin('base','rlon');
llon = evalin('base','llon');
llat = evalin('base','llat');
ulat = evalin('base','ulat');

if strcmp(type,'dd');
    R = georasterref();
    R.RasterSize = [size(I)];
    R.Latlim = [llat ulat];
    R.Lonlim = [llon rlon];
    R.ColumnsStartFrom = 'north';
    geotiffwrite(filename, I, R);

elseif strcmp(type,'utm');
    %% create geotiff file
    R = georasterref();
    R.RasterSize = [size(I)];
    R.Latlim = [llat ulat];
    R.Lonlim = [llon rlon];
    R.ColumnsStartFrom = 'north';

    %% utm
    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);
end