From 42f1b2acb00b3c5f29126aef6579d819a70e6bcd Mon Sep 17 00:00:00 2001 From: dpineiden Date: Thu, 27 Mar 2014 10:19:21 -0300 Subject: [PATCH] =?UTF-8?q?Se=20a=C3=B1aden=20mejoras=20e=20indices=20al?= =?UTF-8?q?=20c=C3=A1lculo,=20se=20programan=20imagenes=20de=20a06?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- AnalisisDirectorio.m | 4 +- buscar_tar.m | 8 ++- corte_imagen.m | 47 +++++++++--- grafica_indices.m | 77 ++++++++++++++++++++ grafica_indices2.m | 77 ++++++++++++++++++++ grafica_indices2009.m | 148 ++++++++++++++++++++++++++++++++++++++ grafica_indices2010.m | 77 ++++++++++++++++++++ grafica_indices_ANG006.m | 78 ++++++++++++++++++++ graficaindices2006.m | 77 ++++++++++++++++++++ graficaindices2007.m | 77 ++++++++++++++++++++ graficaindices2008.m | 77 ++++++++++++++++++++ grilla_ndvi.m | 114 ++++++++++++++++++++++++++++++ indices_vegetacion.m | 115 ++++++++++++++++-------------- reflactancia.m | 15 ++-- test_indices.m | 6 +- test_indices2.m | 11 +++ test_indices2010.m | 11 +++ test_indices_ANG006.m | 149 +++++++++++++++++++++++++++++++++++++++ testindices2006.m | 11 +++ testindices2007.m | 11 +++ testindices2008.m | 11 +++ testindices2009.m | 11 +++ 22 files changed, 1137 insertions(+), 75 deletions(-) create mode 100644 grafica_indices.m create mode 100644 grafica_indices2.m create mode 100644 grafica_indices2009.m create mode 100644 grafica_indices2010.m create mode 100644 grafica_indices_ANG006.m create mode 100644 graficaindices2006.m create mode 100644 graficaindices2007.m create mode 100644 graficaindices2008.m create mode 100644 grilla_ndvi.m create mode 100644 test_indices2.m create mode 100644 test_indices2010.m create mode 100644 test_indices_ANG006.m create mode 100644 testindices2006.m create mode 100644 testindices2007.m create mode 100644 testindices2008.m create mode 100644 testindices2009.m diff --git a/AnalisisDirectorio.m b/AnalisisDirectorio.m index 376dc88..29e649e 100644 --- a/AnalisisDirectorio.m +++ b/AnalisisDirectorio.m @@ -1,7 +1,7 @@ %leer directorio en que se encuentran las imágenes: BaseDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat'; -Directorio='DesAtacama'; -ProjectDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/233_076'; +Directorio='ANG006_IMG'; +ProjectDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/ANG006_IMG/IMG'; cd(ProjectDir); [A lista]=unix(['find -name *.txt']); RutaLog=[BaseDir,'/Logs']; diff --git a/buscar_tar.m b/buscar_tar.m index 9040cee..5ccc1d0 100644 --- a/buscar_tar.m +++ b/buscar_tar.m @@ -1,7 +1,9 @@ BaseDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat'; -Directorio='DesAtacama'; -DirectorioFin='DesAtacama/Imagenes'; -Dir_tar='Imagenes/Comp'; +Directorio='ANG006_IMG'; +%directorio en donde se guardan las imag descomprimidas +DirectorioFin='ANG006_IMG/IMG'; +%directorio de imagenes comprimidas +Dir_tar='IMG_COMP'; Directorio_Padre=[BaseDir,'/',Directorio,'/',Dir_tar]; cd(Directorio_Padre); [A lista]=unix(['find -name *.tar.gz']); diff --git a/corte_imagen.m b/corte_imagen.m index c4175ec..cd4eec7 100644 --- a/corte_imagen.m +++ b/corte_imagen.m @@ -1,20 +1,23 @@ -function [IMc, X ,R , INFO]=corte_imagen(BaseDir,File,UTM_x,UTM_y) -[X, R] = geotiffread(File); +function [IMc, Rc, X ,R , INFO]=corte_imagen(BaseDir,File,UTM_x,UTM_y) %cortar seccion definida + + +[X, R] = geotiffread(File); cd(BaseDir); mkdir('Cortes_Imagenes'); -cd('Cortes_Imagenes') +cd('Cortes_Imagenes'); ruta_actual=pwd; INFO = geotiffinfo(File); UTM.geokey=INFO.GeoTIFFTags.GeoKeyDirectoryTag; UTM.tipo=INFO.Projection; -[n,m]=size(X) +[n,m]=size(X); %pasa coordenadas a deg utm 19 N a deg -[x,y] = R.intrinsicToWorld(R.XLimIntrinsic,R.YLimIntrinsic) +[x,y] = R.intrinsicToWorld(R.XLimIntrinsic,R.YLimIntrinsic); %lognitud pixel largo-ancho [m] largo_pixel= INFO.GeoTIFFTags.ModelPixelScaleTag(1); p=largo_pixel; %en metros +py= INFO.GeoTIFFTags.ModelPixelScaleTag(2); %entregar los dos puntos limite del corte %X_a__dx_1___ % dy_1 @@ -32,14 +35,32 @@ % | % ____________dx_2______X_b + + Xr=[R.XLimWorld(1):p:R.XLimWorld(2)]; + Yr=[R.YLimWorld(1):p:R.YLimWorld(2)]; %valores en coordenads de imagen tiff + +UTM_x=round(UTM_x/p)*p; +UTM_y=round(UTM_y/py)*py; + X_1=UTM_x; X_2=UTM_y; +Rc=R; + +%left=round(abs(x(1)-floor(X_1(1)/p)*p)/p) +%right=round(abs(x(1)-ceil(X_1(2)/p)*p)/p) +%upper=round(abs(y(1)-floor(X_2(1)/p)*p)/p) +%lower=round(abs(y(1)-ceil(X_2(2)/p)*p)/p) + + +left=find(Xr==UTM_x(1),1); +right=find(Xr==UTM_x(2),1); +upper=find(Yr==UTM_y(1),1); +lower=find(Yr==UTM_y(2),1); + +Rc.YLimWorld=[Yr(min(upper,lower)) Yr(max(upper,lower))]; +Rc.XLimWorld=[Xr(min(left,right)) Xr(max(left,right))]; -left=round(abs(x(1)-floor(X_1(1)/p)*p)/p) -right=round(abs(x(1)-ceil(X_1(2)/p)*p)/p) -upper=round(abs(y(1)-floor(X_2(1)/p)*p)/p) -lower=round(abs(y(1)-ceil(X_2(2)/p)*p)/p) % % [left_long, upper_lat] = projfwd(INFO, X_1(2), X_1(1)) % [right_long, lower_lat] = projfwd(INFO, X_2(2), X_2(1)) @@ -50,7 +71,13 @@ % [right_d, lower_d] = R.worldToIntrinsic(right_long, lower_lat); % right = round(abs(int64(right_d))/30); % lower = round(abs(int64(lower_d))/30); -IMc = X( min(upper,lower):max(upper,lower), min(left,right):max(left,right) ); +FILAS=[ min(upper,lower) max(upper,lower)]; +COLUMNAS=[min(left,right) max(left,right) ]; +IMc = X( min(upper,lower):max(upper,lower)-1, min(left,right):max(left,right)-1 ); +Rc.RasterSize=size(IMc); +%filas en y +%columnas en x +Rc = maprasterref('RasterSize',Rc.RasterSize,'YLimWorld', Rc.YLimWorld, 'ColumnsStartFrom','north','XLimWorld', Rc.XLimWorld); end \ No newline at end of file diff --git a/grafica_indices.m b/grafica_indices.m new file mode 100644 index 0000000..d64c53b --- /dev/null +++ b/grafica_indices.m @@ -0,0 +1,77 @@ +%pasar los índices de vegetacion a escala de grises, plotear cada uno en un +%subplot y guardar cada imagen +%Vamos a graficar 'IndiceVegetacion' +Nombre_Proceso='Prueba_Indices'; +Carpeta_output='Output'; +COD_SAT_IMG='LM52330761984359AAA03' +%obtengo los nombres de Indice de Vegetacion +lista=fieldnames(IndiceVegetacion); +[n m]=size(lista); +%n entrega la cantidad de elementos o nombres de lista +%revisar cuales son con valor '0' para omitir su gráfica. +IV=[]; +Nombre_IV=''; +j=1; +for i=1:n + Indice=getfield(IndiceVegetacion,lista{i,1}); + [n_I m_I]=size(Indice); + Uno=any(Indice(:)); +% if n_I0 + for cuadrante=1:4 %for3 + %pasa a darle valores para extraer matriz + switch cuadrante + case 1 + caso=1 + step=[-1 -1]; + Np= i - (step(1) + paso_grilla_fila) + Mq= j - (step(2) + paso_grilla_columna) + case 2 + caso=2 + step=[-1 1]; + Np= i - (step(1) + paso_grilla_fila) + %revisa cuando esté terminando la fila o la columna + if final_fila > 0 && i == filas_G(gp) + paso_grilla_fila = final_fila; + end + if final_columna > 0 && j == columnas_G(gp) + paso_grilla_columna = final_columna; + end + + Mq= j + (step(2) - paso_grilla_columna) + case 3 + caso=3 + step=[1 -1]; + Np= i + (step(1) - paso_grilla_fila) + + %revisa cuando esté terminando la fila o la columna + if final_fila > 0 && i == filas_G(gp) + paso_grilla_fila = final_fila; + end + if final_columna > 0 && j == columnas_G(gp) + paso_grilla_columna = final_columna; + end + Mq= j - (step(2) + paso_grilla_columna) + + case 4 + caso=4 + step=[1 1]; + Np= i + (step(1) - paso_grilla_fila) + + %revisa cuando esté terminando la fila o la columna + if final_fila > 0 && i == filas_G(gp) + paso_grilla_fila = final_fila; + end + if final_columna > 0 && j == columnas_G(gp) + paso_grilla_columna = final_columna; + end + + Mq= j + (step(2) - paso_grilla_columna) + end + + %se procede a recortar la matriz mayor + i + j + + step + + lugar_fila=[i:step(1):Np] + lugar_col=[j:step(2):Mq] + A = M(lugar_fila,lugar_col) + %se suman sus valores + Total_A = sum(A(:)); + %se obtiene el área que cubre la matriz. + [na,nb]= size(A); + ns=na*d; + eo=nb*d; + Area_A = ns*eo; + %se obtiene la densidad del valor en cuadrante. + G_ndvi(i+step(1),j+step(2))=Total_A/Area_A + end %for3 + %se retoma el valor de n en paso + paso_grilla_fila = n; + paso_grilla_columna = n; + + else + msgbox('No se puede realizar la operación, revisa los valores de n, d o la matriz M') + end + end %for2 + end%for1 +end \ No newline at end of file diff --git a/indices_vegetacion.m b/indices_vegetacion.m index 6ad78b5..7670189 100644 --- a/indices_vegetacion.m +++ b/indices_vegetacion.m @@ -4,7 +4,7 @@ %Nombre del Proceso %se desean guardar las variables en .mat?savemat 0:no, 1:si %INDICES={"NDVI","SR","DVI","SAVI","TNDI","NDWI","EVI","TGDVI"} -function [NDVI,SR,DVI,SAVI,TNDI,NDWI,EVI,TGDVI]=indices_vegetacion(BaseDir,MTLDIR,savemat,Nombre_Proceso,CorteX,CorteY,INDICES) +function [IndiceVegetacion, INFO,R,UTM, R_l]=indices_vegetacion(BaseDir,MTLDIR,savemat,Nombre_Proceso,CorteX,CorteY,INDICES,Lsavi) [n m]=size(INDICES); cd(BaseDir); var_mat=savemat; @@ -14,63 +14,83 @@ Carpeta_img=strcat(BaseDir,'/img_',Nombre_Proceso); %Nombres bandas:[blue,green,red,NIR,SWIR,-] hemisferio='S'; -%reflactancia NIR 1 +%reflactancia NIR:R 1 Banda=4; -[R_l(:,:,1),L_l(:,:,1),DN7(:,:,1),UTM,COD_SAT_IMG,INFO.nir]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); -%reflactancia ROJA 2 +[R_l(:,:,1),L_l(:,:,1),DN7(:,:,1),UTM,COD_SAT_IMG,INFO.nir,R.nir]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); +%reflactancia ROJA:R 2 Banda=3; -[R_l(:,:,2),L_l(:,:,2),DN7(:,:,2),UTM,COD_SAT_IMG,INFO.red]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); -%reflactancia AZUL 3 -Banda=1; -[R_l(:,:,3),L_l(:,:,3),DN7(:,:,3),UTM,COD_SAT_IMG,INFO.blue]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); -%reflactancia SWIR 4 -Banda=5; -[R_l(:,:,4),L_l(:,:,4),DN7(:,:,4),UTM,COD_SAT_IMG,INFO.swir]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); -NDVI=0; -SR=0; -DVI=0; -SAVI=0; -TNDI=0; -NDWI=0; -EVI=0; -TGDVI=0; +[R_l(:,:,2),L_l(:,:,2),DN7(:,:,2),UTM,COD_SAT_IMG,INFO.red,R.red]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); +R_l(:,:,3)=zeros(size(R_l(:,:,1))); +R_l(:,:,4)=zeros(size(R_l(:,:,1))); + +for i=1:m + switch INDICES{1,i} + case 'TGDVI' + %reflactancia AZUL:R 3 + Banda=1; + [R_l(:,:,3),L_l(:,:,3),DN7(:,:,3),UTM,COD_SAT_IMG,INFO.blue,R.blue]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); + %reflactancia SWIR:R 4 + case 'NDWI' + Banda=5; + [R_l(:,:,4),L_l(:,:,4),DN7(:,:,4),UTM,COD_SAT_IMG,INFO.swir,R.swir]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); + case 'MSI' + Banda=5; + [R_l(:,:,4),L_l(:,:,4),DN7(:,:,4),UTM,COD_SAT_IMG,INFO.swir,R.swir]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); + case 'II' + Banda=5; + [R_l(:,:,4),L_l(:,:,4),DN7(:,:,4),UTM,COD_SAT_IMG,INFO.swir,R.swir]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); + + end +end + +IndiceVegetacion.NDVI=0; +IndiceVegetacion.SR=0; +IndiceVegetacion.DVI=0; +IndiceVegetacion.SAVI=0; +IndiceVegetacion.TNDI=0; +IndiceVegetacion.NDWI=0; +IndiceVegetacion.EVI=0; +IndiceVegetacion.TGDVI=0; + +IndiceVegetacion.SAVI_03=0; +IndiceVegetacion.SAVI_05=0; +IndiceVegetacion.MSI=0; +IndiceVegetacion.II=0; for i=1:m -INDICES{1,i} switch INDICES{1,i} case 'NDVI' Numerador=R_l(:,:,1)-R_l(:,:,2); Denominador=R_l(:,:,1)+R_l(:,:,2); - NDVI=Numerador./Denominador; + IndiceVegetacion.NDVI=Numerador./Denominador; INFO_x=INFO.nir; - Indice=NDVI; case 'SR' - SR=R_l(:,:,1)./R_l(:,:,2); + IndiceVegetacion.SR=R_l(:,:,1)./(R_l(:,:,2)+1); INFO_x=INFO.nir; - Indice=SR; case 'DVI' - DVI=R_l(:,:,1)-R_l(:,:,2); + IndiceVegetacion.DVI=R_l(:,:,1)-R_l(:,:,2); INFO_x=INFO.nir; - Indice=DVI; case 'SAVI' - L=0.5; + L=0.3; Numerador=R_l(:,:,1)-R_l(:,:,2); Denominador=R_l(:,:,1)+R_l(:,:,2)+L; - SAVI=(Numerador./Denominador)*(1+L); + IndiceVegetacion.SAVI_03=(Numerador./Denominador)*(1+L); INFO_x=INFO.nir; - Indice=SAVI; + L=0.5 + Numerador=R_l(:,:,1)-R_l(:,:,2); + Denominador=R_l(:,:,1)+R_l(:,:,2)+L; + IndiceVegetacion.SAVI_05=(Numerador./Denominador)*(1+L); + INFO_x=INFO.nir; case 'TNDI' Numerador=R_l(:,:,1)-R_l(:,:,2); Denominador=R_l(:,:,1)+R_l(:,:,2); pre_TNDI=(Numerador./Denominador) +1; - TNDI=sqrt(pre_TNDI); + IndiceVegetacion.TNDI=sqrt(pre_TNDI); INFO_x=INFO.nir; - Indice=TNDI; case 'NDWI' Numerador=R_l(:,:,1)-R_l(:,:,4); Denominador=R_l(:,:,1)+R_l(:,:,4); - NDWI=Numerador./Denominador; + IndiceVegetacion.NDWI=Numerador./Denominador; INFO_x=INFO.swir; - Indice=NDWI; case 'EVI' G=2.5; L=1; @@ -78,9 +98,8 @@ C_2=7.5; Numerador=R_l(:,:,1)-R_l(:,:,2); Denominador=R_l(:,:,1)+C_1*R_l(:,:,2)-C_2*R_l(:,:,3)+L; - EVI=G*Numerador./Denominador; + IndiceVegetacion.EVI=G*Numerador./Denominador; INFO_x=INFO.nir; - Indice=EVI; case 'TGDVI' La_1=0.56e-6;%en metros La_2=0.66e-6; @@ -92,29 +111,23 @@ pre_TGDVI=Numerador_1./Denominador_1-Numerador_2./Denominador_2; %pasar a 0 valores menores que 0 pre_TGDVI(pre_TGDVI<0)=0; - TGDVI=pre_TGDVI; + IndiceVegetacion.TGDVI=pre_TGDVI; INFO_x=INFO.blue; - Indice=TGDVI; + + case 'MSI' + IndiceVegetacion.MSI=R_l(:,:,4)./R_l(:,:,1); + case 'II' + IndiceVegetacion.II=(R_l(:,:,1)-R_l(:,:,4))./(R_l(:,:,1)+R_l(:,:,4)); + end -% I=mat2gray(Indice,[min(min(Indice)) max(max(Indice))]); -% [I_16b,Map_16b]=gray2ind(I,65536); -% Histograma=imhist(I_16b); -% Bit_depth=16;%'BitDepth'= 16 bits -% NDVI_16=uint16(NDVI*65535); -% savefile_TIF=strcat(Carpeta_output,'/',Nombre_Proceso,'/',COD_SAT_IMG,'_',INDICES{i},'.TIF'); -% geotiffwrite(savefile_TIF,I_16b,INFO_x,'GeoKeyDirectoryTag',UTM.geokey); -% imwrite(I_16b,savefile_TIF,'tif'); -% Info_tif=imfinfo(savefile_TIF); -% savefile_PNG=strcat(Carpeta_output,'/',Nombre_Proceso,'/',COD_SAT_IMG,'_',INDICES{i},'.PNG'); -% imwrite(I_16b,savefile_PNG,'png'); -% Info_png=imfinfo(savefile_PNG); + end if var_mat==1 -savefile=strcat(Carpeta_output,'/',Nombre_Proceso,'/',COD_SAT_IMG,'.mat'); -save(savefile,'R_l','L_l','DN7','UTM'); +savefile=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'/',COD_SAT_IMG,'.mat'); +save(savefile,'R_l','L_l','DN7','UTM','IndiceVegetacion'); end diff --git a/reflactancia.m b/reflactancia.m index 456c194..7d8c5ea 100644 --- a/reflactancia.m +++ b/reflactancia.m @@ -7,7 +7,7 @@ %entregar Path %Entregar hemisferio -function [R_l,L_l,DN7,UTM,COD_SAT_IMG,R]=reflactancia(Banda,Landsat,Path,MTLDIR,Hemisferio,Nombre_Proceso,CorteX,CorteY) +function [R_l,L_l,DN7,UTM,COD_SAT_IMG,INFO,R]=reflactancia(Banda,Landsat,Path,MTLDIR,Hemisferio,Nombre_Proceso,CorteX,CorteY) cd(Path); cd('..'); BaseDir=pwd; @@ -55,7 +55,7 @@ if Landsat==7 %Para DN7 Nombres_Valores={'SUN_ELEVATION','DATE_HOUR_CONTACT_PERIOD'}; -elseif Landsat==5 +elseif Landsat==5 || Landsat == 8 %Para DN5 Nombres_Valores={'SUN_ELEVATION','LANDSAT_SCENE_ID'}; end @@ -99,8 +99,11 @@ FilePath=char(strcat(Carpeta_img,'/',archivo)); INFO = geotiffinfo(FilePath); -if ~strcmp(CorteX,'') | ~strcmp(CorteY,'') -[IMc X R INFO]=corte_imagen(BaseDir,File,CorteX,CorteY); +if ~strcmp(CorteX,'') || ~strcmp(CorteY,'') +[IMc,Rc, X, R, INFOx]=corte_imagen(BaseDir,FilePath,CorteX,CorteY); +X=IMc; +[a v]=size(X); +R=Rc; else [X, R] = geotiffread(FilePath); end @@ -135,8 +138,8 @@ if Landsat==5 %los valores de DN7 DN7(i,j)=(Slope(Banda)*X(i,j))+Intercep(Banda); - else - DN7(i,j)=X(i,j); + elseif Landsat == 7 || Landsat == 8 + DN7(i,j)=X(i,j); end %la radianza es entonces L_l(i,j)=Gain_radiance(Banda)*DN7(i,j)+Bias_radiance(Banda); diff --git a/test_indices.m b/test_indices.m index 07550f4..14b3326 100644 --- a/test_indices.m +++ b/test_indices.m @@ -1,5 +1,5 @@ -INDICES={'NDVI','SR','DVI','SAVI','TNDI','NDWI','EVI','TGDVI'}; -%INDICES={'EVI','TGDVI'}; +INDICES={'NDVI','SR','NDWI'}; +%INDICES={'NDVI','NDWI'}; BaseDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat' File='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2001/11_noviembre/LM52330761984359AAA03_B1.TIF'; Ruta='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2001/11_noviembre'; @@ -7,5 +7,5 @@ MTLDIR='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2001/11_noviembre'; savemat=0; tic -[NDVI,SR,DVI,SAVI,TNDI,NDWI,EVI,TGDVI]=indices_vegetacion(BaseDir,MTLDIR,savemat,Nombre_Proceso,'','',INDICES); +[IndiceVegetacion, INFO,UTM,R_l]=indices_vegetacion(BaseDir,MTLDIR,savemat,Nombre_Proceso,'','',INDICES); toc \ No newline at end of file diff --git a/test_indices2.m b/test_indices2.m new file mode 100644 index 0000000..e7fa39f --- /dev/null +++ b/test_indices2.m @@ -0,0 +1,11 @@ +INDICES={'NDVI','SR','DVI','SAVI','TNDI','NDWI','EVI','TGDVI'}; +%INDICES={'NDVI','NDWI'}; +BaseDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat' +File='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2000/11_noviembre/LT52330762000323CUB00_B1.TIF'; +Ruta='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2000/11_noviembre'; +Nombre_Proceso='Prueba_Indices'; +MTLDIR='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2000/11_noviembre'; +savemat=0; +tic +[IndiceVegetacion, INFO,UTM,R_l]=indices_vegetacion(BaseDir,MTLDIR,savemat,Nombre_Proceso,'','',INDICES); +toc \ No newline at end of file diff --git a/test_indices2010.m b/test_indices2010.m new file mode 100644 index 0000000..7597b90 --- /dev/null +++ b/test_indices2010.m @@ -0,0 +1,11 @@ +INDICES={'NDVI','SR','DVI','SAVI','TNDI','NDWI','EVI','TGDVI'}; +%INDICES={'NDVI','NDWI'}; +BaseDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat' +File='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2010/11_noviembre/LT52330762010318CUB00_B1.TIF'; +Ruta='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2010/11_noviembre'; +Nombre_Proceso='Prueba_Indices'; +MTLDIR='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2010/11_noviembre'; +savemat=0; +tic +[IndiceVegetacion, INFO,UTM,R_l]=indices_vegetacion(BaseDir,MTLDIR,savemat,Nombre_Proceso,'','',INDICES); +toc \ No newline at end of file diff --git a/test_indices_ANG006.m b/test_indices_ANG006.m new file mode 100644 index 0000000..be68737 --- /dev/null +++ b/test_indices_ANG006.m @@ -0,0 +1,149 @@ +INDICES={'NDVI','SR','NDWI','SAVI','MSI','II'}; +Carpeta_output='Output'; +%INDICES={'NDVI','NDWI'}; +%programa que lee línea a línea el LOG MTL +BaseDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat' +LogTXT='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/Logs/archivos_MTL_ANG006_IMG.txt'; +FID=fopen(LogTXT,'r'); +UTM_x=[372560,387500]; +%si entregan datos en UTM sur se deben convertir a norte ya que landsat +%entrega en UTM norte +UTM_y=[6328630,6307300]-10000000; +I=1; +tic +while ~feof(FID) + leer_linea=fgetl(FID); +File='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/img_ANG006_IMG/IMG'; +Ruta='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/img_ANG006_IMG/IMG'; +Nombre_Proceso='ANG006_IMG'; +des=desmembrar(leer_linea,'/'); +[dn,dm]=size(des); +MTLDIR=['/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/img_ANG006_IMG/IMG/' char(des(dn-1,1))]; +savemat=1; +[IndiceVegetacion, INFO,R,UTM,R_l]=indices_vegetacion(BaseDir,MTLDIR,savemat,Nombre_Proceso,UTM_x,UTM_y,INDICES); + + +IM=[]; +Ib=[]; +UTM.geokey.GTRasterTypeGeoKey=1; + +nombres=desmembrar(INFO.nir.Filename,'/'); +COD_SAT_IMG=char(nombres(length(nombres)-1)); +BitDepth=16; +%guarda imagen NDVI +Im=1 + IM=mat2gray(IndiceVegetacion.NDVI,[min(IndiceVegetacion.NDVI(:)) max(IndiceVegetacion.NDVI(:))]); +Lista_NDVI{I,:}={COD_SAT_IMG,',',min(IndiceVegetacion.NDVI(:)),',',max(IndiceVegetacion.NDVI(:))}; + %IM=mat2gray(IndiceVegetacion.NDVI,[-1 1]); + [Ib,Mapb]=gray2ind(IM,2^BitDepth); + savefile_TIF=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'/TIFF/',COD_SAT_IMG,'_',INDICES{1}); + geotiffwrite(savefile_TIF,Ib,R.nir,'GeoKeyDirectoryTag',UTM.geokey); + +IM=[]; +Ib=[]; + %guarda imagen SR + Im=2 + IM=mat2gray(IndiceVegetacion.SR,[min(IndiceVegetacion.SR(:)) max(IndiceVegetacion.SR(:))]); + Lista_SR{I,:}={COD_SAT_IMG,',',min(IndiceVegetacion.SR(:)),',',max(IndiceVegetacion.SR(:))}; + %IM=mat2gray(IndiceVegetacion.SR,[0 2^BitDepth]); + [Ib,Mapb]=gray2ind(IM,2^BitDepth); + savefile_TIF=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'/TIFF/',COD_SAT_IMG,'_',INDICES{2}); + geotiffwrite(savefile_TIF,Ib,R.nir,'GeoKeyDirectoryTag',UTM.geokey); + + + IM=[]; +Ib=[]; + %guarda imagen NDWI + Im=3 + IM=mat2gray(IndiceVegetacion.NDWI,[min(IndiceVegetacion.NDWI(:)) max(IndiceVegetacion.NDWI(:))]); + Lista_NDWI{I,:}={COD_SAT_IMG,',',min(IndiceVegetacion.NDWI(:)),',',max(IndiceVegetacion.NDWI(:))}; + %IM=mat2gray(IndiceVegetacion.NDWI,[-1 1]); + [Ib,Mapb]=gray2ind(IM,2^BitDepth); + savefile_TIF=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'/TIFF/',COD_SAT_IMG,'_',INDICES{3}); + geotiffwrite(savefile_TIF,Ib,R.nir,'GeoKeyDirectoryTag',UTM.geokey); + + IM=[]; +Ib=[]; + %guarda imagen SAVI 0.3 + Im=4 + IM=mat2gray(IndiceVegetacion.SAVI_03,[min(IndiceVegetacion.SAVI_03) max(IndiceVegetacion.SAVI_03)]); + Lista_SAVI_03{I,:}={COD_SAT_IMG,',',min(IndiceVegetacion.SAVI_03),',',max(IndiceVegetacion.SAVI_03)}; + %IM=mat2gray(IndiceVegetacion.NDWI,[-1 1]); + [Ib,Mapb]=gray2ind(IM,2^BitDepth); + savefile_TIF=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'/TIFF/',COD_SAT_IMG,'_',INDICES{4},'_03'); + geotiffwrite(savefile_TIF,Ib,R.nir,'GeoKeyDirectoryTag',UTM.geokey); + + IM=[]; +Ib=[]; + %guarda imagen SAVI 0.5 + Im=5 + IM=mat2gray(IndiceVegetacion.SAVI_05,[min(IndiceVegetacion.SAVI_05) max(IndiceVegetacion.SAVI_05)]); + Lista_SAVI_05{I,:}={COD_SAT_IMG,',',min(IndiceVegetacion.SAVI_05),',',max(IndiceVegetacion.SAVI_05)}; + %IM=mat2gray(IndiceVegetacion.NDWI,[-1 1]); + [Ib,Mapb]=gray2ind(IM,2^BitDepth); + savefile_TIF=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'/TIFF/',COD_SAT_IMG,'_',INDICES{4},'_05'); + geotiffwrite(savefile_TIF,Ib,R.nir,'GeoKeyDirectoryTag',UTM.geokey); + + IM=[]; +Ib=[]; + %guarda imagen MSI + Im=6 + IM=mat2gray(IndiceVegetacion.MSI,[min(IndiceVegetacion.MSI(:)) max(IndiceVegetacion.MSI(:))]); + Lista_MSI{I,:}={COD_SAT_IMG,',',min(IndiceVegetacion.MSI(:)),',',max(IndiceVegetacion.MSI(:))}; + %IM=mat2gray(IndiceVegetacion.NDWI,[-1 1]); + [Ib,Mapb]=gray2ind(IM,2^BitDepth); + savefile_TIF=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'/TIFF/',COD_SAT_IMG,'_',INDICES{5}); + geotiffwrite(savefile_TIF,Ib,R.nir,'GeoKeyDirectoryTag',UTM.geokey); + + IM=[]; +Ib=[]; + %guarda imagen II + Im=7 + IM=mat2gray(IndiceVegetacion.II,[min(IndiceVegetacion.II(:)) max(IndiceVegetacion.II(:))]); + Lista_II{I,:}={COD_SAT_IMG,',',min(IndiceVegetacion.II(:)),',',max(IndiceVegetacion.II(:))}; + %IM=mat2gray(IndiceVegetacion.NDWI,[-1 1]); + [Ib,Mapb]=gray2ind(IM,2^BitDepth); + savefile_TIF=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'/TIFF/',COD_SAT_IMG,'_',INDICES{6}); + geotiffwrite(savefile_TIF,Ib,R.nir,'GeoKeyDirectoryTag',UTM.geokey); + + I=I+1; + + +end + + savefile_ndvi=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'_ndvi'); + savefile_sr=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'_sr'); + savefile_ndwi=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'_ndwi'); + savefile_savi_03=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'_savi_03'); + savefile_savi_05=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'_savi_05'); + savefile_msi=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'_msi'); + savefile_ii=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'_ii'); + +for i=1:I-1 +FID_ndvi=fopen(savefile_ndvi,'a+'); +FID_sr=fopen(savefile_sr,'a+'); +FID_ndwi=fopen(savefile_ndwi,'a+'); +FID_savi_03=fopen(savefile_savi_03,'a+'); +FID_savi_05=fopen(savefile_savi_05,'a+'); +FID_msi=fopen(savefile_msi,'a+'); +FID_ii=fopen(savefile_ii,'a+'); + +fprintf(FID_ndvi,'%s\n',[Lista_NDVI{i,1}{1} Lista_NDVI{i,1}{2} num2str(Lista_NDVI{i,1}{3}) Lista_NDVI{i,1}{4} num2str(Lista_NDVI{i,1}{5})]); +fprintf(FID_sr,'%s\n',[Lista_SR{i,1}{1} Lista_SR{i,1}{2} num2str(Lista_SR{i,1}{3}) Lista_SR{i,1}{4} num2str(Lista_SR{i,1}{5})]); +fprintf(FID_ndwi,'%s\n',[Lista_NDWI{i,1}{1} Lista_NDWI{i,1}{2} num2str(Lista_NDWI{i,1}{3}) Lista_NDWI{i,1}{4} num2str(Lista_NDWI{i,1}{5})]); +fprintf(FID_savi_03,'%s\n',[Lista_SAVI_03{i,1}{1} Lista_SAVI_03{i,1}{2} num2str(Lista_SAVI_03{i,1}{3}) Lista_SAVI_03{i,1}{4} num2str(Lista_SAVI_03{i,1}{5})]); +fprintf(FID_savi_05,'%s\n',[Lista_SAVI_05{i,1}{1} Lista_SAVI_05{i,1}{2} num2str(Lista_SAVI_05{i,1}{3}) Lista_SAVI_05{i,1}{4} num2str(Lista_SAVI_05{i,1}{5})]); +fprintf(FID_msi,'%s\n',[Lista_MSI{i,1}{1} Lista_MSI{i,1}{2} num2str(Lista_MSI{i,1}{3}) Lista_MSI{i,1}{4} num2str(Lista_MSI{i,1}{5})]); +fprintf(FID_ii,'%s\n',[Lista_II{i,1}{1} Lista_II{i,1}{2} num2str(Lista_II{i,1}{3}) Lista_II{i,1}{4} num2str(Lista_II{i,1}{5})]); + +fclose(FID_ndvi); +fclose(FID_sr); +fclose(FID_ndwi); +fclose(FID_savi_03); +fclose(FID_savi_05); +fclose(FID_msi); +fclose(FID_ii); +end + + +toc diff --git a/testindices2006.m b/testindices2006.m new file mode 100644 index 0000000..a671398 --- /dev/null +++ b/testindices2006.m @@ -0,0 +1,11 @@ +INDICES={'NDVI','SR','DVI','SAVI','TNDI','NDWI','EVI','TGDVI'}; +%INDICES={'NDVI','NDWI'}; +BaseDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat' +File='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2007/10_octubre/LT52330762006307COA00_B1.TIF'; +Ruta='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2007/10_octubre'; +Nombre_Proceso='Prueba_Indices'; +MTLDIR='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2007/10_octubre'; +savemat=0; +tic +[IndiceVegetacion, INFO,UTM,R_l]=indices_vegetacion(BaseDir,MTLDIR,savemat,Nombre_Proceso,'','',INDICES); +toc \ No newline at end of file diff --git a/testindices2007.m b/testindices2007.m new file mode 100644 index 0000000..aaf19b8 --- /dev/null +++ b/testindices2007.m @@ -0,0 +1,11 @@ +INDICES={'NDVI','SR','DVI','SAVI','TNDI','NDWI','EVI','TGDVI'}; +%INDICES={'NDVI','NDWI'}; +BaseDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat' +File='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2007/10_octubre/LT52330762007278CUB00_B1.TIF'; +Ruta='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2007/10_octubre'; +Nombre_Proceso='Prueba_Indices'; +MTLDIR='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2007/10_octubre'; +savemat=0; +tic +[IndiceVegetacion, INFO,UTM,R_l]=indices_vegetacion(BaseDir,MTLDIR,savemat,Nombre_Proceso,'','',INDICES); +toc \ No newline at end of file diff --git a/testindices2008.m b/testindices2008.m new file mode 100644 index 0000000..5d96aa9 --- /dev/null +++ b/testindices2008.m @@ -0,0 +1,11 @@ +INDICES={'NDVI','SR','DVI','SAVI','TNDI','NDWI','EVI','TGDVI'}; +%INDICES={'NDVI','NDWI'}; +BaseDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat' +File='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2008/11_noviembre/LT52330762008313CUB00_B1.TIF'; +Ruta='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2008/11_noviembre'; +Nombre_Proceso='Prueba_Indices'; +MTLDIR='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2008/11_noviembre'; +savemat=0; +tic +[IndiceVegetacion, INFO,UTM,R_l]=indices_vegetacion(BaseDir,MTLDIR,savemat,Nombre_Proceso,'','',INDICES); +toc \ No newline at end of file diff --git a/testindices2009.m b/testindices2009.m new file mode 100644 index 0000000..6e70d54 --- /dev/null +++ b/testindices2009.m @@ -0,0 +1,11 @@ +INDICES={'NDVI','SR','DVI','SAVI','TNDI','NDWI','EVI','TGDVI'}; +%INDICES={'NDVI','NDWI'}; +BaseDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat' +File='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2009/11_noviembre/LT52330762009315CUB01_B1.TIF'; +Ruta='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2009/11_noviembre'; +Nombre_Proceso='Prueba_Indices'; +MTLDIR='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2009/11_noviembre'; +savemat=0; +tic +[IndiceVegetacion, INFO,UTM,R_l]=indices_vegetacion(BaseDir,MTLDIR,savemat,Nombre_Proceso,'','',INDICES); +toc \ No newline at end of file