function [Z R] = ImportAsciiRaster(nodata, ref, type);
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% %
% Produced by Giuliano Langella %
% e-mail:gyuliano@libero.it %
% March 2008 %
% %
% Last Updated: 08 October, 2008 %
% %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%
%
% ----- SYNTAX -----
% Z = ImportAsciiRaster ==> nodata:NaN; R: []; type:'s';
% [Z R] = ImportAsciiRaster ==> nodata:NaN; ref:'r'; type:'s';
% [Z R] = ImportAsciiRaster(nodata) ==> nodata:nodata; ref:'r'; type:'s';
% [Z R] = ImportAsciiRaster(nodata, ref); ==> nodata:nodata; ref:ref; type:'s';
% [Z R] = ImportAsciiRaster(nodata, ref, type); ==> nodata:nodata; ref:ref; type:type;
%
%
% ----- OPTIONS -----
% +'nodata'(optional) = 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): Geospatial Reference Matrix
% --> 'r' (default) for spatial-referencing matrix R (Mapping Toolbox);
% --> 'h' for header matrix (as the .hdr file in Arc-Info Binary Raster).
%
% +'type'(optional): class type of Z output in multiselection (R is a single reference matrix)
% --> '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.
%
%
% ----- 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; % 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
%
%
% -----------------------------------
% LAST USED:
% [Z1 hdr] = ImportAsciiRaster(-1,'h');
% -----------------------------------
if nargin == 0, nodata = NaN; end
% 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:size(FileName,2)
FilePath{NumFile} = strcat(PathName, FileName(NumFile));
end
else
FilePath{1} = strcat(PathName, FileName);
end
% $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
%'''''''''''''''''
clc
disp('ELABORATING...')
%'''''''''''''''''
% START LOOP FOR EACH FILE TO BE IMPORTED
for NumFile = 1:size(FilePath,2)
disp(['...' strvcat(FilePath{NumFile}) '...'])
disp(['No-Data Value: ' num2str(nodata)])
start = 0;
eval(['fid = fopen(''', strvcat(FilePath{NumFile}), ''',''r'');']);
% ******************** Import HEADER ****************************
while start == 0; % to be sure to read HEADER of file
Read_header = fscanf(fid,'%s',1);
switch Read_header
case {'ncols', 'NCOLS'}
ncols = fscanf(fid,'%f',1);
case {'nrows', 'NROWS'}
nrows = fscanf(fid,'%f',1);
case {'xllcorner', 'XLLCORNER'}
xllcorner = fscanf(fid,'%f',1);
case {'yllcorner', 'YLLCORNER'}
yllcorner = fscanf(fid,'%f',1);
case {'cellsize', 'CELLSIZE'}
cellsize = fscanf(fid,'%f',1);
case {'NODATA_value', 'nodata_value', 'NODATA_VALUE'}
NODATA_value = fscanf(fid,'%f',1);
start = 1;
end %switch
end %while
% *****************************************************************
% ----------- Import RASTER -----------
import = fscanf(fid,'%f',[ncols nrows]);
fclose(fid);
% -----------------------------------------
% +++++++++++++ Write No-Data Value +++++++++++++++++++
import(find(import == NODATA_value)) = nodata;
% +++++++++++++++++++++++++++++++++++++++++++++++++++++++
% @@@@@@@@@@@@ Write RASTER @@@@@@@@@@@@@@@
if iscell(FileName) && (nargin < 3 || type == 's') %structure array
eval(['Z.' strvcat(FileName{1,NumFile}(1:length(FileName{1,NumFile})-4)) ' = import'';']);
elseif iscell(FileName) && type == 'd' %double array stack
Z(:,:,NumFile) = import';
else %single layer
Z = import';
end
% @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
% ----------------- Write HEADER -------------------------
if nargout == 2
% Reference Matrix
if nargin < 2 | (nargin > 1 & 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 + cellsize /2;
ulc_Y = yllcorner - cellsize /2 + nrows*cellsize;
% This is the R spatial-referencing matrix used in Mapping Toolbox
temp_R = makerefmat(ulc_X,ulc_Y,cellsize,-cellsize);
% Header Matrix
elseif nargin > 1 & ref == 'h'
% This is the common header used in Arc/Info ASCII-Raster
temp_R = [ncols; nrows; xllcorner; yllcorner; cellsize; NODATA_value];
end %if
% MultiSelection
if iscell(FileName) && (nargin < 3 || type == 's')
eval(['R.' strvcat(FileName{1,NumFile}(1:length(FileName{1,NumFile})-4)) ' = temp_R;']);
else
R = temp_R;
end
end %nargout
% -----------------------------------------------------
end %NumFile
disp('...COMPLETED!')
