|
Hi
I worked with NDVI time series with MODIS sensor, My recommendation would be to first make a Layer Stacking (ENVI) of 24 images, sort by date, if you want to create a statistic for your images...with a classified image... You can modify this script ....
close,/all
year=['2000','2001','2002','2003','2004','2005','2006','2007','2008']
month=['01','02','03','04','05','06','07','08','09','10','11','12']
;define dimension
dim=[197,387]
;--------------------------------
;Crea los vectores de salida
;imprime los resultados de la media en un archivo .txt
filetxt1= '/tesis/ndvi_avg/ndvista/sta/mean.txt'
filetxt2= '/tesis/ndvi_avg/ndvista/sta/median.txt'
filetxt3= '/tesis/ndvi_avg/ndvista/sta/min.txt'
filetxt4= '/tesis/ndvi_avg/ndvista/sta/max.txt'
filetxt5= '/tesis/ndvi_avg/ndvista/sta/area.txt'
;defino formato
fmt='(20(f16.3))'
openw,3,filetxt1
openw,4,filetxt2
openw,5,filetxt3
openw,6,filetxt4
openw,7,filetxt5
printf,3,'Year','Month','Mean.Cl1','Mean.Cl.2','Mean.Cl3','Mean.Cl4','Mean.Cl5','Mean.Cl6','Mean.Cl7','Mean.Cl8','Mean.Cl9'
printf,4,'Year','Month','Median.Cl1','Median.Cl.2','Median.Cl3','Median.Cl4','Median.Cl5','Median.Cl6','Median.Cl7','Median.Cl8','Median.Cl9'
printf,5,'Year','Month','Min.Cl1','Min.Cl.2','Min.Cl3','Min.Cl4','Min.Cl5','Min.Cl6','Min.Cl7','Min.Cl8','Min.Cl9'
printf,6,'Year','Month','Max.Cl1','Max.Cl.2','Max.Cl3','Max.Cl4','Max.Cl5','Max.Cl6','Max.Cl7','Max.Cl8','Max.Cl9'
printf,7,'Year','Month','Area.Cl1','Area.Cl.2','Area.Cl3','Area.Cl4','Area.Cl5','Area.Cl6','Area.Cl7','Area.Cl8','Area.Cl9'
;Define matriz avg y clasif
;--------------------------------
;Abre imagene de clasificacion
;--------------------------------
for i = 0,8 do begin
file='/tesis/ndvi_avg/nd/ndn/'+year[i]+'_clasif' ;Se debe cambiar el path a lo que corresponda
clsf=make_array(dim,/BYTE) ; columna,fila, tipo
openr,1,file
readu,1,clsf
close, 1
;abro im?genes promedio
for mes = 0,11 do begin
filavg ='/tesis//ndvi_avg/ndvista/MOD13A3.B0.'+year[i]+month[mes] ;Se debe cambiar el path a lo que corresponda
im_avg=make_array(dim,/float)
close, 1
openr,1,filavg
readu,1,im_avg
;--------------------------------
;--------------------------------
;calcula media
;nombre de los pixeles que corresponden a los climas
cl1=where(clsf eq 1,n1)
cl2=where(clsf eq 2,n2)
cl3=where(clsf eq 3,n3)
cl4=where(clsf eq 4,n4)
cl5=where(clsf eq 5,n5)
cl6=where(clsf eq 6,n6)
cl7=where(clsf eq 7,n7)
cl8=where(clsf eq 8,n8)
cl9=where(clsf eq 9,n9)
;calcula el area por zona
area=make_array(9,/float)
area[0]=n1*0.926^2
area[1]=n2*0.926^2
area[2]=n3*0.926^2
area[3]=n4*0.926^2
area[4]=n5*0.926^2
area[5]=n6*0.926^2
area[6]=n7*0.926^2
area[7]=n8*0.926^2
area[8]=n9*0.926^2
;Extrae valores por clase
;Define matriz resultado
results_mean=make_array(9,/float)
results_median=make_array(9,/float)
results_max=make_array(9,/float)
results_min=make_array(9,/float)
results_zero=make_array(9,/float)
clima=cl1
;chapito=mean(im_avg[clima])
for kg=0,8 do begin
case kg of
0: clima=cl1
1: clima=cl2
2: clima=cl3
3: clima=cl4
4: clima=cl5
5: clima=cl6
6: clima=cl7
7: clima=cl8
8: clima=cl9
endcase
if (clima[0] eq -1) then continue
results_mean[kg]=MEAN(im_avg[clima])
results_median[kg]=median(im_avg[clima])
results_max[kg]=max(im_avg[clima])
results_min[kg]=min(im_avg[clima])
results_zero[kg]=(histogram(im_avg[clima],min=0))[0]
wait, 0.0001
endfor
printf,3,year[i],month[mes],results_mean[0],results_mean[1],results_mean[2],results_mean[3],results_mean[4],results_mean[5],results_mean[6],results_mean[7],results_mean[8],format=fmt
printf,4,year[i],month[mes],results_median[0],results_median[1],results_median[2],results_median[3],results_median[4],results_median[5],results_median[6],results_median[7],results_median[8],format=fmt
printf,5,year[i],month[mes],results_min[0],results_min[1],results_min[2],results_min[3],results_min[4],results_min[5],results_min[6],results_min[7],results_min[8],format=fmt
printf,6,year[i],month[mes],results_max[0],results_max[1],results_max[2],results_max[3],results_max[4],results_max[5],results_max[6],results_max[7],results_max[8],format=fmt
printf,7,year[i],month[mes],area[0],area[1],area[2],area[3],area[4],area[5],area[6],area[7],area[8],format=fmt
close, 1
endfor
endfor
close,/all
end
|