How to put conditions to find peaks in NDVI time series?

Question

Devendra 2024-4-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2103601-how-to-put-conditions-to-find-peaks-in-ndvi-time-series

评论： gauri 2024-4-26

ndvi_time_series_input 2.csv

I am using following matlab code findpeaks(data,days'MinPeakProminence',60) xlabel('days’) ylabel('NDVI') title('Find Prominent Peaks') here the data is NDVI time series attached here.

9 个评论
显示 7更早的评论隐藏 7更早的评论

Manikanta Aditya 2024-4-6

编辑：Manikanta Aditya 2024-4-6

在 MATLAB Online 中打开

To find peaks in an NDVI time series data using MATLAB, you can use the findpeaks function. The findpeaks function allows you to set various parameters to control the peak detection process.

% Load the NDVI time series data
ndvi_data = [
    % add your NDVI data
];
% Define the date vector (assuming the dates are in the order provided)
dates = datetime(2023, 4, 9) + calmonths(0:size(ndvi_data, 2)-1);
% Loop over each location
for i = 1:size(ndvi_data, 1)
    % Extract the NDVI time series for the current location
    location_data = ndvi_data(i, :);
    
    % Find prominent peaks with a minimum prominence of 0.1
    [pks, locs] = findpeaks(location_data, dates, 'MinPeakProminence', 0.1);
    
    % Plot the NDVI time series and the detected peaks
    figure;
    plot(dates, location_data);
    hold on;
    plot(dates(locs), pks, 'r*');
    
    % Add labels and title
    xlabel('Date');
    ylabel('NDVI');
    title(sprintf('Find Prominent Peaks for Location %d', i));
    
    % Adjust the x-axis ticks for better visibility
    xtickangle(45);
    datetick('x', 'mmm yyyy', 'keeplimits');
end

Manikanta Aditya 2024-4-6

在 MATLAB Online 中打开

Hi,

To modify the code to find peaks with NDVI values greater than 0.4 and a minimum distance of 60 days between peaks, and to return the number of peaks as a vector, you can make the following changes:

% Load the NDVI time series data
ndvi_data = [
    % add your NDVI data
];
% Define the date vector (assuming the dates are in the order provided)
dates = datetime(2023, 4, 9) + calmonths(0:size(ndvi_data, 2)-1);
% Preallocate a vector to store the number of peaks for each location
num_peaks = zeros(size(ndvi_data, 1), 1);
% Loop over each location
for i = 1:size(ndvi_data, 1)
    % Extract the NDVI time series for the current location
    location_data = ndvi_data(i, :);
    
    % Find peaks with a minimum prominence of 0.1, NDVI value greater than 0.4,
    % and a minimum distance of 60 days between peaks
    [pks, locs, ~, prominences] = findpeaks(location_data, dates, ...
        'MinPeakProminence', 0.1, ...
        'MinPeakHeight', 0.4, ...
        'MinPeakDistance', 60);
    
    % Count the number of peaks for the current location
    num_peaks(i) = numel(pks);
    
    % Plot the NDVI time series and the detected peaks (optional)
    figure;
    plot(dates, location_data);
    hold on;
    plot(dates(locs), pks, 'r*');
    xlabel('Date');
    ylabel('NDVI');
    title(sprintf('Find Prominent Peaks for Location %d', i));
    xtickangle(45);
    datetick('x', 'mmm yyyy', 'keeplimits');
end

The findpeaks function is used with specific parameters to detect peaks in NDVI values greater than 0.4 and at least 60 days apart. The number of these peaks is counted and stored.

Hope this helps, let me know, I will post as an answer, then you can accept the answer.

Devendra 2024-4-6

编辑：Devendra 2024-4-7

Actually, right now I am traveling and don’t have my laptop. I will be back to Delhi my home and will run the code. But I am pretty sure it will work.

Further May I request you to please help me in finding some additional informations from the same NDVI file which is as follows;

POS value; MAX sum; AVG sum; SOS; POS; EOS; POS value; Average sum; Maximum sum; Base; Growth rate; Senescence rate;

I want to keep all these parameters in one array in addition to peaks information.

The description of above mentioned parameters are as follows; Start of season (SOS): the date of growth start and crop emergence • Peak of season (POS): the date of the maximum NDVI value in the time series • End of season (EOS): the date of the harvesting• POS value: maximum NDVI value in the time series • Average sum: sum of average NDVI values in the time series • Maximum sum: sum of maximum NDVI values in the time series • Base: the range between minimum and maximum NDVI. The SOS is given as the rst date when the NDVI crosses a threshold of 0.2 and continues in a monotonic growth. The POS is the absolute NDVI maximum within the season, whereas the EOS is the first date after the POS when the NDVI time series crosses a threshold of 0.2 and continues to decrease. The growth and senescence rates are the gradients between SOS and POS and EOS respectively.

I am also attaching the output file as a sample output file.

I am very sorry for causing inconvenience to you.

I will appreciate your kind help.

Devendra

Manikanta Aditya 2024-4-7

编辑：Manikanta Aditya 2024-4-7

在 MATLAB Online 中打开

Hi,

Check this, based on your explaination this is what I was able to come up with:

Replace with correct file, and if any small issues you should be able to fix them.

% Load the NDVI time series data
ndvi_data = readmatrix('ndvi_time_series_input 2.csv');
% Define the date vector (assuming the dates are in the order provided)
dates = datetime(2023, 4, 9) + calmonths(0:size(ndvi_data, 2)-1);
% Initialize an array to store the phenology parameters
phenology_params = nan(size(ndvi_data, 1), 10);  % Only numeric values
phenology_dates = NaT(size(ndvi_data, 1), 3);  % Only datetime values
% Loop over each location
for i = 1:size(ndvi_data, 1)
    % Extract the NDVI time series for the current location
    location_data = ndvi_data(i, :);
    
    % Find peaks with a minimum prominence of 0.1, NDVI value greater than 0.4,
    % and a minimum distance of 10 days between peaks
    [pks, locs, ~, prominences] = findpeaks(location_data, 'MinPeakProminence', 0.1, 'MinPeakHeight', 0.4, 'MinPeakDistance', 10);
    
    % Count the number of peaks for the current location
    num_peaks = numel(pks);
    
    % Determine SOS, POS, and EOS
    sos_idx = find(location_data >= 0.2, 1, 'first');
    pos_idx = locs(location_data(locs) == max(pks));
    eos_idx = find(location_data(pos_idx:end) < 0.2, 1, 'first') + pos_idx - 1;
    
    % Calculate other phenology parameters
    pos_value = location_data(pos_idx);
    max_sum = sum(max(location_data));
    avg_sum = sum(location_data);
    base = max(location_data) - min(location_data);
    duration = eos_idx - sos_idx + 1;
    first_half = pos_idx - sos_idx;
    second_half = eos_idx - pos_idx;
    growth_rate = (pos_value - location_data(sos_idx)) / first_half;
    senescence_rate = (location_data(eos_idx) - pos_value) / second_half;
    
    % Store the phenology parameters
    phenology_params(i, :) = [pos_value, max_sum, avg_sum, base, duration, first_half, second_half, growth_rate, senescence_rate, num_peaks];
    phenology_dates(i, :) = [dates(sos_idx), dates(pos_idx), dates(eos_idx)];
end

The code calculates various phenology parameters such as POS value, MAX sum, AVG sum, SOS, POS, EOS, Base, Duration, First half, Second half, Growth rate, Senescence rate, and Number of peaks for each location. These parameters are stored in the phenology_params array for future use or analysis.

Hope it helps @Devendra!

Manikanta Aditya 2024-4-7

Thank you! Feel free to reach out if needed any help!

Devendra 2024-4-25

Thank you very much for your kind help. I will try to understand it and will try to resolve it. Certainly I will approach you at the critical time.

Once again I am very much thankful to you for your generous and kind help.

Devendra

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Manikanta Aditya 2024-4-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2103601-how-to-put-conditions-to-find-peaks-in-ndvi-time-series#answer_1437586

在 MATLAB Online 中打开

Hi,

Check this, based on your explaination this is what I was able to come up with:

Replace with correct file, and if any small issues you should be able to fix them.

% Load the NDVI time series data
ndvi_data = readmatrix('ndvi_time_series_input 2.csv');
% Define the date vector (assuming the dates are in the order provided)
dates = datetime(2023, 4, 9) + calmonths(0:size(ndvi_data, 2)-1);
% Initialize an array to store the phenology parameters
phenology_params = nan(size(ndvi_data, 1), 10);  % Only numeric values
phenology_dates = NaT(size(ndvi_data, 1), 3);  % Only datetime values
% Loop over each location
for i = 1:size(ndvi_data, 1)
    % Extract the NDVI time series for the current location
    location_data = ndvi_data(i, :);
    
    % Find peaks with a minimum prominence of 0.1, NDVI value greater than 0.4,
    % and a minimum distance of 10 days between peaks
    [pks, locs, ~, prominences] = findpeaks(location_data, 'MinPeakProminence', 0.1, 'MinPeakHeight', 0.4, 'MinPeakDistance', 10);
    
    % Count the number of peaks for the current location
    num_peaks = numel(pks);
    
    % Determine SOS, POS, and EOS
    sos_idx = find(location_data >= 0.2, 1, 'first');
    pos_idx = locs(location_data(locs) == max(pks));
    eos_idx = find(location_data(pos_idx:end) < 0.2, 1, 'first') + pos_idx - 1;
    
    % Calculate other phenology parameters
    pos_value = location_data(pos_idx);
    max_sum = sum(max(location_data));
    avg_sum = sum(location_data);
    base = max(location_data) - min(location_data);
    duration = eos_idx - sos_idx + 1;
    first_half = pos_idx - sos_idx;
    second_half = eos_idx - pos_idx;
    growth_rate = (pos_value - location_data(sos_idx)) / first_half;
    senescence_rate = (location_data(eos_idx) - pos_value) / second_half;
    
    % Store the phenology parameters
    phenology_params(i, :) = [pos_value, max_sum, avg_sum, base, duration, first_half, second_half, growth_rate, senescence_rate, num_peaks];
    phenology_dates(i, :) = [dates(sos_idx), dates(pos_idx), dates(eos_idx)];
end

The code calculates various phenology parameters such as POS value, MAX sum, AVG sum, SOS, POS, EOS, Base, Duration, First half, Second half, Growth rate, Senescence rate, and Number of peaks for each location. These parameters are stored in the phenology_params array for future use or analysis.

Thanks @Devendra!

24 个评论
显示 22更早的评论隐藏 22更早的评论

Manikanta Aditya 2024-4-20

在 MATLAB Online 中打开

% You can add your files here like .tif
files = dir('*.tif');
numFiles = length(files);
% Number of bands to stack
numBands = 6;
% Preallocate a 3D matrix to hold the data
stacked_image = zeros(200, 200, numBands, 'uint16');
% Loop over the files in groups of six
for i = 1:numBands:numFiles
    % Reset the stacked image for each group of six
    stacked_image = zeros(200, 200, numBands, 'uint16');
    for j = 0:(numBands-1)
        if (i + j) <= numFiles
            % Read each file
            [img, R] = geotiffread(files(i + j).name);
            
            % Resample the 10m bands to 20m if necessary
            if ismember(files(i + j).name, {'green.tif', 'red.tif', 'nir.tif'})
                img = imresize(img, 0.5, 'Method', 'bilinear');
            end
            
            % Stack the image
            stacked_image(:,:,j+1) = img;
        end
    end
    % Now stacked_image contains your layer stacked image for each group of six
    % You can save each stacked_image to a file here if needed
end

Manikanta Aditya 2024-4-20

在 MATLAB Online 中打开

The code is escaping the resampling part because the condition in the if statement is not being met. The ismember function checks if any of the file names in the cell array {'*_band3_10m.jp2', '*_band4_10m.jp2', *_band8_10m.jp2'} match the current file name files(i + j).name.

However, based on the example file names you provided earlier (T43RGN_20210119T053141_B03_10m.jp2, T43RGN_20210119T053141_B04_10m.jp2, etc.), it seems that the pattern does not match the specified pattern in the if condition.

To fix this, you should modify the cell array in the ismember function to match the naming convention of your actual file names.

For example, if your file names follow the pattern T43RGN_20210119T053141_B03_10m.jp2, you can use the following condition:

if ismember(files(i + j).name, {'*_B03_10m.jp2', '*_B04_10m.jp2', '*_B08_10m.jp2'})
    img = imresize(img, 0.5, 'Method', 'bilinear');
end

This should match the file names with the pattern *_B03_10m.jp2, *_B04_10m.jp2, and *_B08_10m.jp2, and resample those bands to 20m resolution using imresize.

Make sure that the dimensions of the stacked_image variable are correct for the desired output resolution. If you want the output to be at 20m resolution, you should preallocate stacked_image with the appropriate dimensions (5490 x 5490 x 6 in your case).

Devendra 2024-4-24

Hello Aditya,

I am bit hesitant to request you once again, however, I seek your help to write 6bands layer stacked generated by your code to ENVI file. I have tried it very hard but could not succeed. I humbly request you to have a look on customized enviwrite matlab code and suggest me how to make it workable something like similar to attached envi header file. I am using following matlab lines for writing envi format file;

stacked_image(:,:,j+1) = data;

fname = convertStringsToChars(string(img_names) + '_' + '6Bands.dat');

coordinate_information = R; % This assumes R has the needed referencing info

wavelength_data = ['560.00', '665.00', '705.00', '740.00', '783.00', '842.00']

%Writing the ENVI file

custom_enviwrite(stacked_image, fname, coordinate_information, wavelength_data);

but this code is giving following errors

Error in custom_enviwrite (line 69)

map_info.projection, map_info.pixel_size(1), ...

I am attaching custom_enviwrite code alongwith coordinate_information and want to generate envi header file similar to attached envi header file.

Can you please suggest me how to customize custom_enviwrite matlab code to read coordinate_information and generate header file something like attached envi header file?

I am very much grateful to you for your generous support and kind help.

devendra

Manikanta Aditya 2024-4-25

在 MATLAB Online 中打开

@Devendra, I feel that the issue is related to the structure of the 'coordinate_information' variable, which is expected to have specific fields by the 'custom_enviwrite' function.

The 'custom_enviwrite' function expects the 'coordinate_information' to be a struct with two fields: 'map_info' and 'projection_info'. The 'map_info' field should be a string containing the map information in a specific format, and the 'projection_info' field should be a string containing the projection information.

In your code, you're passing R as the coordinate_information argument. Assuming R is a matlab.maps.rasterref.MapCellsReference object, you can extract the required information from it and create the 'coordinate_information' struct as follows:

% Extract map information
map_info = sprintf('UTM, %f, %f, %f, %f, %f, %f, %d, %s, %s, units=Meters', ...
    R.XIntrinsicLimits(1), R.YIntrinsicLimits(1), ...
    R.XWorldLimits(1), R.YWorldLimits(2), ...
    R.CellExtentInWorldX, R.CellExtentInWorldY, ...
    R.ZoneNumber, R.ColumnsStartFrom, R.Datum.Name);
% Extract projection information
projection_info = sprintf('%s', R.ProjectedCRS);
% Create coordinate_information struct
coordinate_information = struct('map_info', map_info, 'projection_info', projection_info);

After creating the 'coordinate_information' struct, you can pass it to the 'custom_enviwrite' function along with the 'stacked_image' and 'wavelength_data' arguments.

Additionally, you'll need to convert the 'wavelength_data' from a string array to a numeric array before passing it to the 'custom_enviwrite' function:

codewavelength_data = str2double(wavelength_data);

With these changes, the code should work as expected and generate the ENVI file with the appropriate header information.

Here's the updated code that should help you understand:

codestacked_image(:,:,j+1) = data;
fname = convertStringsToChars(string(img_names) + '_' + '6Bands.dat');
% Extract coordinate information from R
map_info = sprintf('UTM, %f, %f, %f, %f, %f, %f, %d, %s, %s, units=Meters', ...
    R.XIntrinsicLimits(1), R.YIntrinsicLimits(1), ...
    R.XWorldLimits(1), R.YWorldLimits(2), ...
    R.CellExtentInWorldX, R.CellExtentInWorldY, ...
    R.ZoneNumber, R.ColumnsStartFrom, R.Datum.Name);
projection_info = sprintf('%s', R.ProjectedCRS);
coordinate_information = struct('map_info', map_info, 'projection_info', projection_info);
wavelength_data = str2double({'560.00', '665.00', '705.00', '740.00', '783.00', '842.00'});
% Writing the ENVI file
custom_enviwrite(stacked_image, fname, coordinate_information, wavelength_data);

Please note that it assumes R is a valid 'matlab.maps.rasterref.MapCellsReference' object and 'img_names' is a valid string or character array representing the base filename for the output ENVI file.

Manikanta Aditya 2024-4-25

在 MATLAB Online 中打开

Based on the error messages and the properties you've provided, it seems that R does not have the ZoneNumber and Datum.Name properties. To resolve these issues, we can modify the code to extract the available information from R and format the map_info string accordingly. Additionally, we can handle the projection_info string in a different way.

Check this:

% Extract map information
map_info = sprintf('UTM, %f, %f, %f, %f, %f, %f, %s, %s, units=Meters', ...    
R.XIntrinsicLimits(1), R.YIntrinsicLimits(1), ...    
R.XWorldLimits(1), R.YWorldLimits(2), ...    
R.CellExtentInWorldX, R.CellExtentInWorldY, ...    
R.ColumnsStartFrom, R.CoordinateSystemType);
% Extract projection information
projection_info = strtrim(mat2str(R.ProjectedCRS.Properties.VectorValues));
% Create coordinate_information struct
coordinate_information = struct('map_info', map_info, 'projection_info', projection_info);
% Write the ENVI file
stacked_image(:,:,j+1) = data;
fname = convertStringsToChars(string(img_names) + '_' + '6Bands.dat');
wavelength_data = str2double({'560.00', '665.00', '705.00', '740.00', '783.00', '842.00'});
custom_enviwrite(stacked_image, fname, coordinate_information, wavelength_data);

Manikanta Aditya 2024-4-25

编辑：Manikanta Aditya 2024-4-25

在 MATLAB Online 中打开

% Extract map information
map_info = sprintf('UTM, %f, %f, %f, %f, %f, %f, %s, %s, units=Meters', ...
    R.XIntrinsicLimits(1), R.YIntrinsicLimits(1), ...
    R.XWorldLimits(1), R.YWorldLimits(2), ...
    R.CellExtentInWorldX, R.CellExtentInWorldY, ...
    R.ColumnsStartFrom, R.CoordinateSystemType);
% Extract projection information
projection_info = char(R.ProjectedCRS);
% Create coordinate_information struct
coordinate_information = struct('map_info', map_info, 'projection_info', projection_info);
% Write the ENVI file
stacked_image(:,:,j+1) = data;
fname = convertStringsToChars(string(img_names) + '_' + '6Bands.dat');
wavelength_data = str2double({'560.00', '665.00', '705.00', '740.00', '783.00', '842.00'});
custom_enviwrite(stacked_image, fname, coordinate_information, wavelength_data);

It seems that the ProjectedCRS object in your R variable does not have a Properties field with VectorValues. Tried a different approach to get the projection information from R.ProjectedCRS.

Devendra 2024-4-25

Aditya,

I could retrieve the projection_info as follows

projection_info: "PROJCRS["WGS 84 / UTM zone 43N",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["UTM zone 43N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",75,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Engineering survey, topographic mapping."],AREA["Between 72°E and 78°E, northern hemisphere between equator and 84°N, onshore and offshore. China. India. Kazakhstan. Kyrgyzstan. Maldives. Pakistan. Russian Federation. Tajikistan."],BBOX[0,72,84,78]],ID["EPSG",32643]]"

But I need only following information from above mentioned projection_info:

coordinate system string = {PROJCS["WGS_1984_UTM_Zone_43N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",75.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]}

Now I request you to please help me in retrieving the coordinate system string out of projection_info.

Thanks a lot for helping me.

Devendra

Manikanta Aditya 2024-4-25

在 MATLAB Online 中打开

% Extract map information
map_info = sprintf('UTM, %f, %f, %f, %f, %f, %f, %s, %s, units=Meters', ...
    R.XIntrinsicLimits(1), R.YIntrinsicLimits(1), ...
    R.XWorldLimits(1), R.YWorldLimits(2), ...
    R.CellExtentInWorldX, R.CellExtentInWorldY, ...
    R.ColumnsStartFrom, R.CoordinateSystemType);
% Extract projection information
projection_info = char(R.ProjectedCRS);
% Extract the required coordinate system string
pattern = 'PROJCS\["(.*?)"\]';
matches = regexp(projection_info, pattern, 'tokens');
if ~isempty(matches)
    coordinate_system_string = ['coordinate system string = {' matches{1}{1} '}'];
else
    coordinate_system_string = '';
end
% Create coordinate_information struct
coordinate_information = struct('map_info', map_info, 'projection_info', coordinate_system_string);
% Write the ENVI file
stacked_image(:,:,j+1) = data;
fname = convertStringsToChars(string(img_names) + '_' + '6Bands.dat');
wavelength_data = str2double({'560.00', '665.00', '705.00', '740.00', '783.00', '842.00'});
custom_enviwrite(stacked_image, fname, coordinate_information, wavelength_data);

gauri 2024-4-26

Aditya, I want to extract full file name using following lines

for example file name is T43RGN_20210119T053141_B03_10m.jp2

fn = {files(i + j).name};

img_names = regexp(fn,'_([\dT]+)_','tokens','once');

img_names = [img_names{:}];

and want to store in a variable named as band_names. I request you to please suggest me how to extract full name of files and define the variable band_names to store the full file names.

Thank you for your support and help.

Devendra

请先登录，再进行评论。