How do I interpolate and smoothen ndvi time series?

Question

Devendra 2024-4-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2109056-how-do-i-interpolate-and-smoothen-ndvi-time-series

评论： Aksh 2024-6-13

sample_ndvi.csv

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

2 个评论
显示无隐藏无

John D'Errico 2024-4-19

编辑：John D'Errico 2024-4-19

在 MATLAB Online 中打开

sample_ndvi.csv

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 2024-4-19

编辑：Devendra 2024-4-19

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 2024-4-22

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2109056-how-do-i-interpolate-and-smoothen-ndvi-time-series#answer_1445766

在 MATLAB Online 中打开

sample_ndvi.csv

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'

39 个评论
显示 37更早的评论隐藏 37更早的评论

Mathieu NOE 2024-4-22

编辑：Mathieu NOE 2024-4-22

在 MATLAB Online 中打开

Test_NaN_fill.csv

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>

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>

figure(1)

surf(zz)

xlabel 'Time'

ylabel 'Series'

% "view" from the top

figure(2)

imagesc(zz)

Mathieu NOE 2024-4-22

在 MATLAB Online 中打开

Test_NaN_fill.csv

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>

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>

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>

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 2024-4-25

编辑：Devendra 2024-4-25

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 2024-4-25

在 MATLAB Online 中打开

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 2024-4-25

编辑：Devendra 2024-4-25

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 2024-4-25

在 MATLAB Online 中打开

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 2024-4-28

编辑：gauri 2024-4-28

在 MATLAB Online 中打开

test_phenology.m

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 2024-5-6

ok let's do that !

gauri 2024-5-13

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 2024-5-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2109056-how-do-i-interpolate-and-smoothen-ndvi-time-series#answer_1453551

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

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

35 个评论
显示 33更早的评论隐藏 33更早的评论

Mathieu NOE 2024-5-14

在 MATLAB Online 中打开

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

Mathieu NOE 2024-5-15

在 MATLAB Online 中打开

env_secant.m

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

Mathieu NOE 2024-5-15

在 MATLAB Online 中打开

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

Mathieu NOE 2024-5-22

good news !

Aksh 2024-6-13

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.

请先登录，再进行评论。

How do I interpolate and smoothen ndvi time series?

2 个评论
显示无隐藏无

采纳的回答

39 个评论
显示 37更早的评论隐藏 37更早的评论

更多回答（1 个）

35 个评论
显示 33更早的评论隐藏 33更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

How do I interpolate and smoothen ndvi time series?

2 个评论 显示 无隐藏 无

采纳的回答

39 个评论 显示 37更早的评论隐藏 37更早的评论

更多回答（1 个）

35 个评论 显示 33更早的评论隐藏 33更早的评论

2 个评论
显示无隐藏无

39 个评论
显示 37更早的评论隐藏 37更早的评论

35 个评论
显示 33更早的评论隐藏 33更早的评论