Se añade proceso para indices de vegetacion INDICES={'NDVI','SR','DVI…

…','SAVI','TNDI','NDWI','EVI','TGDVI'}

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='img_proyecto';
-ProjectDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/img_proyecto';
+Directorio='DesAtacama';
+ProjectDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/233_076';
 cd(ProjectDir);
 [A lista]=unix(['find -name *.txt']);
 RutaLog=[BaseDir,'/Logs'];
@@ -27,24 +27,17 @@
     %se descomprime aqui el archivo (dentro de carpeta)
 end
 fclose(fid);
-archivos_MTL='archivos_MTL.txt';
+archivos_MTL=['archivos_MTL_' Directorio '.txt'];
 [a b]=size(MTLset);
 FID=fopen([RutaLog,'/',archivos_MTL],'w+');
 for i=1:a
 fprintf(FID,MTLset{i,1});
 end
 fclose(FID);
-
 FID=fopen([RutaLog,'/',archivos_MTL],'r');
-
 savemat=0;
 Proceso='proyecto'
-
-
-  str='/';
-
-
-
+  str='/';  
 while ~feof(fid)
 MTL_lin=fgetl(fid);
 patron=MTL_lin;

buscar_tar.m

@@ -1,8 +1,10 @@
 BaseDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat';
-Directorio='img_ZonaControl';
-DirectorioFin='img_proyecto';
-Directorio_Padre=[BaseDir,'/',Directorio];
-[A lista]=unix(['find ',Directorio_Padre,' -name *.tar.gz']);
+Directorio='DesAtacama';
+DirectorioFin='DesAtacama/Imagenes';
+Dir_tar='Imagenes/Comp';
+Directorio_Padre=[BaseDir,'/',Directorio,'/',Dir_tar];
+cd(Directorio_Padre);
+[A lista]=unix(['find -name *.tar.gz']);
 %guardar lista en archivos
 archivos_targz='archivos_targz.txt';
 RutaLog=[BaseDir,'/Logs'];
@@ -15,7 +17,6 @@
 unix(['mkdir ',BaseDir,'/',NewDir]);
 fid=fopen([RutaLog,'/',archivos_targz],'r');
 while ~feof(fid)
-    cd(BaseDir);
     leer_linea=fgetl(FID); 
     indice=1;
     nombre=desmembrar(leer_linea,'/');
@@ -25,7 +26,7 @@
             indice=j;
         end
     end
-    cd(NewDir)
+    cd([BaseDir,'/',DirectorioFin]);
     for j=indice+1:a-1        
         unix(['mkdir ',nombre{j,1}]);
         cd(nombre{j,1});
@@ -34,8 +35,9 @@
         unix(['mkdir ',nombre_toma(1:length(nombre_toma)-7)]);
         cd(nombre_toma(1:length(nombre_toma)-7));
     %se han creado las carpetas
+    line_tar=[BaseDir,'/',Directorio,'/',Dir_tar,leer_linea(2:length(leer_linea))];
     %se descomprime aqui el archivo (dentro de carpeta)
-    unix(['tar -xvf ',leer_linea]); 
-
+    unix(['tar -xf ',line_tar]); 
+    unix(['rm ',line_tar]);
 end
 fclose(fid);

corte_imagen.m

@@ -1,5 +1,5 @@
 
-function [IMc X R INFO]=corte_imagen(BaseDir,Ruta,File,nombre_corte,UTM_x,UTM_y,largo_pixel)
+function [IMc, X ,R , INFO]=corte_imagen(BaseDir,File,UTM_x,UTM_y)
 [X, R] = geotiffread(File);
 %cortar seccion definida
 cd(BaseDir);
@@ -13,6 +13,7 @@
 %pasa coordenadas a deg utm 19 N a deg
 [x,y] = R.intrinsicToWorld(R.XLimIntrinsic,R.YLimIntrinsic)
 %lognitud pixel largo-ancho [m]
+largo_pixel= INFO.GeoTIFFTags.ModelPixelScaleTag(1);
 p=largo_pixel; %en metros
 %entregar los dos puntos limite del corte
 %X_a__dx_1___
@@ -35,10 +36,10 @@
 X_1=UTM_x;
 X_2=UTM_y;
 
-left=abs(x(1)-X_1(1))/p
-right=abs(x(1)-X_1(2))/p
-upper=abs(y(1)-X_2(1))/p
-lower=abs(y(1)-X_2(2))/p
+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))

cuadros_geograficos_zc.m

@@ -0,0 +1,74 @@
+cuadro_zc.nombre={
+    'Río Loa Alto';
+    'Salar de Ascotan';
+    'Geyser del Tatio';
+    'Oasis de Atacama';
+    'Volcán Licancabur';
+    'Ayllus de San Pedro de Atacama';
+    'Salar de Atacama';
+    'Salar de Tara';
+    'Salar de aguas Calientes';
+    'Sistema de Lagunas Miscanti y Miñiques';
+    'Salar de Punta Negra';
+    'Salar de Pujsa';
+    'Laguna Lejía';
+    'Sistema Hidrológico Soncor'
+    };
+
+%valores en wgs85
+cuadro_zc.x={
+[-68.85320,-68.09895]
+[-68.26799,-67.85977]
+[-68.08994,-67.81551]
+[-68.30257,-68.15484]
+[-67.78503,-67.60751]
+[-68.65887,-68.02238]
+[-69.06626,-68.88556]
+[-68.97034,-68.86352] 
+[-70.60807,-70.36595] 
+[-70.58620,-70.44180] 
+[-70.63696,-70.48724] 
+[-70.09172,-69.93866] 
+[-69.56745,-69.49855] 
+[-68.40515,-68.18795] 
+[-68.67859,-68.53386] 
+};
+
+
+cuadro_zc.y={
+[-21.08443,-21.88897] 
+[-22.26893,-22.59543] 
+[-22.48025,-22.97415] 
+[-22.83610,-22.98443] 
+[-23.34728,-23.54745]  
+[-22.93933,-23.86970]  
+[-24.34012,-24.75150]
+[-22.41437 ,-22.50059]
+[-24.67875 ,-25.22364]
+[-24.45511 ,-24.59138]
+[-23.00353 ,-23.54435]
+[-21.40138 ,-21.54488]
+[-21.61380 ,-21.73473]
+[-21.41117 ,-21.72944] 
+[-24.88029 ,-25.12285]
+};
+
+[n m]=size(cuadro_zc.x);
+utmzone=19;
+hemi='N';
+for i=1:n
+    Lat=cuadro_zc.y{i}(1);
+    Lon=cuadro_zc.x{i}(1);
+[x,y,utmzone,utmhemi] = wgs2utm(Lat,Lon,utmzone,hemi);
+cuadro_zc.UTMx{i}(1)=x;
+cuadro_zc.UTMy{i}(1)=y;
+    Lat=cuadro_zc.y{i}(2);
+    Lon=cuadro_zc.x{i}(2);
+[x,y,utmzone,utmhemi] = wgs2utm(Lat,Lon,utmzone,hemi);
+cuadro_zc.UTMx{i}(2)=x;
+cuadro_zc.UTMy{i}(2)=y;
+end
+%pasar información de recuadros de zona de control a database
+%diccionario METADATA LANDSAT:
+%http://glovis.usgs.gov/ImgViewer/landsat_dictionary.html
+

indices_landsat.m

@@ -0,0 +1,41 @@
+%buscar cada cuadro en la lista de imagenes de proyecto
+%entregar cuatro ´indices en los que encuentra cada punto.
+%esto indica bajo que imagenes landsat se encuentra el recuadro
+%Puntos_x es un par de valores que indican minimo y maximo horizontal
+%Puntos_y un par de valores que indican minimo maximo vertical
+%Arreglo_x es una lista de valores de dos columnas que contiene los valores
+%minimo y maximo horizontal de las imagenes landsat en UTM 19N
+%Arreglo_y es una lista de valores de dos columnas que contiene los valores
+%minimo y maximo vertical de las imagenes landsat en UTM 19N
+function [x1 x2 y1 y2]=indices_landsat(Puntos_x,Puntos_y,Arreglo_x,Arreglo_y)
+[n m]=size(Arreglo_x);%el mismo que Arreglo_y n filas y m=2 columnas
+%tamaño de Puntos_- es de 1 fila, 2 columnas
+j=1;
+l=1;
+h=1;
+k=1;
+x1=0;
+x2=0;
+y1=0;
+y2=0;
+for i=1:n
+   if and(Arreglo_y(i,1)<=Puntos_y(1),Puntos_y(1)<=Arreglo_y(i,2))
+       y1(h)=i;
+       h=h+1;
+   end
+   if and(Arreglo_y(i,1)<=Puntos_y(2),Puntos_y(2)<=Arreglo_y(i,2))
+       y2(k)=i;
+       k=1+k;
+   end    
+   if and(Arreglo_x(i,1)<=Puntos_x(1),Puntos_x(1)<=Arreglo_x(i,2))
+       x1(j)=i;
+       j=1+j;
+   end
+   if and(Arreglo_x(i,1)<=Puntos_x(2),Puntos_x(2)<=Arreglo_x(i,2))
+       x2(l)=i;
+       l=1+l;
+   end
+
+
+end
+end

indices_vegetacion.m

@@ -0,0 +1,121 @@
+
+%definir directorio de trabajo
+%setear en este directorio
+%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)
+[n m]=size(INDICES);
+cd(BaseDir);
+var_mat=savemat;
+Carpeta_output='Output';%salida, guardar resultados
+i=1;
+Landsat=5;
+Carpeta_img=strcat(BaseDir,'/img_',Nombre_Proceso);
+%Nombres bandas:[blue,green,red,NIR,SWIR,-]
+hemisferio='S';
+%reflactancia NIR 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
+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;
+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;
+            INFO_x=INFO.nir;
+            Indice=NDVI;
+       case 'SR'
+           SR=R_l(:,:,1)./R_l(:,:,2);
+           INFO_x=INFO.nir;
+            Indice=SR;
+       case 'DVI'
+           DVI=R_l(:,:,1)-R_l(:,:,2);
+           INFO_x=INFO.nir;
+           Indice=DVI;
+       case 'SAVI'
+           L=0.5;
+            Numerador=R_l(:,:,1)-R_l(:,:,2);
+            Denominador=R_l(:,:,1)+R_l(:,:,2)+L;
+           SAVI=(Numerador./Denominador)*(1+L); 
+           INFO_x=INFO.nir;
+            Indice=SAVI;
+       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);
+            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;
+            INFO_x=INFO.swir;
+            Indice=NDWI;
+       case 'EVI'
+           G=2.5;
+           L=1;
+           C_1=6;
+           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;
+           INFO_x=INFO.nir;
+            Indice=EVI;
+       case 'TGDVI'            
+            La_1=0.56e-6;%en metros
+            La_2=0.66e-6;
+            La_3=0.83e-6;            
+            Numerador_1=R_l(:,:,1)-R_l(:,:,2);
+            Denominador_1=La_3-La_2;            
+            Numerador_2=R_l(:,:,2)-R_l(:,:,3);
+            Denominador_2=La_2-La_1;            
+            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;
+            INFO_x=INFO.blue;
+            Indice=TGDVI;
+   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');
+end
+
+
+end