image thumbnail

ImportAsciiRaster

by

 

12 Mar 2008 (Updated )

Import Arc/Info AsciiRaster

ImportAsciiRaster(varargin)
function [Z R] = ImportAsciiRaster(varargin)
% [Z R] = ImportAsciiRaster(...)
% 
%     % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%     %                                                                 %
%     %                 Produced by Giuliano Langella                   %
%     %                   e-mail:gyuliano@libero.it                     %
%     %                           March 2008                            %
%     %                                                                 %
%     %                 Last Updated: 20 April, 2010                    %
%     %                                                                 %
%     % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%
%
%
%   -----   TODO   -----
%   1.\ The function should take care of xll-center and yll-center when
%       they are found instead of -corner conditions.
% 
%   -----  SYNTAX  -----
%   Z = ImportAsciiRaster;                  ==> nodata:NaN; ref:'h'; type:'s'; FilePath:uigetfile
%   [Z h] = ImportAsciiRaster;              ==> nodata:NaN; ref:'h'; type:'s'; FilePath:uigetfile
%   [Z R] = ImportAsciiRaster(varargin);    ==> see OPTIONS
%
%
%   -----  OPTIONS  -----
%  +'nodata'(optional) [numeric] = output No-Data Value.
%   Default output NoData is NaN. When specified it substitutes the input
%   no-data value, which is stored in 'NODATA_value' from header.
%
%  +'ref'(optional) [string]: Geospatial Reference Matrix.
%   --> 'h' (default) for header matrix (as the .hdr file in Arc-Info Binary Raster).
%   --> 'r' for spatial-referencing matrix R (Mapping Toolbox);
%
%  +'type'(optional) [string]: class type of Z output in multiselection case only.
%   --> 's' (default) gets a structure array with field name equal to
%       imported filename without extension;
%   --> 'd' gets a double array with size (nrows, ncols, N. files); it is
%       useful in case of several layers with same reference matrix.
%       NOTE: be careful to multi-select grids with same spatial extent.
%
%  +'FilePath'(optional) [char]: Complete path + filename + file extension.
%       If it is not given the uigetfile will ask for file(s).
%
%   -----  DESCRIPTION  -----
% This function imports Arc-Info ASCII Raster (*.asc, *.txt) files. You are
% allowed to import single or multiple files by multi-selection (CRTL +
% Mouse_Left_Click). In the latter case by default a structure array is
% created in order to store the several rasters; items of the structure
% array are called as the imported raster. A facultative option is now
% added to create a stack of double 3-D array instead of a structure array
% in case more grid files are selected. This might be useful to handle for
% instance multi-temporal remote sensing datasets (such as NDVI), or
% multiple realizations of Geostatistical Simulation. You can create a
% stack with dimension compatibly with you computer memory.
% 
%   -----  EXAMPLE_1  -----
%   [Z R] = ImportAsciiRaster('r');           % import DEM and Aspect
%   Exaggeration_Factor = 1.6;                % use an exaggeration factor
%   figure(gcf); mapshow(Z.DEM*Exaggeration_Factor,R.DEM,'DisplayType','surface');
%   colormap(demcmap(Z.DEM))
%   view(3); grid off; axis off
%   axis on; xlabel('Easting'); ylabel('Northing'); zlabel('Elevation')
%
%   -----  EXAMPLE_2  -----
%   [Z R] = ImportAsciiRaster(NaN, 'r', 'd');       %import 15 Geostatistical Simulations
%   M = mean(Z,3);      %compute mean along 3rd dimension
%   figure(gcf); h=pcolor(M); axis off;             %pcolor plot
%   set(h,'LineStyle','none')
%   colormap(jet), colorbar('Location', 'SouthOutside')
%   saveas(gcf, 'clay_simulation.eps', 'psc2')      %save to color EPS for LaTex
%

%% MAIN

% Read varargins [given by user]
FilePath = cell(0);
FileName = cell(0);
for i = 1:size(varargin,2)
    if iscell(varargin{i})
        FilePath = varargin{i};
        for i = 1:length(FilePath)
            FileName{end+1} = path2file(FilePath{i});
        end
    else
        switch varargin{i}
            case{'h','r'}                                       % header
                ref = varargin{i};
            case{'s','d'}                                       % var type
                type = varargin{i};
            otherwise        
                if isnan(varargin{i}) | isnumeric(varargin{i})  % nodata
                    nodata = varargin{i};
                elseif ischar(varargin{i})                      % filepath
                    FilePath{end+1} = varargin{i};
                    FileName{end+1} = path2file(FilePath{end});
                    if ~exist(FilePath{end},'file')
                        error(['File ',varargin{i},' does not exist!'])
                    end
                end
        end
    end
end

% get file(s) [use UI if any file is given as varargin]
if ~sum(size(FilePath))
    % UIGETFILE - .asc/.txt - Multiselect on
    [FileName, PathName] = uigetfile({'*.asc','Arc-Info ASCII Raster (*.asc)'; ...
       '*.txt','Text ASCII Raster (*.txt)'; ...
       '*.*',  'All Files (*.*)'}, ...
       'Pick Raster(s)', ...
       pwd, ...
       'MultiSelect', 'on');
    % FILEPATH
    if iscell(FileName)
        for NumFile = 1:length(FileName)
            FilePath{NumFile} = strcat(PathName, FileName(NumFile));
        end
    else
        FilePath{1} = strcat(PathName, FileName);
        FileName = {FileName};
    end
end

% set undefined varargins
if ~exist('nodata', 'var'),     nodata  = NaN;      end     % default
if ~exist('ref', 'var'),        ref     = 'h';      end     % default
if ~exist('type', 'var'),       type    = 's';      end     % default
if length(FileName)==1,         type    =  '';      end     % double

%'''''''''''''''''
% clc
fprintf('\nReading...\n')
%'''''''''''''''''

% START LOOP FOR EACH FILE TO BE IMPORTED
for NumFile = 1:size(FilePath,2)

    fprintf('\t-->[%s]\n', strvcat(FilePath{NumFile}))
    fprintf('\tOutput NODATA-Value:\t %f\n', nodata)
    fprintf('\tSpatial Reference Mat:\t %s\n', ref)
    if ~isempty(type),fprintf('\tGrid class type:\t %s\n', type),end
    
    start = 0;
    eval(['fid = fopen(''', strvcat(FilePath{NumFile}), ''',''r'');']);
    if fid==-1, error('Invalid file name!!'), end

    % ********************  Import HEADER  ****************************
    while start == 0;   % to be sure to read HEADER of file

        Read_header = fscanf(fid,'%s',1);
        % x- and y- llcenter are for instance used in ISATIS
        switch Read_header
            case {'ncols', 'NCOLS'}
                ncols = fscanf(fid,'%f',1);
            case {'nrows', 'NROWS'}
                nrows = fscanf(fid,'%f',1);
            case {'xllcorner', 'XLLCORNER', 'xllcenter'}
                xllcorner = fscanf(fid,'%f',1);
            case {'yllcorner', 'YLLCORNER', 'yllcenter'}
                yllcorner = fscanf(fid,'%f',1);
            case {'cellsize', 'CELLSIZE'}
                cellsize = fscanf(fid,'%f',1);
            case {'xdim'}
                xdim = fscanf(fid,'%f',1);
            case {'ydim'}
                ydim = fscanf(fid,'%f',1);                
            case {'NODATA_value', 'nodata_value', 'NODATA_VALUE'}
                NODATA_value = fscanf(fid,'%f',1);
                start = 1;
            otherwise
                error('Unrecognized %s in header',Read_header)
        end     %switch

    end         %while
    % *****************************************************************

    % ----------- pixel geometry -----------
    if exist('cellsize','var') % ==> the pixel is squared
        xdim = cellsize;
        ydim = cellsize;
    end
    % --------------------------------------  
    
    % -----------   Import RASTER   -----------
    import = fscanf(fid,'%f',[ncols nrows]);
    fclose(fid);
    % -----------------------------------------

    
    % +++++++++++++  Write No-Data Value  +++++++++++++++++++
    import(find(import == NODATA_value)) = nodata;
    % +++++++++++++++++++++++++++++++++++++++++++++++++++++++

    %mulsel = @(O,I,inv) [O '.' strrep(strvcat(FileName{NumFile}{:}(1:length(FileName{NumFile}{:})-4)),'-','_') ' = ' I inv ';'];
    mulsel = @(O,I,inv) [O '.' strrep(strvcat(FileName{NumFile}(1:length(FileName{NumFile})-4)),'-','_') ' = ' I inv ';'];
    % @@@@@@@@@@@@   Write RASTER   @@@@@@@@@@@@@@@
    if type == 's'                  % structure array of 2-D layers
        %eval(['Z.' strrep(strvcat(FileName{NumFile}{:}(1:length(FileName{NumFile}{:})-4)),'-','_') ' = import'';']);
        eval(mulsel('Z','import',''''))
    elseif type == 'd'              % double 3-D array stack
        Z(:,:,NumFile) = import';
    else                            % single 2-D layer
        Z = import';
    end
    % @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

    
    % -----------------  Write HEADER  -------------------------
    if nargout == 2
        % Reference Matrix
        if ref == 'r'
            % ulc: Upper Left Centre point, that is the point in the centre of
            % the first upper-left pixel. This point is requested by
            % 'makerefmat.m' function.
            ulc_X = xllcorner + xdim /2;
            ulc_Y = yllcorner - ydim /2 + nrows*ydim;
            % This is the R spatial-referencing matrix used in Mapping Toolbox
            temp_R = makerefmat(ulc_X,ulc_Y,xdim,-ydim);
        % Header Matrix
        elseif ref == 'h'
            if ~exist('cellsize','var')
                warning('Unable to make canonical header')
                temp_R = [ncols; nrows; xllcorner; yllcorner; xdim; ydim; NODATA_value];
            else
                % This is the common header used in Arc/Info ASCII-Raster
                temp_R = [ncols; nrows; xllcorner; yllcorner; cellsize; NODATA_value];
            end
        end     %if
        % MultiSelection
        if type == 's'
            %eval(['R.' strrep(strvcat(FileName{NumFile}{:}(1:length(FileName{NumFile}{:})-4)),'-','_') ' = temp_R;']);
            eval(mulsel('R','temp_R',''))
        else
            R = temp_R;
        end
    end         %nargout
    % -----------------------------------------------------
    
    
end             %NumFile


fprintf('...Completed!\n')

return

function [filename] = path2file(filepath)

fnd_b_slash = strfind(char(filepath), filesep);

if fnd_b_slash 
    STR = char(filepath);
    filename = cellstr(STR(fnd_b_slash(end)+1:end));
else
    filename = filepath;
end
return

Contact us