Se añaden algunas mejoas en procesamiento por paquetes de imagenes

corte_imagen.m

@@ -2,10 +2,12 @@
 function [IMc, Rc, X ,R , INFO]=corte_imagen(BaseDir,File,UTM_x,UTM_y)
 %cortar seccion definida
 
-
 [X, R] = geotiffread(File);
 cd(BaseDir);
+pwd;
+if ~isdir('Cortes_Imagenes')
 mkdir('Cortes_Imagenes');
+end
 cd('Cortes_Imagenes');
 ruta_actual=pwd;
 INFO = geotiffinfo(File);
@@ -58,8 +60,8 @@
 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))];
+Rc.YLimWorld=[Yr(min(upper,lower)) Yr(max(upper,lower))]
+Rc.XLimWorld=[Xr(min(left,right)) Xr(max(left,right))]
 
 % 
 % [left_long, upper_lat] = projfwd(INFO, X_1(2), X_1(1))

desmembrar.m

@@ -5,9 +5,9 @@
   str=Signo;
   k_str=strfind(patron,str);
   n= length(k_str);
+
   for i=1:n-1
      nombre{i,1}=Ruta(k_str(i)+1:k_str(i+1)-1); 
   end
   nombre{n,1}=Ruta(k_str(n)+1:length(Ruta));
-
-end
+  end

geotiff2shape.m

@@ -0,0 +1,15 @@
+tic
+%comprimir a solo valores <0 en indice NDVI
+ [I,J]=find(v_iv<0);%usar J que entrega todos los índices en que ocurre
+ r=1;
+ for i=1:length(v_iv)
+     if isempty(find(J==i))
+        short_v_iv(r)=v_iv(i);
+        short_v_lon(r)=b_lon(i);
+        short_v_lat(r)=b_lat(i);
+        r=1+r;
+     end
+ end
+     test2 = struct('Geometry', 'Multipoint', 'id', num2cell(v_iv), 'X', num2cell(b_lon), 'Y', num2cell(b_lat));
+     %shapewrite(test2, 'test2'); 
+toc

get_list_files.m

@@ -8,3 +8,4 @@
 list_dir=dir(fullfile(path,type));
 list_dir={list_dir.name};
 out=list_dir;
+end

img_lands_one.m

@@ -0,0 +1,104 @@
+Ruta='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/img_ANG006_IMG/IMG/'
+cd(Ruta);
+pwd
+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
+Rc=[];
+while ~feof(FID)    
+leer_linea=fgetl(FID); 
+File='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/img_ANG006_IMG/IMG';
+TOTDIR='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/img_ANG006_IMG/TOT/';
+Ruta='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/img_ANG006_IMG/IMG';
+Nombre_Proceso='ANG006_IMG';
+dese=desmembrar(leer_linea,'/');
+[dn,dm]=size(dese);
+name_des=char(cellstr(dese(dn-1:1)))
+set_DIR=['/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/img_ANG006_IMG/IMG/',name_des]
+cd(set_DIR)
+set_DIR
+pwd
+%'listar tif en carpeta
+[A,L]=unix('find -name "*.TIF"');
+des=desmembrar(L,'./');
+[nd,md]=size(des);
+%usar L{i}(2:end) para extraer nombre de archivo
+%nd=1
+I=zeros(711,498,nd);
+for i=1:nd
+
+%'leer cada imagen tif
+%'cortar cada imagen tif
+name=char(cellstr(des{i}(2:end)));
+FileName=[set_DIR '/' name];
+if ~strcmp(UTM_x,'') || ~strcmp(CorteY,'')    
+%filtrar y omitir banda 6 que es mas densa
+Rx=Rc
+[IMc,Rc, X, R, INFOx]=corte_imagen(BaseDir,FileName,UTM_x,UTM_y);
+[n,m]=size(IMc);
+end
+bla=char(desmembrar(name,'_'));
+banda=bla(1:end-4);
+switch banda
+    case 'B1'
+        capa=1;
+    case 'B2'
+                capa=2;
+    case 'B3'
+                capa=3;
+    case 'B4'
+                capa=4;
+    case 'B5'
+                capa=5;
+    case 'B6'
+                capa=6;
+    case 'B7'
+                capa=7;
+    case 'B8'
+                capa=8;
+    case 'B9'
+                capa=9;
+    case 'B10'
+                capa=10;
+    case 'B11'
+                capa=11;
+    case 'BQA'
+                capa=12;
+end
+
+
+[n,m]=size(IMc)
+ele_banda=banda
+comparacion=~strcmp(char(banda),'B8')
+%si n!=711 && m!=498
+%cambiar el condicional, ya que no se sabe cuales son las capas que
+%corresponden,esto en caso de que salga error de calclula en esta ocasion
+if ~strcmp(char(banda),'B8')    
+I(:,:,capa)=IMc;
+else    
+%'guardar nueva imagen cortada
+Rc=Rx;
+end
+
+end
+
+
+%guardar imagen en TIFF con geodatos
+filename = [TOTDIR name];
+%imwrite(I,filename,'tif')
+%K = mat2gray(IMc);
+INFOx.GeoTIFFTags.GeoKeyDirectoryTag.GTRasterTypeGeoKey=1
+geotiffwrite(filename, I, Rc,'GeoKeyDirectoryTag', INFOx.GeoTIFFTags.GeoKeyDirectoryTag);
+
+end
+
+toc

indices2txt_ascii.m

@@ -0,0 +1,47 @@
+clear
+%leer directorio en que se encuentran las imágenes:
+BaseDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat';
+MatDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/Output/ANG006_IMG/Mat_Matlab';
+XYZdata='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/Output/ANG006_IMG'
+cd(MatDir);
+Proyecto='ANG006';
+[A lista]=unix(['find -name "*.mat"']);
+RutaLog=[BaseDir,'/Logs'];
+archivos_MAT=['archivos_MAT_' Proyecto '.txt'];
+%se crea archivo txt con listado de .mat del grupo
+FID=fopen([RutaLog,'/',archivos_MAT],'w+');
+fprintf(FID,lista);
+fclose(FID);
+%se lee el archivo anterior
+charMAT='_.mat';
+fid=fopen([RutaLog,'/',archivos_MAT],'r');
+MAT=0
+MATset={};
+while ~feof(fid)
+    leer_linea=fgetl(fid);
+    cola=leer_linea(end-3:end);
+    ArchivoMat=leer_linea(3:end);
+    %cargar archivo .mat
+    PathMat=[MatDir '/' ArchivoMat];
+    load(PathMat);
+    names = fieldnames(IndiceVegetacion);
+    [n,m]=size(names);
+    %ver que no tenga ceros all(indice(:))==1
+    for i=1:n
+        indicename= char( names( i ) );
+        indice=getfield(IndiceVegetacion, indicename);
+        if all(indice(:))==1
+       % FileName= INFO.nir.Filename
+        NameFile=[ArchivoMat(1:end-4) '_' indicename]; 
+        p=R.nir.DeltaX;
+        XLimWorld=R.nir.XLimWorld;
+        YLimWorld=R.nir.YLimWorld;
+        [v_iv,v_lat,v_lon]=map2vector(indice,p,XLimWorld,YLimWorld);
+        GEO_data=map_vector2ascii(v_iv,v_lat,v_lon,XYZdata,NameFile);
+        end
+    end
+
+
+    %se descomprime aqui el archivo (dentro de carpeta)
+end
+fclose(fid);

indices_vegetacion.m

@@ -127,7 +127,7 @@
 
 if var_mat==1
 savefile=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'/',COD_SAT_IMG,'.mat');
-save(savefile,'R_l','L_l','DN7','UTM','IndiceVegetacion');
+save(savefile,'R_l','L_l','DN7','UTM','IndiceVegetacion','R','INFO');
 end
 
 

map_vector2ascii.m

@@ -0,0 +1,8 @@
+function GEO_data=map_vector2ascii(v_iv,v_lat,v_lon,Dir,NameFile)
+v_iv=v_iv';
+v_lat=v_lat';
+v_lon=v_lon';
+ID=[1:length(v_iv)]';
+GEO_data=[ID,v_iv,v_lat,v_lon];
+dlmwrite([Dir '/XYZDAT/' NameFile '.dat'],GEO_data);
+end

test_indices_ANG006.m

@@ -31,7 +31,7 @@
 COD_SAT_IMG=char(nombres(length(nombres)-1));
 BitDepth=16;
 %guarda imagen NDVI
-Im=1
+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]); 
@@ -42,7 +42,7 @@
 IM=[];
 Ib=[];
  %guarda imagen SR
- Im=2 
+ 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]); 
@@ -54,7 +54,7 @@
   IM=[];
 Ib=[];
   %guarda imagen NDWI
-  Im=3
+  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]); 
@@ -65,7 +65,7 @@
   IM=[];
 Ib=[];
   %guarda imagen SAVI 0.3
-  Im=4
+  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]); 
@@ -76,7 +76,7 @@
   IM=[];
 Ib=[];
   %guarda imagen SAVI 0.5
-  Im=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]); 
@@ -87,7 +87,7 @@
     IM=[];
 Ib=[];
   %guarda imagen MSI
-  Im=6
+  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]); 
@@ -98,7 +98,7 @@
     IM=[];
 Ib=[];
   %guarda imagen II
-  Im=7
+  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]); 

vector2shape.m

@@ -0,0 +1,18 @@
+
+function S_shape=vector2shape(v_iv,v_lon,v_lat,Namefile,Corte)
+
+%comprimir a solo valores <0 en indice NDVI
+ [I,J]=find(v_iv<Corte);%usar J que entrega todos los índices en que ocurre
+ r=1;
+ for i=1:length(v_iv)
+     if isempty(find(J==i))
+        short_v_iv(r)=v_iv(i);
+        short_v_lon(r)=v_lon(i);
+        short_v_lat(r)=v_lat(i);
+        r=1+r;
+     end
+ end
+     S_shape = struct('Geometry', 'Multipoint', 'id', num2cell(v_iv), 'X', num2cell(v_lon), 'Y', num2cell(v_lat));
+     shapewrite(S_shape, Namefile); 
+
+end