classdef GRIDobj %GRIDobj Create instance of a GRIDobj % % Syntax % % DEM = GRIDobj(X,Y,dem) % DEM = GRIDobj('ESRIasciiGrid.txt') % DEM = GRIDobj('GeoTiff.tif') % DEM = GRIDobj(); % DEM = GRIDobj([]); % DEM = GRIDobj(FLOWobj or GRIDobj or STREAMobj,class) % % % Description % % GRIDobj creates an instance of the grid class, which contains a % numerical or logical matrix and information on georeferencing. When a % GRIDobj is created from a file, the number format of the data in % GRIDobj is either single or double. Unsigned and signed integers are % converted to single. For unsigned integers, missing values are % assumed to be denoted as intmax(class(input)). For signed integers, % missing values are assumed to be intmin(class(input)). Please check, % that missing values in your data have been identified correctly % before further analysis. % % Note that while throughout this help text GRIDobj is associated with % gridded digital elevation models, instances of GRIDobj can contain % other gridded, single band, datasets such as flow accumulation grids, % gradient grids etc. % % DEM = GRIDobj(X,Y,dem) creates a DEM object from the coordnate % matrices or vectors X and Y and the matrix dem. The elements of dem % refer to the elevation of each pixel. % % DEM = GRIDobj('ESRIasciiGrid.txt') creates a DEM object from an ESRI % Ascii grid exported from other GI systems. % % DEM = GRIDobj('GeoTiff.tif') creates a DEM object from a Geotiff. % % DEM = GRIDobj() opens a dialog box to read either an ESRI Ascii Grid % or a Geotiff. % % DEM = GRIDobj([]) creates an empty instance of GRIDobj % % DEM = GRIDobj(FLOWobj or GRIDobj or STREAMobj,class) creates an % instance of GRIDobj with all common properties (e.g., spatial % referencing) inherited from another instance of a FLOWobj, GRIDobj % or STREAMobj class. DEM.Z is set to all zeros where class can be % integer classes or double or single. By default, class is double. % % Example % % % Load DEM % DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif'); % % Display DEM % imageschs(DEM) % % See also: FLOWobj, STREAMobj, GRIDobj/info % % Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de) % Date: 17. August, 2017 properties %Public properties %Z matrix with elevation values % The Z property contains the elevation values in a 2D matrix. % % See also GRIDobj Z %CELLSIZE cellsize of the grid (scalar) % The cellsize property specifies the spacing of the grid in x and y % directions. Note that TopoToolbox requires grids to have square cells, % e.g., dx and dy are the same. % % See also GRIDobj cellsize %REFMAT 3-by-2 affine transformation matrix % The refmat property specifies a 3-by-2 affine transformation matrix as % used by the mapping toolbox. % % See also GRIDobj, makerefmat refmat %SIZE size of the grid (two element vector) % The cellsize property is a two element vector that contains the number % of rows and columns of the Z matrix. % % See also GRIDobj, size size %NAME optional name (string) % The name property allows to specify a name of the grid. By default and % if the constructor is called with a filename, the name property is set % to the name of the file. % % See also GRIDobj name %ZUNIT unit of grid values (string) % The zunit is optional and is used to store the physical unit (e.g. m) % of an instance of GRIDobj. This property is currently not fully % supported and TopoToolbox functions usually assume that the unit is in % meters and equals the xyunit property. % % See also GRIDobj zunit %XYUNIT unit of the coordinates (string) % The xyunit is optional and is used to store the physical unit (e.g. m) % of the coordinates. This property is currently not fully % supported and TopoToolbox functions usually assume that the unit is in % meters and equals the zunit property. % % See also GRIDobj xyunit %GEOREF additional information on spatial referencing (structure array) % The georef property stores an instance of map.rasterref.MapCellsReference, % a structure array GeoKeyDirectoryTag, and a mapping structure % (mstruct). % % See also GRIDobj, geotiffinfo georef end methods function DEM = GRIDobj(varargin) % GRIDobj constructor if nargin == 3 %% GRIDobj is created from three matrices % GRIDobj(X,Y,dem) X = varargin{1}; Y = varargin{2}; if min(size(X)) > 1 X = X(1,:); end if min(size(Y)) > 1 Y = Y(:,1); end DEM.Z = varargin{3}; DEM.size = size(DEM.Z); if numel(X) ~= DEM.size(2) || numel(Y) ~= DEM.size(1) error('TopoToolbox:GRIDobj',... ['Coordinate matrices/vectors don''t fit the size of the \n'... 'the grid']); end if (Y(2)-Y(1)) > 0 % the y coordinate vector must be monotonically % decreasing, so that the left upper edge of the DEM is % north-west (on the northern hemisphere). DEM.Z = flipud(DEM.Z); Y = Y(end:-1:1); end dy = Y(2)-Y(1); dx = X(2)-X(1); if abs(abs(dx)-abs(dy))>1e-9 error('TopoToolbox:GRIDobj',... 'The resolution in x- and y-direction must be the same'); end DEM.refmat = double([0 dy;... dx 0;... X(1)-dx Y(1)-dy]); DEM.cellsize = dx; DEM.georef = []; DEM.name = []; elseif nargin <= 2 if nargin == 0 %% No input arguments. File dialog box will open and ask % for a txt or tiff file as input FilterSpec = {'*.txt;*.asc;*.tif;*.tiff','supported file types (*.txt,*.asc,*.tif,*.tiff)';... '*.txt', 'ESRI ASCII grid (*.txt)';... '*.asc', 'ESRI ASCII grid (*.asc)';... '*.tif', 'GeoTiff (*.tif)';... '*.tiff', 'GeoTiff (*.tiff)';... '*.*', 'all files (*.*)'}; DialogTitle = 'Select ESRI ASCII grid or GeoTiff'; [FileName,PathName] = uigetfile(FilterSpec,DialogTitle); if FileName == 0 error('TopoToolbox:incorrectinput',... 'no file was selected') end filename = fullfile(PathName, FileName); elseif nargin > 0 % One input argument if isempty(varargin{1}) % if empty array than return empty GRIDobj return elseif isa(varargin{1},'GRIDobj') || ... isa(varargin{1},'FLOWobj') || ... isa(varargin{1},'STREAMobj') % empty GRIDobj DEM = GRIDobj([]); % find common properties of F and G and from F to G pg = properties(DEM); pf = properties(varargin{1}); p = intersect(pg,pf); for r = 1:numel(p) DEM.(p{r}) = varargin{1}.(p{r}); end if nargin == 1 cl = 'double'; else cl = varargin{2}; end if strcmp(cl,'logical') DEM.Z = false(DEM.size); else DEM.Z = zeros(DEM.size,cl); end DEM.name = ''; return end % GRIDobj is created from a file filename = varargin{1}; end % check if file exists if exist(filename,'file')~=2 error('File doesn''t exist') end % separate filename into path, name and extension [pathstr,DEM.name,ext] = fileparts(filename); if any(strcmpi(ext,{'.tif', '.tiff'})) % it is a GeoTiff try % try to read using geotiffread (requires mapping % toolbox) [DEM.Z, DEM.refmat, ~] = geotiffread(filename); gtiffinfo = geotiffinfo(filename); DEM.georef.SpatialRef = gtiffinfo.SpatialRef; DEM.georef.GeoKeyDirectoryTag = gtiffinfo.GeoTIFFTags.GeoKeyDirectoryTag; georef_enabled = true; catch ME % mapping toolbox is not available. Will try to % read the tif file together with the tfw file georef_enabled = false; % the tfw file has the same filename but a tfw % extension tfwfile = fullfile(pathstr,[DEM.name '.tfw']); % check whether file exists. If it exists then read % it using worldfileread or own function tfwfile_exists = exist(tfwfile,'file'); if tfwfile_exists try % prefer builtin worldfileread, if % available DEM.refmat = worldfileread(tfwfile); catch ME W = dlmread(tfwfile); DEM.refmat(2,1) = W(1,1); DEM.refmat(1,2) = W(4,1); DEM.refmat(3,1) = W(5,1)-W(1); DEM.refmat(3,2) = W(6,1)-W(4); end DEM.Z = imread(filename); else if ~tfwfile_exists error('TopoToolbox:GRIDobj:read',... 'GRIDobj cannot read the TIF-file because it does not have a tfw-file.'); else throw(ME) end end end % Unless any error occurred, we now attempt to generate % an mapping projection structure. This will not work % if the DEM is in a geographic coordinate system or if % the projection is not supported by mstruct. if georef_enabled try DEM.georef.mstruct = geotiff2mstruct(gtiffinfo); catch DEM.georef.mstruct = []; warning('TopoToolbox:GRIDobj:projection',... ['GRIDobj cannot derive a map projection structure. This is either\n' ... 'because the grid is in a geographic coordinate system or because\n' ... 'geotiff2mstruct cannot identify the projected coordinate system used.\n' ... 'TopoToolbox assumes that horizontal and vertical units of DEMs are \n'... 'the same. It is recommended to use a projected coordinate system,\n' ... 'preferably UTM WGS84. Use the function GRIDobj/reproject2utm\n' ... 'to reproject your grid.']) end end % Finally, check whether no_data tag is available. This tag is % not accessible using geotiffinfo (nice hack by Simon % Riedl) tiffinfo = imfinfo(filename); if isfield(tiffinfo,'GDAL_NODATA') nodata_val = str2double(tiffinfo.GDAL_NODATA); end else [DEM.Z,R] = rasterread(filename); DEM.refmat = R; DEM.georef = []; end DEM.size = size(DEM.Z); DEM.cellsize = abs(DEM.refmat(2)); % remove nans demclass = class(DEM.Z); nodata_val_exists = exist('nodata_val','var'); switch demclass case {'uint8','uint16','uint32'} % unsigned integer DEM.Z = single(DEM.Z); if nodata_val_exists nodata_val = single(nodata_val); DEM.Z(DEM.Z == nodata_val) = nan; else DEM.Z(DEM.Z==intmax(demclass)) = nan; end case {'int8','int16','int32'} % signed integer DEM.Z = single(DEM.Z); if nodata_val_exists nodata_val = single(nodata_val); DEM.Z(DEM.Z == nodata_val) = nan; else DEM.Z(DEM.Z==intmin(demclass)) = nan; end case {'double','single'} if nodata_val_exists DEM.Z(DEM.Z == cast(nodata_val,class(DEM.Z))) = nan; end case 'logical' otherwise error('TopoToolbox:GRIDobj','unrecognized class') end end end end end % Subfunction for ASCII GRID import function [Z,refmat] = rasterread(file) fid=fopen(file,'r'); % loop through header header = struct('ncols',[],... 'nrows',[],... 'xllcorner',[],... 'yllcorner',[],... 'cellsize',[],... 'nodata',[]); names = fieldnames(header); nrnames = numel(names); try fseek(fid,0,'bof'); for r = 1:nrnames ; headertext = fgetl(fid); [headertext, headernum] = strtok(headertext,' '); I = cellfun(@(x,y) strcmpi(x(1:4),y(1:4)),names,repmat({headertext},nrnames,1)); header.(names{I}) = str2double(headernum); end catch ME1 error('header can not be read') end % read raster data Z = fscanf(fid,'%lg',[header.ncols header.nrows]); fclose(fid); Z(Z==header.nodata) = NaN; Z = Z'; % create X and Y using meshgrid refmat = [0 -header.cellsize;... header.cellsize 0;... header.xllcorner+(0.5*header.cellsize) - header.cellsize ... (header.yllcorner+(0.5*header.cellsize))+((header.nrows)*header.cellsize)]; end