%% Intensity Analysis of Spatial Point Patterns
% *Phaedon Kyriakidis* 
%
% Department of Geography

%% Import Lansing tree data
load lansingTrees.mat
who

%% Contents of file Lansingtrees
lansInfo

%% 
disp( lansTrees(1:10,:) )

%% Some variable assignment to make our life easier
x = lansTrees(:,1);
y = lansTrees(:,2);
xy = [x y];
c = lansTrees(:,3); % 3rd column = integer code (1,2,...,6) with tree type

%%
plot(x,y,'*');
axis equal square

%% Group-specific scatter plots
h = gscatter(x,y,c);
set(h,'MarkerSize',5);
axis equal square
legend(lansLabels);

%% Declare variables for hickory and maple trees
xH = x(c==2);  % x-coords of hickory trees
xM = x(c==3);  % x-coords of maple trees
yH = y(c==2);  % y-coords of hickory trees
yM = y(c==3);  % y-coords of maple trees

%% 
plot(xH,yH,'r*',xM,yM,'b+');
axis equal square
legend('Hickory','Maple');
title('Locations of Hickory (*) and Maple (+) trees');


%% Scatter plot and marginal histograms for Hickory
nBins = [50 40];
scatterhist(xH,yH,nBins);

%% Scatter plot and marginal histograms for Maple
nBins = [50 40];
scatterhist(xM,yM,nBins);

%% QUADRAT COUNTS

%% 2D frequency histogram of Hickory counts
nBins = [30 20];
figure
hist3([xH yH],nBins)
xlabel('x')
ylabel('y')
zlabel('freq')
title('Hickory intensity');

%% 2D frequency histogram of Maple counts
nBins = [20 10];
figure
hist3([xM yM],nBins)
xlabel('x')
ylabel('y')
zlabel('freq')
title('Maple intensity')

%% Specifying a grid of quadrats
figure
QuadGridSpecs = quadratgridspecs(xy,'dataplot');
% QuadGridSpecs is a (nDim x 3) array with quadrat grid specifications:
% [nRecX xminRec xsizeRec; nRecY yminRec ysizeRec]
axis square
title('Quadrats')
Q = quadratgridspecs(xy,'dataplot',[0 1; 0 1],[10 10]); % more options

%% Points in quadrats (map) for Hickory
QuadGridSpecs = [10 0 0.1; 10 0 0.1] % 
QuadCountsH = point2quadrat([xH yH],'freq',QuadGridSpecs,'dataplot');
title('Hickory quadrat counts');

%% Frequency histogram of quadrats (Hickory)
TabH = tabulate(QuadCountsH) 
% 1st column = point count = # of points in quadrat
% 2nd column = # of quadrats with corresponding point count
% 3rd column = same as column #2, but in terms of percentages
stem(TabH(:,1),TabH(:,2),'filled')
title('Hickory: Quadrat frequency per point count');
xlabel('point count')
ylabel('number of quadrats')

%% Points in quadrats (map) for Maple
QuadGridSpecs = [10 0 0.1; 10 0 0.1]
QuadCountsM = point2quadrat([xM yM],'freq',QuadGridSpecs,'dataplot');
title('Maple quadrat counts');

%% Frequency histogram of quadrats (Maple)
figure
TabM = tabulate(QuadCountsM) 
% 1st column = point count = # of points in quadrat
% 2nd column = # of quadrats with corresponding point count
% 3rd column = same as column #2, but in terms of percentages
stem(TabM(:,1),TabM(:,2),'filled')
title('Maple: Quadrat frequency per point count');
xlabel('point count')
ylabel('number of quadrats')

%% Frequency histogram of quadrats (Hickory) -- custom function
[CountH,CountStatsH] = point2quadratcount([xH yH],QuadGridSpecs);
% Outputs:
% CountH(:,1) = count number (0,1,2,...)
% CountH(:,2) = # of quadrats with corresponding count
% CountStatsH = [mQ sQ indexDispers chi2val] 
% [mQ = mean quadrat count, sQ = std deviation of quadrat counts,
% indexDispers = varQ/mQ, chi2val = (N-1)*varQ/mQ]
stem(CountH(:,1),CountH(:,2),'filled')
title('Hickory: Quadrat frequency histogram')

%% Frequency histogram of quadrats (Maple) -- custom function
[CountM,CountStatsM] = point2quadratcount([xM yM],QuadGridSpecs);
% Outputs:
% CountM(:,1) = count number (0,1,2,...)
% CountM(:,2) = # of quadrats with corresponding count
% CountStatsM = [mQ sQ indexDispers chi2val] 
% [mQ = mean quadrat count, sQ = std deviation of quadrat counts,
% indexDispers = varQ/mQ, chi2val = (N-1)*varQ/mQ]
stem(CountM(:,1),CountM(:,2),'filled')
title('Maple: Quadrat frequency histogram')

%% INTENSITY ESTIMATION

%% Specification of grid of points for estimating intensity
figure
minmax = []; % min/max coordinates of bounding box 
% An empty array implies these will be taken from the data
GridSpecs = makegridspecs(xy,minmax,[10 10],'dataplot')

%% A more dense grid
GridSpecs = makegridspecs(xy,minmax,[100 100],'dataplot')
%GridSpecs = makegridspecs(xy,[0.01 99.99; 0.01 99.99],[100 100],'dataplot',[0.01 0.01])

%% Available kernel types and corresponding integer codes
kerneltypes

%% Intensity estimation for Hickory
w = []; % No weights = equally important tree locations
kernekType = 3; % Gaussian kernel
bandWidths = []; % Let the program choose default bandwidths
% 'raw' -> output density not standardized
% 'one' -> output density integrates to 1
% 'ndata' -> output density integrates to number of events
K = kerneldensity([xH yH],[1 2],w,GridSpecs,3,[0.05 0.05],'raw','dataplot');

%% Intensity estimation for Maple
w = []; % equally weighted/important tree locations
figure
K = kerneldensity([xM yM],[1 2],w,GridSpecs,3,bandWidths,'raw','dataplot');