How do I interpolate and smoothen ndvi time series?

I have ununiformly distributed noisey ndvi time series. I want to interpolate and smoothen it using matlab. I am attaching the sample ndvi file and request you to please suggest me how to achieve it using matlab. I would appreciate your kind cooperation and help.

Jyoti

[John D'Errico]

Assuming you mean each row of that table is a distinct time series, sampled as a function of the first row,

The signal to noise ratio appears to be low, at least if you look at any individual series.

The data series are VERY short, which makes it difficult to extract any signal.

There is one data point far away from the rest in x.

At the same time, ALL of those series have a very similar character.

data = readtable('sample_ndvi.csv');
t = table2array(data(1,2:end));
z = table2array(data(2:end,2:end));
surf(z)
xlabel 'Time'
ylabel 'Series'

I've left the time axis as a simple index, since that would just make it all far more dificult to visualize.

The point is, there appears to be a lot of signal, once you look at all of the sereis at once. They all have the same shapes, apparently very noisy. So exactly what smoothing would you want to do?

[Devendra]

Thanks a lot for your suggestions. I want to interpolate the time series over cloudy pixels. Please suggest me how to do interpolation over cloudy pixels. There is a clear dip in the month of august due to cloud. so I want to fill it with neighoubring pixel values.

Jyoti

Answer 1

[Mathieu NOE]

hello

so I simply replace the 8th column with NaN and replaced them using fillmissing2

data = readtable('sample_ndvi.csv');
month = table2array(data(1,2:end));
z = table2array(data(2:end,2:end));
% surf(z)
% xlabel 'Time'
% ylabel 'Series'

% replace august (8th col of z) using fill missing 
z(:,8) = NaN;
zz = fillmissing2(z,'cubic');
surf(zz)
xlabel 'Time'
ylabel 'Series'

[Devendra]

Thank you very much for your help. It serves my purpose. But I have another case where NaN values are palced scattered. How to interploate over these NaN values? I am attaching the file and request you to please suggest me how to fill NaN values in this file.

Devendra

[Mathieu NOE]

the fillmissing2 function handles this case without issues

data = readtable('Test_NaN_fill.csv')

data = 18x21 table
              Var1                 Var2         Var3         Var4         Var5          Var6          Var7          Var8          Var9         Var10         Var11         Var12         Var13         Var14         Var15         Var16         Var17        Var18        Var19        Var20        Var21  
    _________________________    _________    _________    _________    _________    __________    __________    __________    __________    __________    __________    __________    __________    __________    __________    __________    _________    _________    _________    _________    _________

    {0x0 char               }    2.022e+07    2.022e+07    2.022e+07    2.022e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.023e+07    2.023e+07    2.023e+07    2.023e+07    2.023e+07
    {0x0 char               }    2.022e+07    2.022e+07    2.022e+07    2.022e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.0221e+07    2.023e+07    2.023e+07    2.023e+07    2.023e+07    2.023e+07
    {'BAGHPAT.shp'          }          NaN          NaN          NaN          NaN           NaN           NaN           NaN           NaN           NaN           NaN           NaN       0.57919           NaN           NaN           NaN      0.58079      0.55945        0.567      0.59541      0.43577
    {'BAGHPAT.shp'          }          NaN          NaN          NaN          NaN           NaN           NaN           NaN           NaN           NaN           NaN           NaN       0.57919           NaN           NaN           NaN      0.58079      0.55945        0.567      0.59541      0.43577
    {'Bali.shp'             }      0.42776      0.42282      0.44293      0.44103       0.42003        0.4488       0.45565       0.48285       0.43426       0.53449       0.52491       0.52155           NaN       0.47044       0.47355      0.54379      0.52442      0.53578      0.54182      0.47255
    {'Bali.shp'             }      0.42776      0.42282      0.44293      0.44103       0.42003        0.4488       0.45565       0.48285       0.43426       0.53449       0.52491       0.52155           NaN       0.47044       0.47355      0.54379      0.52442      0.53578      0.54182      0.47255
    {'Faizullapur.shp'      }      0.41732          NaN      0.42442       0.4252           NaN       0.43723       0.44024       0.46608       0.42663       0.51161       0.48422        0.5215           NaN       0.43098       0.44231      0.51243      0.49729      0.50531      0.50813      0.43763
    {'Faizullapur.shp'      }      0.41732          NaN      0.42442       0.4252           NaN       0.43723       0.44024       0.46608       0.42663       0.51161       0.48422        0.5215           NaN       0.43098       0.44231      0.51243      0.49729      0.50531      0.50813      0.43763
    {'Gandhi_Urf_Gyasri.shp'}      0.41803      0.42586      0.43649      0.44353           NaN       0.44155       0.44126       0.45301       0.42011       0.51484       0.53512       0.53423       0.40837       0.45641       0.45891       0.5148      0.49989      0.50717      0.51661      0.45093
    {'Gandhi_Urf_Gyasri.shp'}      0.41803      0.42586      0.43649      0.44353           NaN       0.44155       0.44126       0.45301       0.42011       0.51484       0.53512       0.53423       0.40837       0.45641       0.45891       0.5148      0.49989      0.50717      0.51661      0.45093
    {'Kanoli.shp'           }      0.42262      0.42636      0.41273      0.40906           NaN       0.43418       0.43915        0.4634       0.43108       0.54004       0.57082        0.5348           NaN       0.46595       0.47865      0.53911      0.51253      0.52782      0.53632      0.47902
    {'Kanoli.shp'           }      0.42262      0.42636      0.41273      0.40906           NaN       0.43418       0.43915        0.4634       0.43108       0.54004       0.57082        0.5348           NaN       0.46595       0.47865      0.53911      0.51253      0.52782      0.53632      0.47902
    {'Kheriki.shp'          }      0.41379      0.43573      0.43066      0.41393           NaN       0.44298       0.43735       0.44581       0.40941       0.50032       0.47033       0.52479           NaN       0.44489       0.45295      0.52205      0.50205      0.51063      0.51688      0.43939
    {'Kheriki.shp'          }      0.41379      0.43573      0.43066      0.41393           NaN       0.44298       0.43735       0.44581       0.40941       0.50032       0.47033       0.52479           NaN       0.44489       0.45295      0.52205      0.50205      0.51063      0.51688      0.43939
    {'Meetli.shp'           }      0.43394      0.44313       0.4502        0.451           NaN       0.46058       0.45813       0.47756       0.43085       0.52023       0.54378       0.52733       0.41655       0.45888       0.47546      0.54075      0.52637      0.53585      0.54656      0.47659
    {'Meetli.shp'           }      0.43394      0.44313       0.4502        0.451           NaN       0.46058       0.45813       0.47756       0.43085       0.52023       0.54378       0.52733       0.41655       0.45888       0.47546      0.54075      0.52637      0.53585      0.54656      0.47659

z = table2array(data(3:end,2:end))

z = 16x20

NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.5792 NaN NaN NaN 0.5808 0.5595 0.5670 0.5954 0.4358 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.5792 NaN NaN NaN 0.5808 0.5595 0.5670 0.5954 0.4358 0.4278 0.4228 0.4429 0.4410 0.4200 0.4488 0.4557 0.4828 0.4343 0.5345 0.5249 0.5216 NaN 0.4704 0.4736 0.5438 0.5244 0.5358 0.5418 0.4725 0.4278 0.4228 0.4429 0.4410 0.4200 0.4488 0.4557 0.4828 0.4343 0.5345 0.5249 0.5216 NaN 0.4704 0.4736 0.5438 0.5244 0.5358 0.5418 0.4725 0.4173 NaN 0.4244 0.4252 NaN 0.4372 0.4402 0.4661 0.4266 0.5116 0.4842 0.5215 NaN 0.4310 0.4423 0.5124 0.4973 0.5053 0.5081 0.4376 0.4173 NaN 0.4244 0.4252 NaN 0.4372 0.4402 0.4661 0.4266 0.5116 0.4842 0.5215 NaN 0.4310 0.4423 0.5124 0.4973 0.5053 0.5081 0.4376 0.4180 0.4259 0.4365 0.4435 NaN 0.4415 0.4413 0.4530 0.4201 0.5148 0.5351 0.5342 0.4084 0.4564 0.4589 0.5148 0.4999 0.5072 0.5166 0.4509 0.4180 0.4259 0.4365 0.4435 NaN 0.4415 0.4413 0.4530 0.4201 0.5148 0.5351 0.5342 0.4084 0.4564 0.4589 0.5148 0.4999 0.5072 0.5166 0.4509 0.4226 0.4264 0.4127 0.4091 NaN 0.4342 0.4392 0.4634 0.4311 0.5400 0.5708 0.5348 NaN 0.4659 0.4786 0.5391 0.5125 0.5278 0.5363 0.4790 0.4226 0.4264 0.4127 0.4091 NaN 0.4342 0.4392 0.4634 0.4311 0.5400 0.5708 0.5348 NaN 0.4659 0.4786 0.5391 0.5125 0.5278 0.5363 0.4790

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

zz = fillmissing2(z,'v4')  % need the v4  option to get rid of all NaN's

zz = 16x20

0.3941 0.4086 0.4197 0.4265 0.4327 0.4446 0.4599 0.4755 0.4966 0.5313 0.5642 0.5792 0.5724 0.5617 0.5706 0.5808 0.5595 0.5670 0.5954 0.4358 0.4142 0.4213 0.4322 0.4336 0.4306 0.4443 0.4593 0.4708 0.4755 0.5250 0.5568 0.5792 0.5480 0.5188 0.5340 0.5808 0.5595 0.5670 0.5954 0.4358 0.4278 0.4228 0.4429 0.4410 0.4200 0.4488 0.4557 0.4828 0.4343 0.5345 0.5249 0.5216 0.5076 0.4704 0.4736 0.5438 0.5244 0.5358 0.5418 0.4725 0.4278 0.4228 0.4429 0.4410 0.4200 0.4488 0.4557 0.4828 0.4343 0.5345 0.5249 0.5216 0.4976 0.4704 0.4736 0.5438 0.5244 0.5358 0.5418 0.4725 0.4173 0.4198 0.4244 0.4252 0.4267 0.4372 0.4402 0.4661 0.4266 0.5116 0.4842 0.5215 0.4845 0.4310 0.4423 0.5124 0.4973 0.5053 0.5081 0.4376 0.4173 0.4212 0.4244 0.4252 0.4331 0.4372 0.4402 0.4661 0.4266 0.5116 0.4842 0.5215 0.4609 0.4310 0.4423 0.5124 0.4973 0.5053 0.5081 0.4376 0.4180 0.4259 0.4365 0.4435 0.4439 0.4415 0.4413 0.4530 0.4201 0.5148 0.5351 0.5342 0.4084 0.4564 0.4589 0.5148 0.4999 0.5072 0.5166 0.4509 0.4180 0.4259 0.4365 0.4435 0.4422 0.4415 0.4413 0.4530 0.4201 0.5148 0.5351 0.5342 0.4084 0.4564 0.4589 0.5148 0.4999 0.5072 0.5166 0.4509 0.4226 0.4264 0.4127 0.4091 0.4254 0.4342 0.4392 0.4634 0.4311 0.5400 0.5708 0.5348 0.4695 0.4659 0.4786 0.5391 0.5125 0.5278 0.5363 0.4790 0.4226 0.4264 0.4127 0.4091 0.4220 0.4342 0.4392 0.4634 0.4311 0.5400 0.5708 0.5348 0.4939 0.4659 0.4786 0.5391 0.5125 0.5278 0.5363 0.4790

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure(1)
surf(zz)
xlabel 'Time'
ylabel 'Series'

% "view" from the top
figure(2)
imagesc(zz)

[Mathieu NOE]

for those who have not access to fillmissing2, there is an alternative on the Fex :

Inpaint over missing data in 1-D, 2-D, 3-D,... ND arrays - File Exchange - MATLAB Central (mathworks.com)

data = readtable('Test_NaN_fill.csv');
z = table2array(data(5:end,2:end));
% zz = fillmissing2(z,'cubic');
zz = inpaintn(z); % FEX : https://fr.mathworks.com/matlabcentral/fileexchange/27994-inpaint-over-missing-data-in-1-d-2-d-3-d-nd-arrays?s_tid=srchtitle
figure(1)
surf(zz)
xlabel 'Time'
ylabel 'Series'

% "view" from the top
figure(2)
imagesc(zz)

[Mathieu NOE]

and you could use also smoothn (remove's NaN's and smooth the data - as it names suggest ! )

smoothn - File Exchange - MATLAB Central (mathworks.com)

data = readtable('Test_NaN_fill.csv');
z = table2array(data(3:end,2:end));
zz = smoothn(z); 
figure(1)
surf(zz)
xlabel 'Time'
ylabel 'Series'

% "view" from the top
figure(2)
imagesc(zz)

[Devendra]

Thank you so much for your generous support and help. This is what I was looking for. Only one last thing is to write interpolated data in excel csv file with same header (first and second rows) and same first column (name of shape files) as attached input csv file. I request you to please suggest me how to put this data into a excel csv file similar to attached input file. I really appriciate your kind help.

Devendra

[Mathieu NOE]

before I do that I have one more suggestion (John would have been unhappy if I had forgotten about his work !!)

inpaint_nans - File Exchange - MATLAB Central (mathworks.com)

data = readtable('Test_NaN_fill.csv');
z = table2array(data(3:end,2:end));
zz = inpaint_nans(z,3);
figure(1)
surf(zz)
xlabel 'Time'
ylabel 'Series'

% "view" from the top
figure(2)
imagesc(zz)

[Devendra]

Thanks a lot for more functions to be used. I will download them and try them one by one.

Devendra

[Mathieu NOE]

here's the full code with export to csv file

Maybe you noticed I have made some corrections above in my previous comments , because I was initially only taking rows 5 to end , but the correct rows index is from 3 to end for this csv file.

after this modification, I was surprised to see that fillmissing2 seemed not able to replace some NaN's values on the first lines (boundaries) , whatever the interpolation method I pick amon linear, cubic, natural.

Had to go to "v4" method to get all NaN's replaced. Something I have to remember for the future (or use the alternatives I mentionned above).

data = readtable('Test_NaN_fill.csv');
[m,n] = size(data);
ind_rows = 3:m;
ind_cols = 2:n;
z = table2array(data(ind_rows,ind_cols))

z = 16x20

NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.5792 NaN NaN NaN 0.5808 0.5595 0.5670 0.5954 0.4358 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.5792 NaN NaN NaN 0.5808 0.5595 0.5670 0.5954 0.4358 0.4278 0.4228 0.4429 0.4410 0.4200 0.4488 0.4557 0.4828 0.4343 0.5345 0.5249 0.5216 NaN 0.4704 0.4736 0.5438 0.5244 0.5358 0.5418 0.4725 0.4278 0.4228 0.4429 0.4410 0.4200 0.4488 0.4557 0.4828 0.4343 0.5345 0.5249 0.5216 NaN 0.4704 0.4736 0.5438 0.5244 0.5358 0.5418 0.4725 0.4173 NaN 0.4244 0.4252 NaN 0.4372 0.4402 0.4661 0.4266 0.5116 0.4842 0.5215 NaN 0.4310 0.4423 0.5124 0.4973 0.5053 0.5081 0.4376 0.4173 NaN 0.4244 0.4252 NaN 0.4372 0.4402 0.4661 0.4266 0.5116 0.4842 0.5215 NaN 0.4310 0.4423 0.5124 0.4973 0.5053 0.5081 0.4376 0.4180 0.4259 0.4365 0.4435 NaN 0.4415 0.4413 0.4530 0.4201 0.5148 0.5351 0.5342 0.4084 0.4564 0.4589 0.5148 0.4999 0.5072 0.5166 0.4509 0.4180 0.4259 0.4365 0.4435 NaN 0.4415 0.4413 0.4530 0.4201 0.5148 0.5351 0.5342 0.4084 0.4564 0.4589 0.5148 0.4999 0.5072 0.5166 0.4509 0.4226 0.4264 0.4127 0.4091 NaN 0.4342 0.4392 0.4634 0.4311 0.5400 0.5708 0.5348 NaN 0.4659 0.4786 0.5391 0.5125 0.5278 0.5363 0.4790 0.4226 0.4264 0.4127 0.4091 NaN 0.4342 0.4392 0.4634 0.4311 0.5400 0.5708 0.5348 NaN 0.4659 0.4786 0.5391 0.5125 0.5278 0.5363 0.4790

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

zz = fillmissing2(z,"cubic") % KO

zz = 16x20

NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.5792 0.5807 0.5783 0.5768 0.5808 0.5595 0.5670 0.5954 0.4358 NaN NaN NaN NaN NaN NaN 0.4974 0.5011 0.5162 0.5255 0.5558 0.5792 0.5259 0.5177 0.5200 0.5808 0.5595 0.5670 0.5954 0.4358 0.4278 0.4228 0.4429 0.4410 0.4200 0.4488 0.4557 0.4828 0.4343 0.5345 0.5249 0.5216 0.4938 0.4704 0.4736 0.5438 0.5244 0.5358 0.5418 0.4725 0.4278 0.4228 0.4429 0.4410 0.4200 0.4488 0.4557 0.4828 0.4343 0.5345 0.5249 0.5216 0.4977 0.4704 0.4736 0.5438 0.5244 0.5358 0.5418 0.4725 0.4173 0.4203 0.4244 0.4252 0.4294 0.4372 0.4402 0.4661 0.4266 0.5116 0.4842 0.5215 0.4797 0.4310 0.4423 0.5124 0.4973 0.5053 0.5081 0.4376 0.4173 0.4211 0.4244 0.4252 0.4316 0.4372 0.4402 0.4661 0.4266 0.5116 0.4842 0.5215 0.4712 0.4310 0.4423 0.5124 0.4973 0.5053 0.5081 0.4376 0.4180 0.4259 0.4365 0.4435 0.4434 0.4415 0.4413 0.4530 0.4201 0.5148 0.5351 0.5342 0.4084 0.4564 0.4589 0.5148 0.4999 0.5072 0.5166 0.4509 0.4180 0.4259 0.4365 0.4435 0.4421 0.4415 0.4413 0.4530 0.4201 0.5148 0.5351 0.5342 0.4084 0.4564 0.4589 0.5148 0.4999 0.5072 0.5166 0.4509 0.4226 0.4264 0.4127 0.4091 0.4203 0.4342 0.4392 0.4634 0.4311 0.5400 0.5708 0.5348 0.4920 0.4659 0.4786 0.5391 0.5125 0.5278 0.5363 0.4790 0.4226 0.4264 0.4127 0.4091 0.4204 0.4342 0.4392 0.4634 0.4311 0.5400 0.5708 0.5348 0.4991 0.4659 0.4786 0.5391 0.5125 0.5278 0.5363 0.4790

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

zz = fillmissing2(z,"v4") % OK

zz = 16x20

0.3941 0.4086 0.4197 0.4265 0.4327 0.4446 0.4599 0.4755 0.4966 0.5313 0.5642 0.5792 0.5724 0.5617 0.5706 0.5808 0.5595 0.5670 0.5954 0.4358 0.4142 0.4213 0.4322 0.4336 0.4306 0.4443 0.4593 0.4708 0.4755 0.5250 0.5568 0.5792 0.5480 0.5188 0.5340 0.5808 0.5595 0.5670 0.5954 0.4358 0.4278 0.4228 0.4429 0.4410 0.4200 0.4488 0.4557 0.4828 0.4343 0.5345 0.5249 0.5216 0.5076 0.4704 0.4736 0.5438 0.5244 0.5358 0.5418 0.4725 0.4278 0.4228 0.4429 0.4410 0.4200 0.4488 0.4557 0.4828 0.4343 0.5345 0.5249 0.5216 0.4976 0.4704 0.4736 0.5438 0.5244 0.5358 0.5418 0.4725 0.4173 0.4198 0.4244 0.4252 0.4267 0.4372 0.4402 0.4661 0.4266 0.5116 0.4842 0.5215 0.4845 0.4310 0.4423 0.5124 0.4973 0.5053 0.5081 0.4376 0.4173 0.4212 0.4244 0.4252 0.4331 0.4372 0.4402 0.4661 0.4266 0.5116 0.4842 0.5215 0.4609 0.4310 0.4423 0.5124 0.4973 0.5053 0.5081 0.4376 0.4180 0.4259 0.4365 0.4435 0.4439 0.4415 0.4413 0.4530 0.4201 0.5148 0.5351 0.5342 0.4084 0.4564 0.4589 0.5148 0.4999 0.5072 0.5166 0.4509 0.4180 0.4259 0.4365 0.4435 0.4422 0.4415 0.4413 0.4530 0.4201 0.5148 0.5351 0.5342 0.4084 0.4564 0.4589 0.5148 0.4999 0.5072 0.5166 0.4509 0.4226 0.4264 0.4127 0.4091 0.4254 0.4342 0.4392 0.4634 0.4311 0.5400 0.5708 0.5348 0.4695 0.4659 0.4786 0.5391 0.5125 0.5278 0.5363 0.4790 0.4226 0.4264 0.4127 0.4091 0.4220 0.4342 0.4392 0.4634 0.4311 0.5400 0.5708 0.5348 0.4939 0.4659 0.4786 0.5391 0.5125 0.5278 0.5363 0.4790

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure(1)
surf(zz)
xlabel 'Time'
ylabel 'Series'

% "view" from the top
figure(2)
imagesc(zz)

% create new table and save to csv file
out = data;
out{ind_rows,ind_cols} = zz;
writetable(out,"file_out.csv", "WriteVariableNames", 0);
type file_out.csv

,20220404,20220409,20220414,20220419,20220504,20220514,20220519,20220603,20220613,20220623,20220708,20220812,20221031,20221125,20221205,20230203,20230208,20230213,20230223,20230315
,20220404,20220409,20220414,20220419,20220504,20220514,20220519,20220603,20220613,20220623,20220708,20220812,20221031,20221125,20221205,20230203,20230208,20230213,20230223,20230315
BAGHPAT.shp,0.394120111555345,0.408639344817413,0.419692857853629,0.426475071497293,0.432721203362728,0.444615248599096,0.459896693508927,0.475512581971739,0.496589202846398,0.531329944271754,0.564231755930429,0.5791876,0.572428826521747,0.561685851267669,0.570554011779578,0.5807873,0.5594513,0.5670022,0.5954056,0.4357717
BAGHPAT.shp,0.414215965590067,0.421340854618457,0.432235430336945,0.433568889318686,0.430578937373586,0.444301288238815,0.459294249750126,0.470772737737899,0.475485503221496,0.524980588649825,0.556804548381212,0.5791876,0.547956631144423,0.518846618762357,0.534043819205039,0.5807873,0.5594513,0.5670022,0.5954056,0.4357717
Bali.shp,0.4277567,0.4228248,0.4429285,0.4410278,0.420029,0.4487961,0.455652,0.4828474,0.434257,0.53449,0.5249099,0.5215526,0.507562787879291,0.4704439,0.4735548,0.5437946,0.5244235,0.5357844,0.5418203,0.4725451
Bali.shp,0.4277567,0.4228248,0.4429285,0.4410278,0.420029,0.4487961,0.455652,0.4828474,0.434257,0.53449,0.5249099,0.5215526,0.497631810437344,0.4704439,0.4735548,0.5437946,0.5244235,0.5357844,0.5418203,0.4725451
Faizullapur.shp,0.4173222,0.419780720062812,0.424421,0.4252048,0.426731349026019,0.437231,0.4402393,0.4660785,0.4266327,0.5116112,0.4842201,0.521501,0.484490441172325,0.4309783,0.4423086,0.5124264,0.49729,0.5053054,0.5081333,0.4376335
Faizullapur.shp,0.4173222,0.421162905132942,0.424421,0.4252048,0.433089138844667,0.437231,0.4402393,0.4660785,0.4266327,0.5116112,0.4842201,0.521501,0.460890338587612,0.4309783,0.4423086,0.5124264,0.49729,0.5053054,0.5081333,0.4376335
Gandhi_Urf_Gyasri.shp,0.4180337,0.4258612,0.436492,0.4435274,0.443948523541826,0.4415459,0.4412603,0.4530149,0.4201115,0.5148381,0.5351193,0.5342322,0.4083741,0.4564052,0.4589096,0.5148034,0.4998898,0.5071719,0.516614,0.4509309
Gandhi_Urf_Gyasri.shp,0.4180337,0.4258612,0.436492,0.4435274,0.442186019485276,0.4415459,0.4412603,0.4530149,0.4201115,0.5148381,0.5351193,0.5342322,0.4083741,0.4564052,0.4589096,0.5148034,0.4998898,0.5071719,0.516614,0.4509309
Kanoli.shp,0.4226173,0.4263566,0.4127327,0.4090571,0.425410129124295,0.4341846,0.4391508,0.4634013,0.4310781,0.540037,0.5708241,0.5348016,0.469474836180193,0.4659472,0.4786469,0.5391057,0.512527,0.527822,0.5363246,0.4790224
Kanoli.shp,0.4226173,0.4263566,0.4127327,0.4090571,0.422009921672363,0.4341846,0.4391508,0.4634013,0.4310781,0.540037,0.5708241,0.5348016,0.493944186020775,0.4659472,0.4786469,0.5391057,0.512527,0.527822,0.5363246,0.4790224
Kheriki.shp,0.4137874,0.4357251,0.4306581,0.4139348,0.426151175466462,0.4429752,0.4373541,0.4458094,0.4094112,0.5003205,0.4703299,0.524792,0.493180506523894,0.4448929,0.4529536,0.5220464,0.5020531,0.5106305,0.5168774,0.4393871
Kheriki.shp,0.4137874,0.4357251,0.4306581,0.4139348,0.432442993682095,0.4429752,0.4373541,0.4458094,0.4094112,0.5003205,0.4703299,0.524792,0.472629698421247,0.4448929,0.4529536,0.5220464,0.5020531,0.5106305,0.5168774,0.4393871
Meetli.shp,0.4339411,0.4431321,0.4502024,0.4510042,0.45191527722411,0.4605763,0.4581292,0.4775634,0.4308531,0.5202342,0.5437839,0.527333,0.4165452,0.4588833,0.4754621,0.5407523,0.5263727,0.535851,0.5465643,0.4765886
Meetli.shp,0.4339411,0.4431321,0.4502024,0.4510042,0.446425211968403,0.4605763,0.4581292,0.4775634,0.4308531,0.5202342,0.5437839,0.527333,0.4165452,0.4588833,0.4754621,0.5407523,0.5263727,0.535851,0.5465643,0.4765886
Sankrod.shp,0.4312419,0.4419311,0.4436681,0.4427217,0.4155471,0.4609338,0.4449305,0.4631734,0.471258968663044,0.5220414,0.55035,0.5251712,0.459506601714459,0.4484476,0.4802725,0.5940056,0.5455006,0.5608983,0.546536,0.4597271
Sankrod.shp,0.4312419,0.4419311,0.4436681,0.4427217,0.4155471,0.4609338,0.4449305,0.4631734,0.488253838867251,0.5220414,0.55035,0.5251712,0.476746091080244,0.4484476,0.4802725,0.5940056,0.5455006,0.5608983,0.546536,0.4597271

[Devendra]

Thank you very much for your kind help. I am very much thankful to you for your kind cooperation and support.

Devendra

[Mathieu NOE]

as always, my pleasure !

[Devendra]

May I request you a little more help in the following error;

Actually I have defined the

coordinate_information = struct('map_info', map_info, 'projection_info', projection_info);

I am writing it as follows

fprintf(fid, coordinate_information);

than it is giving following error as follows

Error in custom_enviwrite (line 68)

Error using fprintf

Invalid format.

fprintf(fid, coordinate_information);

I request you to please suggest me how to write it in correct way. I would appreciate your kind cooperation and help.

Devendra

[Mathieu NOE]

hello again

I am not sure what you are trying to do here - can you explain a bit more ?

do you want to display or store your coordinate_information struct to a file ?

why not export as a cell array using writecell ?

[Mathieu NOE]

I don't know if this is your goal, but if you want to export nicely formated tables to txt files , I have this (a bit old but still working)

see the 2 functions in attachment

results in txt file containing :

COne CTwo CThree CFour CFive

ROne 9.5929E-01 8.4072E-01 3.4998E-01 3.5166E-01 2.8584E-01

RTwo 5.4722E-01 2.5428E-01 1.9660E-01 8.3083E-01 7.5720E-01

RThree 1.3862E-01 8.1428E-01 2.5108E-01 5.8526E-01 7.5373E-01

RFour 1.4929E-01 2.4352E-01 6.1604E-01 5.4972E-01 3.8045E-01

RFive 2.5751E-01 9.2926E-01 4.7329E-01 9.1719E-01 5.6782E-01

%demo code for MyWriteTable
A = rand(5,5);
CHeaders = 'COne CTwo CThree CFour CFive'; % columns headers 
RHeaders = 'ROne RTwo RThree RFour RFive'; % rows headers
MyWriteTable('TableData.txt',A,'13.4E','w',Makemat(CHeaders),Makemat(RHeaders))

[Devendra]

I want to store coordinate_information struct to a file. I am using the following matlab line to write

fprintf(fid,'%s \n','interleave = bsq');

% Additional coordinate information

fprintf(fid,'%g \n',coordinate_information); % to write struct coordinate_information to a file.

Please suggest me how do I write it in a file.

Thank you very much for your kind help.

Devendra

[Mathieu NOE]

ok but I need to have map_info and projection_info (what is in your struct array)

can you share it ? you can save it in a mat file

[Devendra]

coordinate_information = struct('map_info', map_info, 'projection_info', projection_info);

disp(coordinate_information);

map_info: 'UTM, 0.500000, 0.500000, 699960.000000, 3300000.000000, 10.000000, 10.000000, north, planar, units=Meters'

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

Actually 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]]}

I want to write only above mentioned coordinate system string out of projection_info.

This is the information to be written in the file.

Thanks a lot for helping me.

Devendra

[Mathieu NOE]

sorry but I am still confused

so "coordinate system string" is comint out of "projection_info" ? I can't find anything in "coordinate system string" that matches with "projection_info"

the final goal is just to write strings (of whatever length ?) to a txt file ?

[Devendra]

yes please. my final goal is just to write strings (of whatever length ) to a txt file. I request you to please suggest me to write coordinate_information = struct('map_info', map_info, 'projection_info', projection_info); in a txt file.

Devendra

[Mathieu NOE]

I don't see the need to create a structure here , if the goal is simply to write strings or char arrays to a txt file

see demo code below

the resulting txt file should look like this

-------

map info :

UTM, 0.500000, 0.500000, 699960.000000, 3300000.000000, 10.000000, 10.000000, north, planar, units=Meters

-------

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

-------

coordinate system :

{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]]}

code :

% some char or string arrays

map_info = 'UTM, 0.500000, 0.500000, 699960.000000, 3300000.000000, 10.000000, 10.000000, north, planar, units=Meters';

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]]"';
 
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]]}';
 
 % Let's say you have a cell array of strings
listOfStrings = {'-------';'map info : ';map_info;'-------';'projection info : '; projection_info;'-------';'coordinate system : '; coordinate_system_string};

% Open a new text file for writing
fileID = fopen('output.txt', 'w');

% Check if the file was opened successfully
if fileID == -1
    error('File could not be opened.');
end

% Write each string in the cell array to the file on a new line
for i = 1:length(listOfStrings)
    fprintf(fileID, '%s\n', listOfStrings{i});
end

% Close the file
fclose(fileID);

[Devendra]

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

Devendra

[Mathieu NOE]

as always , my pleasure !!

[Devendra]

There is a little issue in writing the header files. I am attaching two header files one is created by MATLAB code and second is generated by ENVI. When I open the header file created by MATLAB code in ENVI software a portion of it gets deleted. I am attaching that header file also. What I observed that map info and projection info are in { } brackets for the header file created by ENVI software while the header file created by MATLAB code are not in { }. probably that may be the reason I am still not sure. The header file created by MATLAB is 20210203T053029_6Bands.hdr while header file created by ENVI software is ENVI_header_file

while a portion of header file deleted is 20210119T053141_6Bands.

I request you to please have a look on these files and suggest me how to produce the header file similar to ENVI generated header file.

Thank you very much for your help.

Devendra

[Mathieu NOE]

ok, I believe I understand the problem

so I start with the same 3 strings that do not have any {} brackets , then I add them (around so to say) , then I export to txt file

this is how the result looks like on y side :

result :

-------

map info =

{UTM, 0.500000, 0.500000, 699960.000000, 3300000.000000, 10.000000, 10.000000, north, planar, units=Meters}

-------

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

-------

coordinate system =

{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]]}

code

% some char or string arrays (assuming without {})
map_info = 'UTM, 0.500000, 0.500000, 699960.000000, 3300000.000000, 10.000000, 10.000000, north, planar, units=Meters';
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]]';
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]]';
 
 % Let's say you have a cell array of strings
 % and also I add the {} around each string
listOfStrings = {'-------';'map info = ';['{' map_info '}'];'-------';'projection info = '; ['{' projection_info '}'];'-------';'coordinate system = '; ['{' coordinate_system_string '}']};
% Open a new text file for writing
fileID = fopen('output.txt', 'w');
% Check if the file was opened successfully
if fileID == -1
    error('File could not be opened.');
end
% Write each string in the cell array to the file on a new line
for i = 1:length(listOfStrings)
    fprintf(fileID, '%s\n', listOfStrings{i});
end
% Close the file
fclose(fileID);

[gauri]

I have a MATLAB code for the computation of following phenology parameters;

•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 and chlorophyll ab- sence

•POS value: maximum NDVI value in the time series

•Peaks: number of local maximas

•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

•First half: duration of green-up in days between SOS and POS

•Second half: duration of senescence in days between POS and EOS

•Growth rate: gradient between SOS and POS

•Senescence rate: gradient between POS and EOS

However, I have also put threshold values which are creating problem and code is not executed over entire file. What I wanted is to get rid of these constraints and put some alternative to these constraints.

I want to modify the code for the following things without any numerical values constraints;

First minima before first maxima as SOS

First highest maxima as POS and its value is POS value

Last minima after first maxima and second maxima as EOS OR first minima from the end of NDVI series

Number of peaks are first maxima and may be second maxima

I am attaching MATLAB code and request you to please have a look on it and suggest me how to modify the code for the incorporation of above mentioned conditions.

I would appreciate your kind cooperation and your kind help.

Gauri

[Mathieu NOE]

hello @gauri

it's better not to mix several topics in the same post

by posting seperately you also increases the number of views and potential answers

I'll have a look in the next days , but to work on your code I need some input data (csv file)

all the best

[gauri]

Thank you so much for your kind reply and suggestion. I will keep in mind to post seperately any further questions. However, I am attaching the NDVI file on which I have been working and request you to please have a look on code and data and suggest me how to make code workable over attached ndvi file.

I would be highly obliged to you for your kind help.

Gauri

[Mathieu NOE]

I tried quickly your code but it shows errors on dates

Error in test_phenology2 (line 13)

gdt = datetime(dates,"ConvertFrom","yyyymmdd") % Gregorian

seems the csv file you provided dos not contain any data / info related to dates :

[gauri]

I am sorry for my mistake. I am attaching the ndvi file with date as header.

Gauri

[Mathieu NOE]

hello again

with the provided csv file , I have an error occuring at the 4th iteration of the for loop

Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.

Error in test_phenology (line 59)

phenology_params(i, :) = [sos,pos,eos,...

this is due to findpeaks returning empty output

[pks, locs, ~, prominences] = findpeaks(location_data, 'MinPeakProminence', 0.1, 'MinPeakHeight', 0.5, 'MinPeakDistance', 5);

the corresponding data plot looks like

I need to understand how did you choose the findpeaks parameters and what for ?

what is the target of your code ?

[Mathieu NOE]

if I plot all your data (using imagesc)

imagesc(ndvi_data)
colormap('jet')
colorbar('vert')

we can see either you have one peak (per location) or you have multiple similar amplitude peaks - like in the prevoius comment

what are we supposed to do in that case ?

[gauri]

My target is to extract 3 dates and number of peaks;

First date is SOS(start of season) first minima just before the first maxima

Second date is POS (peak of season) first maxima just after the first minima and POS value is the NDVI value of first maxima

Third date is EOS(end of season) last minima after all maximas (peaks)

Number of peaks : Generally there are two peaks; first peak is first maxima and second peak is second maxima which is seperated by at least 4 to 5 dates in NDVI series.

I hope this will help.

Thanks for your kind help.

Gauri

[Mathieu NOE]

hello again

I looked at all you data , but this file (Sample_ndvi.csv) shows man time curves that are not evident to analyse and will probably fail to go through the tests / criteria you have defined ;

I suppose that "good" curves are when we have two well shaped peaks . Then what do we do in case it's very different ?

[gauri]

Thanks a lot for your valuable observations. I have been now trying to include cloudy images also and interpolating them by using HANTS MATLAB application. I hope this will produce well defined two peaks in my NDVI data. I will request you once again after getting HANTS interpolated NDVI data covering full period of 13 months(two images per months). Now you may please close it for a while.

Once again I am grateful to you for your kind help.

Jyoti

[Mathieu NOE]

ok - let's see what comes next !

have a good day

[gauri]

HANTS MATLAB application is not working on my data due to some defined parameters. I have 26 georeferenced ENVI format satellite image data over a period of 13 months (two image per months). Most of them are clear sky images. However, few of them due to monsson season are cloudy images specially in the month of july, august and september (middle of year). I want to interpolate and smoothen these cloudy images. May I request you to please suggest me how to do it or if you can please provide me some matlab code to do it? I would be highly obliged to you.

The size of images are 1350,1153,7(lines,samples,bands)

Looking forward for your kind suggestions.

Gauri

[Mathieu NOE]

I can try to do something , but you have to find a way to share the images

[Devendra]

I can put the images on google drive and provide you the link. Thanks for your positive response.

[Mathieu NOE]

ok let's do that !

[gauri]

I am pasting the link for data. Please get the data from this link.

https://drive.google.com/drive/folders/10M2JX7bFGrvitfW6kxiKvmfweEDAB-3p?usp=sharing

Answer 2

[gauri]

I am pasting the link for data. Please get the data from this link.

https://drive.google.com/drive/folders/10M2JX7bFGrvitfW6kxiKvmfweEDAB-3p?usp=sharing

[gauri]

Could you please let me know the status of interploation of image data sent to you by above google link?

I appreciate your kind help and support.

[Mathieu NOE]

hello

I have downloaded the files , but how do you load these files in the workspace ?

are you using the Image Toolbox ? Identify Vegetation Regions Using Interactive NDVI Thresholding - MATLAB & Simulink - MathWorks France

[gauri]

I read these files using code as follows;

d='C:\works\' % your directory

files = dir(fullfile(d, '*.dat'));

% Preallocate a 3D matrix to hold the data

stacked_image = zeros(numFiles, 1315, 1153);

numFiles = length(files);

for i = 1:numFiles;

% Read each file

img = hypercube(files(i).name);

% convert integer numbers into real numbers

data = single(img.DataCube);

red = data(:,:,2); % red band

nir = data(:,:,6); % near infrared band

ndvi = (nir - red) ./ (nir + red);

% Stack the image

stacked_image(i, :, :) = ndvi

end

I hope it will help to read these files.

Thanks for your kind help.

[Mathieu NOE]

hello again

that is my trouble, you're using hypercube from the image processing toolbox (I don't have)

maybe you could run your code and save stacked_image array in a mat file and share this one (maybe very big )

[gauri]

I have stored the ndvi data in ndvi.mat file and uploaded it on following link;

https://drive.google.com/drive/folders/10M2JX7bFGrvitfW6kxiKvmfweEDAB-3p?q=after:2024-05-07%20parent:10M2JX7bFGrvitfW6kxiKvmfweEDAB-3p

I request you to please check it and let me know whether it is ok or not.

[Mathieu NOE]

hmmm, it says could not access server , try later ...

does it work on your side ?

[gauri]

Yes it does work on my side. I can download it here.

[Mathieu NOE]

tried again but still have the same error

will try a bit later on

[gauri]

I am again sending the link as follows;

https://drive.google.com/drive/folders/10M2JX7bFGrvitfW6kxiKvmfweEDAB-3p?usp=sharing

I request you to please check this link.

[Mathieu NOE]

it works !

[Mathieu NOE]

so this is a first trial

idea was to replace "bad" images with nans then interpolate the 3D array using either one of these Fex submissions :

inpaint_nans - File Exchange - MATLAB Central (mathworks.com)

Inpaint over missing data in 1-D, 2-D, 3-D,... ND arrays - File Exchange - MATLAB Central (mathworks.com)

for the time being I did a manual selection of which images needed to be interpolated , but maybe a smarter approach is doable

the selection is based on the mean value of each image, we can see the drop in this curve

so we have 24 mean values and the graph looks like :

so I decided I wanted to interpolate image 7, 10,11 and 18

once we have corrected (interpolated) the images, we get a new mean curve :

full code :

load('ndvi.mat')

[m,n,p] = size(stacked_image);
mi = min(stacked_image,[],'all');
ma = max(stacked_image,[],'all');

x_axis = (1:m); % 12 months = 24 images

%% step1 : plot the raw data
for ck = 1:m
    figure
    data2D = squeeze(stacked_image(ck,:,:));
    data2Dmean(ck) = mean(data2D,'all','omitnan');
    imagesc(data2D);
    caxis([mi ma]);
    colorbar('vert');
    
end

% mean 2D raw data plot
figure
plot(x_axis,data2Dmean);

%% step 2 : identify images to be interpolated
ind = [7 10 11 18];% selected image to correct : 

stacked_image2 = stacked_image;
stacked_image2(ind,:,:) = NaN;

% now do interpolation on selected images 
% faster : 
stacked_image2=inpaint_nans3(stacked_image2,1); % see Fex : https://fr.mathworks.com/matlabcentral/fileexchange/21214-inpainting-nan-elements-in-3-d
% slower : 
% stacked_image2=inpaintn(stacked_image2); % see Fex : https://fr.mathworks.com/matlabcentral/fileexchange/27994-inpaint-over-missing-data-in-1-d-2-d-3-d-nd-arrays?s_tid=srchtitle

for ck = 1:m
    figure
    data2D = squeeze(stacked_image2(ck,:,:));
    data2Dmean(ck) = mean(data2D,'all','omitnan');
    imagesc(data2D);
    caxis([mi ma]);
    colorbar('vert');
    
end

% mean 2D interpolated data plot

figure
plot(x_axis,data2Dmean);

[gauri]

Thanks a lot for your kind help. Looking forward for identying the cloudy images automatically. Now I want to write back these corrected image data into a mat file. Shall I write stacked_image2 as final data?.

[Mathieu NOE]

hello again

yes you can now save stacked_image2 as final data (as mat file)

I tried to find a solution for the automatic detection of cloudy image , based on the ratio of mean and std values

seems to work but may be tweaked in other situations

you need the attached function () to run it

new code :

load('ndvi.mat')

[m,n,p] = size(stacked_image);
mi = min(stacked_image,[],'all');
ma = max(stacked_image,[],'all');

x_axis = (1:m); % 12 months = 24 images


%% step1 : plot the raw data
for ck = 1:m
    data2D = squeeze(stacked_image(ck,:,:));
%     figure(ck)
%     imagesc(data2D);
%     caxis([mi ma]);
%     colorbar('vert');
    data2Dmean(ck) = mean(data2D,'all','omitnan');
    data2Dstd(ck) = std(data2D,1,'all','omitnan');
    
end
% normalise mean and std values
data2Dmean = data2Dmean./max(data2Dmean);
data2Dstd = data2Dstd./max(data2Dstd);
r = abs(data2Dstd./data2Dmean); % ratio std over mean

% mean 2D raw data plot
N = m/4; % N must be choosen >= number of consecutive cloudy images (assumed) + 1 
[down] = env_secant(x_axis,r, N, 'bottom');

% figure
% plot(x_axis,r,'-*',x_axis,down,'-*');

% selected indices (one year)
ind = find(r>1.5*down); % selected image to correct (down = threshold curve)

return
%% step 2 : identify images to be interpolated
% ind = [7 10 11 18];% selected image to correct : 

stacked_image2 = stacked_image;
stacked_image2(ind,:,:) = NaN;

% now do interpolation on selected images 
% faster : 
stacked_image2=inpaint_nans3(stacked_image2,1); % see Fex : https://fr.mathworks.com/matlabcentral/fileexchange/21214-inpainting-nan-elements-in-3-d
% slower : 
% stacked_image2=inpaintn(stacked_image2); % see Fex : https://fr.mathworks.com/matlabcentral/fileexchange/27994-inpaint-over-missing-data-in-1-d-2-d-3-d-nd-arrays?s_tid=srchtitle

for ck = 1:m
    figure
    data2D = squeeze(stacked_image2(ck,:,:));
    data2Dmean(ck) = mean(data2D,'all','omitnan');
    imagesc(data2D);
    caxis([mi ma]);
    colorbar('vert');
    
end

% mean 2D interpolated data plot

figure
plot(x_axis,data2Dmean);

[gauri]

I have tried it and following errors have occured;

Error using assert

Parameter <y_data> has to be of vector type, holding finite numeric values!

Error in env_secant (line 25)

assert( isvector(y_data) && isnumeric(y_data) && all( isfinite( y_data ) ), ...

Error in ndvi_smooth (line 61)

[down] = env_secant(x_axis,r, N, 'bottom');

I request to please suggest me how to fix it.

[Mathieu NOE]

hello again

did that occur with the same mat file ? (ndvi.mat) or is with another set of files ?

[gauri]

This occured on a different set of files. I can save it in mat file and put on a google drive if required.

[gauri]

I have put this another set of file (new_ndvi.mat) on following link;

https://drive.google.com/file/d/1cQjW8I4Ak8CwpGpP_hE4W5PvGkofYT7e/view?usp=drive_link

I request you to please have a look on it.

[Mathieu NOE]

hello

yes the problem comes from stacked_image is organized differently in the two provided mat file

ndvi.mat : stacked_image 24x1315x1153

new_ndvi.mat : stacked_image 1315x1153x24

so I had to change the way to index the data vs the image number , like

ndvi.mat : data2D = squeeze(stacked_image(ck,:,:));

new_ndvi.mat : data2D = squeeze(stacked_image(:,:,ck));

updated code :

(I have come back to te manual selection of images as my method is not yet robust enough)

clc
clearvars
close all

load('new_ndvi.mat') % OK

[m,n,p] = size(stacked_image);
mi = min(stacked_image,[],'all');
ma = max(stacked_image,[],'all');

x_axis = (1:p); % 12 months = 24 images


%% step1 : plot the raw data
for ck = 1:p
    data2D = squeeze(stacked_image(:,:,ck));
    figure(ck)
    imagesc(data2D);
    caxis([mi ma]);
    colorbar('vert');
    data2Dmean(ck) = mean(data2D,'all','omitnan');
    data2Dstd(ck) = std(data2D,1,'all','omitnan');
    
end

%% step 2 : identify images to be interpolated
ind = [6];% selected image to correct : 
% ind = [6 7 9 10];% selected image to correct : 
stacked_image2 = stacked_image;
stacked_image2(:,:,ind) = NaN;

% now do interpolation on selected images 
% faster : 
stacked_image2=inpaint_nans3(stacked_image2,1); % see Fex : https://fr.mathworks.com/matlabcentral/fileexchange/21214-inpainting-nan-elements-in-3-d
% slower : 
% stacked_image2=inpaintn(stacked_image2); % see Fex : https://fr.mathworks.com/matlabcentral/fileexchange/27994-inpaint-over-missing-data-in-1-d-2-d-3-d-nd-arrays?s_tid=srchtitle

for ck = 1:m
    figure
    data2D = squeeze(stacked_image2(:,:,ck));
    data2Dmean(ck) = mean(data2D,'all','omitnan');
    imagesc(data2D);
    caxis([mi ma]);
    colorbar('vert');
    
end

% mean 2D interpolated data plot

figure
plot(x_axis,data2Dmean);

[gauri]

Thanks a lot for your kind cooperation. It worked well. I will keep the dimensions of stacked_image unchanged.

I appreciate your kind help.

[Mathieu NOE]

as always, my pleasure !

[gauri]

Thank you very much for your help. Further, once again I need your help to retrieve the data using shape file as follows;

mask = reshape(isinterior(polygon, lon(:), lat(:)), [length(lat), length(lon)]);

% Convert the logical mask to uint16 format

maskConverted = uint16(mask);

% Apply the converted mask to the raster

imgMasked = img .* maskConverted;

The size of imgMasked are same as that of img(5490 5490 7). However, it contains non zero data over shape file area and rest of the data are zeros only. I wanted to retrieve non zeros data out of imgMasked with the size of (1153 1315 7) as size of the shape file. I request you to please suggest me how to retrieve the desired data.

I would be highly obliged to you.

[Mathieu NOE]

hello again

can you share imgMasked as a mat file ?

[gauri]

Thank you very much for your generous support. This is the link for imgMasked.mat file;

https://drive.google.com/file/d/1COeEtDXcX3FW7HXXGzYfLBR0cqrKjRhE/view?usp=drive_link

[Mathieu NOE]

hello

I believe the new size should be (1315 x 1153 x 7) and not (1153 x 1315 x 7) if you want to keep only the non zero elements

the code to crop the image is simple :

imgMasked2 = imgMasked(1:1315,1:1153,:);

[gauri]

Thanks for your kind suggestions. Actually it worked but it is giving some zero values also in last few lines. Therefore, dimensions may be little bit different then 1315 x 1153 x 7. I request you to please suggest me how to figure out the actual dimensions of imgMasked non zero values.

I appreciate your kind help and cooperation.

[Mathieu NOE]

try this instead :

load('imgMasked.mat')
[m,n,p] = size(imgMasked);

m1 = find(mean(imgMasked(:,:,1),1)<eps,1,'first') - 1;
m2 = find(mean(imgMasked(:,:,1),2)<eps,1,'first') - 1;

imgMasked2 = imgMasked(1:m2,1:m1,:); 

[gauri]

Thank you very much for your help. I am using it as follows

% Generate a mask by testing if the raster's geographic grid points are inside the polygon

mask = reshape(isinterior(polygon, lon(:), lat(:)), [length(lat), length(lon)]);

% Convert the logical mask to uint16 format

maskConverted = uint16(mask);

% Apply the converted mask to the raster

imgMasked = stacked_layer .* maskConverted;

[m, n, p] = size(imgMasked);

m1 = find(mean(imgMasked(:,:,1),1)<eps,1,'first') - 1;

m2 = find(mean(imgMasked(:,:,1),2)<eps,1,'first') - 1;

imgMasked2 = imgMasked(1:m2,1:m1,:);

disp(size(imgMasked2));

it is showing 0x0x7 dimensions. I request you to please suggest me how to fix it.

[Mathieu NOE]

hello again

what are the values of m, n, p , m1 and m2 ?

[gauri]

m = 5490 n = 5490 p = 7 m1 = 0 m2 = 0 These are the values I get.

[Mathieu NOE]

funny

is it the exact same imgMasked data that I got from you or are you using another data set ?

if not, can you share again the data ?

[gauri]

It is a different file but generated with same procedure. The link for this file is as follows;

https://drive.google.com/file/d/1JhJF-NHpzZiXLXcN6aKpQvahwqHsh8tv/view?usp=sharing

Thank you so much for your help.

[Mathieu NOE]

ok I understand now why I would not work , here a more robust approach - tested on both data sets

m1 = (mean(imgMasked(:,:,1),1)>eps);
m2 = (mean(imgMasked(:,:,1),2)>eps);

imgMasked2 = imgMasked(m2,m1,:); 

[gauri]

Thank you very much for your kind help. It worked very well for all files.

[Mathieu NOE]

good news !

[Aksh]

I have requested matlab community on my out of memory issue in my matlab code. I request you to please have a look on my question on the following link;

https://in.mathworks.com/matlabcentral/answers/2128211-how-to-resolve-the-out-of-memory-issue-in-shape-file-usage

I would be highly obliged to you for your kind cooperation and help.