Analisis estadistico indices

analisis_estadistico.m

@@ -0,0 +1,139 @@
+%%definir el directorio y  nombres de archivos que contienen el min-max de
+%%cada imagen en particular
+%%se carga como un objeto: Indice.{lista_nombre,lista_min,lista_max}
+%%directorio
+    directorio='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/Output/';
+    proyecto='ANG006_IMG';
+  directorio_tiff= strcat(directorio,proyecto,'/TIFF/');
+    indices={'ii','msi','ndvi','ndwi','savi_03','savi_05','sr'};
+%%indices a analizar
+    ind_analisis=3;%o [3 5 7] por ejemplo
+    minmax_i=1;
+
+for g=1:length(ind_analisis)
+    Atiff=[];
+    Ltidd=[];
+    Lista_tiff=[];%en cada iteracion se vacia
+%%realizar un analisis para cada indice requerido
+    string_indice=indices{g};    
+    %Leer cada archivo de carpeta tiff indice
+    %obtener la lista de archivos .tif en directorio
+    %find  -iname "*.txt"
+    %buscar solo en este directorio:
+    cd(directorio_tiff);
+    [Atiff, Ltiff]=unix(['find  -maxdepth 1 -iname  "*_',upper(string_indice),'.tif"']); 
+    %otra busqueda, con detalle de cuadrante, en caso de trabajar con varios cuadrantes en proyecto:
+    %Nro_cuadrante=233083 %en anglo
+    %[Atiff2, Ltiff2]=unix(['find  -maxdepth 1 -iname  "*',Nro_cuadrante,'*_NDVI.tif"']);
+    files_tiff= strfind(Ltiff,'.tif');
+    %tomar el nombre como string
+    cantidad_archivos_tiff=length(files_tiff);
+    Lista_tiff{1}=L(3:files_tiff(1)-1);    
+    for i=2:cantidad_archivos_tiff
+        Lista_tiff{i}=L(files_tiff(i-1)+7:files_tiff(i)-1);    
+    end
+    [L_fil_tiff,L_col_tiff]=size(Lista_tiff); %entrega largo de lista, tomando L_col
+%ya no es necesario, se selecciona con antelacion < %recortar el valor de
+%indice>
+    %buscar primero '_' de cada indice
+
+%     for j=1:L_col_tiff
+%         posicion_indice = strfind(Lista_tiff{i},'_');
+%         indice_actual =  Lista_tiff{i}(posicion_indice(1)+1:end);
+%         posicion_indice=[];
+%     %pasar string a minuscula
+%         indice_actual_lower=lower(indice_actual);
+    %asignar a cada indice un valor de indice, que seria la posición de cada
+    %indice en la lista objeto
+    %indices={'ii','msi','ndvi','ndwi','savi_03','savi_05','sr'};
+    %end
+    %</%recortar el valor de
+%indice>  
+%componer el nombre del archivo minmax a leer
+archivo_minmax_leer=strcat(directorio,proyecto,'_',string_indice);    
+%buscar para el indice en particular el valor min-max
+       %find  -maxdepth 1 -iname  "*.txt" 
+       %entre lista de archivos en directorio
+% ./ANG006_IMG_msi.txt
+% ./ANG006_IMG_ndwi.txt
+% ./ANG006_IMG_ndvi.txt
+% ./ANG006_IMG_savi_05.txt
+% ./ANG006_IMG_savi_03.txt
+% ./ANG006_IMG_ii.txt
+% ./ANG006_IMG_sr.txt
+cd(directorio);
+%No es necesario, se define en base a cada indice requerido.
+% [A L]=unix(['find  -maxdepth 1 -iname  "',archivo_minmax_leer,'.txt" ']); %entrega 0 y la lista de arriba
+% %guardar lista de nombres de archivo en:
+% files_ind= strfind(L,'.txt');
+% cantidad_archivos=length(files_ind)
+% Lista{1}=L(3:files_ind(1)-1);    
+% for i=2:cantidad_archivos
+%     Lista{i}=L(files_ind(i-1)+7:files_ind(i)-1);    
+% %leer archivos y cargar como estrcutura de listas
+% 
+% %tomar una lista adecuada y buscar dentro el archivo
+%     if strcomp(,archivo_minmax_leer)
+%         %leer archivo
+%         
+%     end    
+% %pasar el valor min y max de string a numero
+% end
+% 
+% 
+% [L_fil,L_col]=size(Lista); %entrega largo de lista, tomando L_col
+% 
+% %usar la formular x(z)=z*(b-a)/(2^n-1)) + a (obtenida en
+% %manual_conversion_indice)
+
+%se debe leer el archivo archivo_minmax_leer,'.txt"
+fid=fopen([archivo_minmax_leer,'.txt'],'r');
+%%cargar datos en listas
+coma=',';%Caracter de separacion CSV
+indice_imagen=[];
+%se genera una estructura que puede contener los valores minmax de todos
+%los indices analizados
+indice_imagen.indice{g}=string_indice;
+indice_imagen.rangos{g}(1)={minmax_i};
+while ~feof(fid)
+    leer_linea=fgetl(fid);
+    puntos_corte=strfind(leer_linea,coma);
+    indice_imagen.ID{minmax_i}=leer_linea(1:puntos_corte(1)-1);
+    indice_imagen.minimo{minmax_i}=leer_linea(puntos_corte(1)+1:puntos_corte(2)-1);
+    indice_imagen.maximo{minmax_i}=leer_linea(puntos_corte(2)+1:end);
+    minmax_i=minmax_i+1;
+end
+indice_imagen.rangos{g}(2)={minmax_i-1};%valores se rescatan con: cell2mat(ans)
+%http://www.mathworks.com/help/matlab/ref/cell2mat.html
+%esto podria ser en parelelo (para cuando se pueda)
+%ciclo de lectura de datos, leer cada imagen, buscar datos minmax de cada
+%imagen
+%usando geotiffreader.---> corte imagen
+%se tiene una lista de puntos ubicados dentro de la imagen a leer.
+%para el caso en que ahora se requieren sacar varias imagenens de una misma
+%foto no es la solucion mas eficiente, ya que cada vez que saca un nuevo
+%corte lee y carga el archivo para entregar cada una de ellas como valor
+%matricial
+%Por lo tanto hay que modificaar la funcion para que UTM_x y UTM_y sean una
+%lista de posiciones y no un solo valor, asi se puede ahorrar calculo. Se
+%propone corte_imagen_multi
+%[IMc, Rc, X ,R , INFO]=corte_imagen(BaseDir,File,UTM_x,UTM_y)
+
+
+%en directorio_tiff
+%tomando nombre de archivo de Lista_tiff
+%generar un valor correpondiente a porcentajes de analisis:
+porcentajes=[0.5,0.6,0.7,0.8];
+
+
+end %obtener nombre archivo
+
+%los valores a transformar son cada a% para hacer histograma de barras
+p_h=[=:.05:1];
+%analizar histograma y seccionar cada intervalor para entregar bloques.
+
+%guardar parametros estadistics
+
+%limpiar imagen cargada
+
+%volver a leer otra imagen

corte_imagen_multi.m

@@ -0,0 +1,113 @@
+%permite recuperar multiples zonas  a partir de una imagen.
+function [corte_multi, X ,R , INFO]=corte_imagen_multi(BaseDir,File,UTM_x,UTM_y)
+%cortar seccion definida
+%UTM_x: lista de valores para posicion X {[X1(1) X1(2)],X2(1),X2(2),....]}
+%UTM_y: lista de valores para posicion Y
+%UTM son valores de posicion reales, no de pixeles
+%obtener largo de lista----[ 1 , n ] =size(UTM_x) tomamos n por la cantidad
+%de elementos apres de posicion
+[fil, cols]=size(UTM_x);
+[X, R] = geotiffread(File);
+cd(BaseDir);
+pwd;
+if ~isdir('Cortes_Imagenes')
+mkdir('Cortes_Imagenes');
+end
+cd('Cortes_Imagenes');
+ruta_actual=pwd;
+INFO = geotiffinfo(File);
+UTM.geokey=INFO.GeoTIFFTags.GeoKeyDirectoryTag;
+UTM.tipo=INFO.Projection;
+[n,m]=size(X);
+%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
+py= INFO.GeoTIFFTags.ModelPixelScaleTag(2);
+
+%entregar los dos puntos limite del corte
+%X_a__dx_1___
+%   dy_1
+%     |     X_1(1).................
+%            .                  .  
+%            .                  .
+%            .                  .
+%            .                  .
+%            .                  .
+%            .                  .
+%            .                  .
+%            .                  .
+%..............................X_2(1)______________
+%                                             dy_2  
+%                                               |
+%                                ____________dx_2______X_b
+
+%entregar los dos puntos limite del corte
+%X_a__dx_1___
+%   dy_1
+%     |X_1(2)........
+%      .            .  
+%      .            .
+%      .            .
+%...................X_2(2)______________
+%                              
+%                              
+%                              
+%                         
+%            
+%                                             dy_2  
+%                                               |
+%                                ____________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
+
+
+%se comienza la iteracion para cada zona. recortando de la imagen inicial
+%las partes de interes
+
+for g=1:n
+    UTM_x_round=round(UTM_x{g}/p)*p;
+    UTM_y_round=round(UTM_y{g}/py)*py;   
+
+    X_1=UTM_x_round;
+    X_2=UTM_y_round;
+
+    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_round(1),1);
+    right=find(Xr==UTM_x_round(2),1);
+    upper=find(Yr==UTM_y_round(1),1);
+    lower=find(Yr==UTM_y_round(2),1);
+
+    corte_multi.Rc.YLimWorld{g}=[Yr(min(upper,lower)) Yr(max(upper,lower))]
+    corte_multi.Rc.XLimWorld{g}=[Xr(min(left,right)) Xr(max(left,right))]
+
+% [left_long, upper_lat] = projfwd(INFO, X_1(2), X_1(1))
+% [right_long, lower_lat] = projfwd(INFO, X_2(2), X_2(1))
+% [left_d, upper_d] = R.worldToIntrinsic(left_long, upper_lat);
+% left = round(abs(int64(left_d))/30);
+% upper = round(abs(int64(upper_d))/30);
+% [right_d, lower_d] = R.worldToIntrinsic(right_long, lower_lat);
+% right = round(abs(int64(right_d))/30);
+% lower = round(abs(int64(lower_d))/30);
+
+    FILAS=[ min(upper,lower) max(upper,lower)];
+    COLUMNAS=[min(left,right) max(left,right) ];
+    corte_multi.imagen{g} = X(  min(upper,lower):max(upper,lower)-1, min(left,right):max(left,right)-1  );
+    corte_multi.Rc.RasterSize{g}=size(corte_multi.imagen{g});
+    %filas en y
+    %columnas en x
+    corte_multi.Rc{g} = maprasterref('RasterSize',    corte_multi.Rc.RasterSize{g},'YLimWorld', corte_multi.Rc.YLimWorld{g}, 'ColumnsStartFrom','north','XLimWorld', corte_multi.Rc.XLimWorld{g});
+end
+
+
+end

grafica_indices_ANG006.m

@@ -58,8 +58,7 @@
 %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
-
+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); 

indices_vegetacion.m

@@ -42,7 +42,7 @@
 
     end
 end
-
+R_l=R_l+1;%% pasar de 0 a 255 a 1 a 256 (para evitar divisiones por 0 proximamente)
 IndiceVegetacion.NDVI=0;
 IndiceVegetacion.SR=0;
 IndiceVegetacion.DVI=0;
@@ -64,7 +64,7 @@
             IndiceVegetacion.NDVI=Numerador./Denominador;
             INFO_x=INFO.nir;
        case 'SR'
-           IndiceVegetacion.SR=R_l(:,:,1)./(R_l(:,:,2)+1);
+           IndiceVegetacion.SR=R_l(:,:,1)./(R_l(:,:,2));
            INFO_x=INFO.nir;
        case 'DVI'
            IndiceVegetacion.DVI=R_l(:,:,1)-R_l(:,:,2);

reflactancia.m

@@ -154,6 +154,8 @@
     end
 end
 
+
+
 %calculo indice de vegetacion NDVI: (B4-B3)/(B4+B3)
 
 end

test_indices_ANG006.m

@@ -111,13 +111,13 @@
 
 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');
+  savefile_ndvi=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'_ndvi','.txt');
+  savefile_sr=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'_sr','.txt');
+  savefile_ndwi=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'_ndwi','.txt');  
+    savefile_savi_03=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'_savi_03','.txt');
+    savefile_savi_05=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'_savi_05','.txt');
+    savefile_msi=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'_msi','.txt');
+  savefile_ii=strcat(BaseDir,'/',Carpeta_output,'/',Nombre_Proceso,'_ii','.txt');
 
 for i=1:I-1
 FID_ndvi=fopen(savefile_ndvi,'a+');
@@ -128,6 +128,8 @@
 FID_msi=fopen(savefile_msi,'a+');
 FID_ii=fopen(savefile_ii,'a+');
 
+%se guardan los valores min max de cada indice para cada imagen
+
 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})]);

test_multicorte_geotiff.m

@@ -0,0 +1 @@
+%dar ruta no

valores_histograma.m

@@ -0,0 +1,29 @@
+function [hist,Ib]=valores_histograma(Imagen,I,b)
+
+[n,m]=size(Imagen);
+ImD=double(Imagen);
+bits=2^b;
+h=zeros(1,bits);
+
+for i=1:n
+    for j=1:m
+        k = ImD(i,j);
+        h(k+1) = h(k+1)+1;
+    end
+end
+
+I_b=[];
+hist_b=[];
+k=1;
+for i=1:length(I)
+if hist(i)>0
+hist_b(k)=hist(i);
+I_b(k)=I(i);
+k=k+1;
+end
+end
+
+hist=[];
+hist=hist_b;
+
+end

vector2shape.m

@@ -1,3 +1,11 @@
+%function S_shape=vector2shape(v_iv,v_lon,v_lat,Namefile,Corte)
+%v_iv: vector de valors
+%v_lon: vector longitud asociado a valores
+%v_lat: vector latitud asociado a valores
+%Namefile: nombre de archivo
+%Corte: valor de corte
+%Esta funcion permite transfornar valores asociados geográficamente a un
+%set de archivos tipo shape para lectura en GIS.
 
 function S_shape=vector2shape(v_iv,v_lon,v_lat,Namefile,Corte)