Calculate pixel values at increasing radius

ndvi_classes.m function could be used to calculate the number and percentage of pixels at increasing radius from a lake following certain criterion’s. The function was originally developed to count the number of pixels falling in specific NDVI classes around lakes but may be useful for different proxies and situations. The function use georeferenced tiff images as input. For details please refer to Aleman et al. (2013) or to the included help. The image used in the example can be downloaded here.

The function:

function [classe,classe_per]=ndvi_classes(ImageName,critere1,critere2)

% INPUT:  ImageName is the name of the geotiff image for imput
%         critere1 (and 2) criterion for classes evaluation
% 
% OUTPUT: classe: vector of NDVI classes at increasing radius from
%                 the lake center (number of pixel falling in class)
%         classe_per: Same as classe but in percentage of pixels
% EXAMPLE: Using the given Image 'classif_GBL_94.tif'
%{
        [classe classe_per]=ndvi_classes('classif_GBL_94.tif','>2','<4');
        %Plot the classe as a function of distance in km 
        %(1pixel=0.0301km, input 200 pixels when asked)
        plot((1:200)*0.0301,classe_per,'k-')
        ylabel('% pixels >2 and <4')
        xlabel('Distance from lake center (km)')
%}         
% Blarquez Olivier 2012
% blarquez@gmail.com

[Imge,dat] = geotiffread(ImageName);

imagesc(double(Imge))

select= input('Select point interactively? Y/N: ', 's');
if select=='Y'
title('Please double clic on lake center:')
p = impoint(gca,[]);
p = wait(p);
p=round(p);
elseif select=='N' 
valstring = input('Enter latitude and longitude: ', 's');
valparts = regexp(valstring, '[ ,]', 'split');
coords = str2double(valparts);
[~,p(2)]=min(abs(linspace(dat.Latlim(1),dat.Latlim(2),dat.RasterSize(1))-...
    coords(1)));
[~,p(1)]=min(abs(linspace(dat.Lonlim(1),dat.Lonlim(2),dat.RasterSize(2))-...
    coords(2)));
end

radius = input('Radius in pixels to be evaluated: ');
Imge(:,((p(1)+(radius)):end))=[];
Imge(:, 1:(p(1)-(radius)))=[];
Imge(((p(2)+(radius)):end),:)=[];
Imge(1:(p(2)-(radius)),:)=[];

n=length(Imge);

classe=zeros(radius,1);
classe_per=zeros(radius,1);
for r=1:radius
    warning('off','all')
imshow(Imge)
%imagesc(double(Imge))

hEllipse = imellipse(gca,[(n/2-(r)) n/2-(r) r*2 r*2]);
maskImage = hEllipse.createMask();
imshow(maskImage);

maskedImge = Imge .* cast(maskImage,class(Imge));
%imshow(maskedImge);
maskedImge(maskedImge==0)=NaN;
pixList=reshape(maskedImge,[],1);
pixList(isnan(pixList),:)=[];
 
classe(r,1)=length(pixList(eval(['pixList',critere1]) & ...
     eval(['pixList',critere2]),:));
classe_per(r,1)=classe(r,1)./length(pixList(:,1));
end
close all
    warning('on','all')
end

Leave a Reply

Your email address will not be published. Required fields are marked *