function [i] = hypsometry(raster, classes, options, Sa, Sb);
PURPOSE: calculate the hypsometric integral of a raster w/ n classes and
plot the hypsometry curve.
Returns scalar [i] with hypsometric integral
-------------------------------------------------------------------
USAGE: [i] = hypsometry(raster, classes, options , Sa, Sb)
where: [raster] is the input (elevation) grid of size (cols,rows)
[classes] (optional) is the number of classes to use, default is
total number of elements
[options] (optional) vector of form [1 1] specifying whether to
normalize hypsometry (first element) and whether to plot
curve (second element)(1=yes,0=no), default is [1 1]
[Sa] (optional) is a character string made from one element
from any or all the following 3 columns, eg 'r.':
b blue . point - solid
g green o circle : dotted
r red x x-mark -. dashdot
c cyan + plus -- dashed
m magenta * star (none) no line
y yellow s square
k black d diamond
v triangle (down)
^ triangle (up)
< triangle (left)
> triangle (right)
p pentagram
h hexagram
[Sb] (optional) is a vector specifying linewidth and markersize in
the form [lw,ms]. To set Sb, all other options have to be set as well.
-------------------------------------------------------------------------
OUTPUTS:
[i] scalar with the hypsometric integral, based on cellsize 1 and
number of classes. For true area, multiply with cellsize.
-------------------------------------------------------------------
EXAMPLE: hypsometry(abs(randn(100,100)),20,[1 1],'ro-',[2 2])
NOTES: NaNs are ignored. If class numbers are provided, mean for equal
size classes are used for hypsometry calculation.
SEE ALSO: plot, numel
Felix Hebeler, Geography Dept., University Zurich, March 20070001 function [i] = hypsometry(raster, classes, options ,Sa, Sb) 0002 % function [i] = hypsometry(raster, classes, options, Sa, Sb); 0003 % PURPOSE: calculate the hypsometric integral of a raster w/ n classes and 0004 % plot the hypsometry curve. 0005 % Returns scalar [i] with hypsometric integral 0006 % ------------------------------------------------------------------- 0007 % USAGE: [i] = hypsometry(raster, classes, options , Sa, Sb) 0008 % where: [raster] is the input (elevation) grid of size (cols,rows) 0009 % [classes] (optional) is the number of classes to use, default is 0010 % total number of elements 0011 % [options] (optional) vector of form [1 1] specifying whether to 0012 % normalize hypsometry (first element) and whether to plot 0013 % curve (second element)(1=yes,0=no), default is [1 1] 0014 % [Sa] (optional) is a character string made from one element 0015 % from any or all the following 3 columns, eg 'r.': 0016 % 0017 % b blue . point - solid 0018 % g green o circle : dotted 0019 % r red x x-mark -. dashdot 0020 % c cyan + plus -- dashed 0021 % m magenta * star (none) no line 0022 % y yellow s square 0023 % k black d diamond 0024 % v triangle (down) 0025 % ^ triangle (up) 0026 % < triangle (left) 0027 % > triangle (right) 0028 % p pentagram 0029 % h hexagram 0030 % 0031 % [Sb] (optional) is a vector specifying linewidth and markersize in 0032 % the form [lw,ms]. To set Sb, all other options have to be set as well. 0033 % ------------------------------------------------------------------------- 0034 % OUTPUTS: 0035 % [i] scalar with the hypsometric integral, based on cellsize 1 and 0036 % number of classes. For true area, multiply with cellsize. 0037 % ------------------------------------------------------------------- 0038 % EXAMPLE: hypsometry(abs(randn(100,100)),20,[1 1],'ro-',[2 2]) 0039 % 0040 % NOTES: NaNs are ignored. If class numbers are provided, mean for equal 0041 % size classes are used for hypsometry calculation. 0042 % 0043 % SEE ALSO: plot, numel 0044 % 0045 % Felix Hebeler, Geography Dept., University Zurich, March 2007 0046 if nargin < 1 0047 error('This function needs some input to work!'); 0048 elseif nargin > 5 0049 error('This is too much information! Please give no more than 5 arguments (raster,classes,plot,Sa,Sb)!'); 0050 end 0051 0052 0053 % sort and eliminate NaNs 0054 raster=sort(reshape(raster,1,numel(raster)),'descend'); 0055 raster(isnan(raster))=[]; 0056 elements=numel(raster); 0057 mr=max(raster); 0058 i=1; 0059 % check if data needs to be resampled to classes 0060 if exist('classes','var') 0061 if ischar(classes) 0062 Sa=classes; 0063 classes=elements; 0064 elseif size(classes,2)<2 % check if classes is not the options flag 0065 isize=floor(elements/classes); 0066 t=nan(1,classes); 0067 for c=1:classes 0068 t(c)=nanmean(raster( ((c-1)*isize)+1:(c*isize) )); 0069 end 0070 raster=t; 0071 clear t; 0072 i=isize; % set i to number of classes for correct integral 0073 else 0074 options = classes; 0075 classes=elements; 0076 end 0077 else 0078 classes=elements; 0079 end 0080 % check if hypsometry curve should be plotted 0081 if ~exist('options','var') 0082 options=[1 1]; 0083 else 0084 if ischar(options) %make sure its not Sa we have here 0085 Sa=options; 0086 options=[1 1]; 0087 elseif max(options)~=1 && min(options)~=0 0088 options=[1 1]; 0089 end 0090 end 0091 0092 i=i*sum(raster); % integral is sum of elevation, times classes, if used 0093 if options(1)==1 0094 i=i/(mr*elements); %normalize 0095 raster=(raster-min(raster))./(max(raster)-min(raster)); 0096 end 0097 if options(2)==1 0098 if ~exist('Sa','var') 0099 Sa='r-'; %set default color here if to override system preferences 0100 end 0101 if ~exist('Sb','var') 0102 Sb=[1,3]; %defaults: LineWidth =1, MarkerSize = 3 0103 end 0104 % converting format strings so plot will accept them 0105 pa=['''LineWidth''',',',num2str(Sb(1)),',','''MarkerSize''',',',num2str(Sb(2))]; 0106 Sa=['''',Sa,'''']; 0107 0108 eval(['plot(linspace(0,100,numel(raster)),raster,',Sa,',',pa,');']); 0109 title(['Hypsometry curve. Integral i=',num2str(i),', number of classes = ',num2str(classes)]); 0110 xlabel('Cumulative area [%]') 0111 if options(1)==1 0112 ylim([0 1]); 0113 ylabel('Normalised y-Value') 0114 else 0115 ylim([min(raster) (max(raster)+max(raster)/20)]) 0116 ylabel('y-Value') 0117 end 0118 end