Se añaden mejoras e indices al cálculo, se programan imagenes de a06

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'];

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

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

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_I<m_I
+%         V=Indice*Indice';
+%     else
+%         V=Indice'*Indice;
+%     end
+%     [L,U]=lu(V);
+%     Det_L(i)=det(L);
+%     Det_U(i)=det(U);
+%     if Det_L(i)!=0 || Det_L(i)!=0
+    if Uno==1
+        IV(:,:,j)=Indice;
+        Nombre_IV{j,1}=lista{i};
+        j=j+1;
+    end
+end
+%con esto, tenemos el set de indices que realmente se calcularon y se
+%obtienen sus nombres para poder graficarlos con un subplot
+%recurrimos nuevamente a obtener la cantidad de los nombre de indices
+%recuperados
+IndiceVegetacion=[];
+[n m]=size(Nombre_IV)
+%definir tamaño de subplot
+fa=factor(n);
+N=length(fa);
+n_a=ceil(N/2);
+%multiplicar de a  a n_a
+g1=prod(fa(1:n_a));
+%multiplicar de n_a+1 a NN
+g2=prod(fa(n_a+1:N));
+%creamos figure
+figure
+scale=0.3;
+for i=1:n
+I(:,:,i)=mat2gray(IV(:,:,i),[min(min(IV(:,:,i))) max(max(IV(:,:,i)))]);    
+B(:,:,i) = imresize(I(:,:,i), scale);
+subplot(g1,g2,i);
+imshow(B(:,:,i));
+title(Nombre_IV{i,1});
+end
+%Guardar, se hace separado para asegurar que se grafiquen las imagenes en
+%primer lugar
+%hora se procede a guardar cada set de imagenes como geotiff
+for i=1:n
+ [I_16b,Map_16b]=gray2ind(I(:,:,i),65536);
+  savefile_TIF=strcat(Carpeta_output,'/',Nombre_Proceso,'/',COD_SAT_IMG,'_',Nombre_IV{i});
+  geotiffwrite(savefile_TIF,I_16b,INFO.nir,'GeoKeyDirectoryTag',UTM.geokey); 
+end
+
+
+% I=mat2gray(Indice,[min(min(Indice)) max(max(Indice))]);
+% [I_16b,Map_16b]=gray2ind(I,65536);
+
+% 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);

grafica_indices2.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='LT52330762000323CUB00'
+%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_I<m_I
+%         V=Indice*Indice';
+%     else
+%         V=Indice'*Indice;
+%     end
+%     [L,U]=lu(V);
+%     Det_L(i)=det(L);
+%     Det_U(i)=det(U);
+%     if Det_L(i)!=0 || Det_L(i)!=0
+    if Uno==1
+        IV(:,:,j)=Indice;
+        Nombre_IV{j,1}=lista{i};
+        j=j+1;
+    end
+end
+%con esto, tenemos el set de indices que realmente se calcularon y se
+%obtienen sus nombres para poder graficarlos con un subplot
+%recurrimos nuevamente a obtener la cantidad de los nombre de indices
+%recuperados
+IndiceVegetacion=[];
+[n m]=size(Nombre_IV)
+%definir tamaño de subplot
+fa=factor(n);
+N=length(fa);
+n_a=ceil(N/2);
+%multiplicar de a  a n_a
+g1=prod(fa(1:n_a));
+%multiplicar de n_a+1 a NN
+g2=prod(fa(n_a+1:N));
+%creamos figure
+figure
+scale=0.3;
+for i=1:n
+I(:,:,i)=mat2gray(IV(:,:,i),[min(min(IV(:,:,i))) max(max(IV(:,:,i)))]);    
+B(:,:,i) = imresize(I(:,:,i), scale);
+subplot(g1,g2,i);
+imshow(B(:,:,i));
+title(Nombre_IV{i,1});
+end
+%Guardar, se hace separado para asegurar que se grafiquen las imagenes en
+%primer lugar
+%hora se procede a guardar cada set de imagenes como geotiff
+for i=1:n
+ [I_16b,Map_16b]=gray2ind(I(:,:,i),65536);
+  savefile_TIF=strcat(Carpeta_output,'/',Nombre_Proceso,'/',COD_SAT_IMG,'_',Nombre_IV{i});
+  geotiffwrite(savefile_TIF,I_16b,INFO.nir,'GeoKeyDirectoryTag',UTM.geokey); 
+end
+
+
+% I=mat2gray(Indice,[min(min(Indice)) max(max(Indice))]);
+% [I_16b,Map_16b]=gray2ind(I,65536);
+
+% 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);