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

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.

[Manikanta Aditya]

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

[Devendra]

I am highly obliged to you for your kind help. May I further request for little modifications in the code. I want to use it with the conditions that height of peaks ( value of NDVI should be greater than 0.4 ) should be greater than 0.4 and distance between peaks should be greater than 60 days. I want output as a vector containing the number of peaks only.

I really appreciate your kind support.

Deva

[Manikanta Aditya]

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]

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

[Devendra]

I am attaching the sample file as an output.

[Manikanta Aditya]

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!

[Devendra]

Thank you so much for your kind and generous support. I will accept the answer as soon as you post it.

I am very grateful to you.🙏🙏

Devendra

[Manikanta Aditya]

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

[Devendra]

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]

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!

[Devendra]

Thank you very much for your help and support. It worked very well.

Devendra

[Manikanta Aditya]

Great @Devendra! Thank you.

[Devendra]

I want to layer stack six sentinel2 bands three 10m resolution (green,red and nir) and three 20m resolution (rededge5, rededge6 and rededge 7) through a metlab code in for loop. Because I have 200 files so I want to layer stack these six bands at 20 m resolution using for loop in matlab code. I request you to please provide me some starting point and suggest me how to do it using matlab code.

I always appreciate your kind cooperation.

Devendra

[Manikanta Aditya]

% You can add yout files here like .tif
files = dir('*.tif');

% Preallocate a 3D matrix to hold the data
stacked_image = zeros(200, 200, 6, 'uint16');

for k = 1:length(files)
    % Read each file
    [img, R] = geotiffread(files(k).name);
    
    % Resample the 10m bands to 20m if necessary
    if ismember(files(k).name, {'green.tif', 'red.tif', 'nir.tif'})
        img = imresize(img, 0.5, 'Method', 'bilinear');
    end
    
    % Stack the image
    stacked_image(:,:,k) = img;
end

% Now stacked_image contains your layer stacked image

A basic example for your requirement, hope this helps to provide initial start.

[Devendra]

Thank you very much for your prompt help. I have run the code and following errors have occured;

'C:\working_bpt\bands\'

T43RGN_20210119T053141_B03_10m.jp2

T43RGN_20210119T053141_B04_10m.jp2

T43RGN_20210119T053141_B05_20m.jp2

Unable to perform assignment because the size of the left side is 10980-by-10980 and the size of the right side is 5490-by-5490.

stacked_image(:,:,k) = img;

Since the three bands (band3,band4 and band8) have 10m resolution (10980,10980) and three bands (band5,band6 and band7) have 20m resolution (5490,5490). I want to resample 10m bands into 20m bands resolution so that staced_image should have dimensions of (5490,5490,6). Also I want to write it as envi file using enviwrite. Please suggest me how to do it.

I am very thankful to your generous support and help.

Devendra

[Manikanta Aditya]

% 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

[Devendra]

It is escaping the following part of code

% Resample the 10m bands to 20m if necessary

if ismember(files(i + j).name, {'green.jp2', 'red.jp2', 'nir.jp2'})

img = imresize(img, 0.5, 'Method', 'bilinear');

end

and giving the following error

'C:\working_bpt\bands\'

Unable to perform assignment because the size of the left side is 200-by-200 and the size of the right side is

10980-by-10980.

Error in layer_stack (line 24)

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

I request you to please fix this error also.

Devendra

[Manikanta Aditya]

@Devendra, You should also try from your end rather than relying completely here for the answers, if you don't try you won't learn.

[Devendra]

I am very sorry for asking too much. I will keep in mind this for future. Perhaps this may be my last request for this code. Only thing I want to request you that why it is escaping the following part of code

% Resample the 10m bands to 20m if necessary

if ismember(files(i + j).name, {'*_band3_10m.jp2', '*_band4_10m.jp2', *_band8_10m.jp2'})

img = imresize(img, 0.5, 'Method', 'bilinear');

end

I am once again very sorry.

Devendra

[Manikanta Aditya]

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]

Thanks a lot for your kind help. I used

if ismember(files(i + j).name, {'*_B03_10m.jp2', '*_B04_10m.jp2', '*_B08_10m.jp2'})

img = imresize(img, 0.5, 'Method', 'bilinear');

end

but It gave the following error

Unable to perform assignment because the size of the left side is 5490-by-5490 and the size of the right side is 10980-by-10980.

However, I did little research and used the following lines

if endsWith(files(i+j).name, '_10m.jp2', IgnoreCase=true);

img = imresize(img, 0.5, 'Method', 'bilinear');

end

it worked well. Thanks a lot for your kind help and sorry for troubling you frequently. I will avoid it in future.

Devendra

[Manikanta Aditya]

No problem @Devendra, Thanks!

[Devendra]

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]

@Devendra, I feel you can create another question with your queries, rather than making this question filled with comments..

[Manikanta Aditya]

@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.

[Devendra]

Thank you so much for your kind help. It has almost done the work. But a small error has occured which is as follows;

Unrecognized method, property, or field 'ZoneNumber' for class '

I have observed that MapCellsReference with properties: does not contain these two variables R.ZoneNumber, R.Datum.Name Therefore error has occured. If I remove these two variables than another error has occured as follows

Error using sprintf

Conversion to text from projcrs is not possible.

Error in layer_stack (line 41)

projection_info = sprintf('%s', R.ProjectedCRS);

May I request you to please suggest me how to resolve these errors. I am very much thankful to you.

Devendra

[Manikanta Aditya]

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);

[Devendra]

Now It has given error on following line

Unrecognized method, property, or field 'Properties' for class 'projcrs'.

Error in layer_stack (line 40)

projection_info = strtrim(mat2str(R.ProjectedCRS.Properties.VectorValues));

I request you to please have look on it and suggest me how to go about it.

Devendra

[Manikanta Aditya]

% 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]

It has again given the error as follows;

Error using char

Conversion to char from projcrs is not possible.

Error in layer_stack (line 40)

projection_info = char(R.ProjectedCRS);

Please leave it as I have already taken your considerable time.

Thank you very much for giving me a workable code. I will try to understand it first and than will try to resolve the issues.

Once again thanks a lot for your kind help.

Devendra

[Manikanta Aditya]

As I do not have all the data, it might not able to resolve all the issues, so please check at your end, you should be able to fix the issues.

[Devendra]

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]

% 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]

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