From 55d0cb29b9374fd9c029bc5fe9ca919541014b15 Mon Sep 17 00:00:00 2001 From: dpineiden Date: Wed, 4 Jun 2014 13:12:20 -0400 Subject: [PATCH] Analisis estadistico indices --- analisis_estadistico.m | 139 ++++++++++++++++++++++++++++++++++++++ corte_imagen_multi.m | 113 +++++++++++++++++++++++++++++++ grafica_indices_ANG006.m | 3 +- indices_vegetacion.m | 4 +- reflactancia.m | 2 + test_indices_ANG006.m | 16 +++-- test_multicorte_geotiff.m | 1 + valores_histograma.m | 29 ++++++++ vector2shape.m | 8 +++ 9 files changed, 304 insertions(+), 11 deletions(-) create mode 100644 analisis_estadistico.m create mode 100644 corte_imagen_multi.m create mode 100644 test_multicorte_geotiff.m create mode 100644 valores_histograma.m diff --git a/analisis_estadistico.m b/analisis_estadistico.m new file mode 100644 index 0000000..c565a28 --- /dev/null +++ b/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 + % +%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 \ No newline at end of file diff --git a/corte_imagen_multi.m b/corte_imagen_multi.m new file mode 100644 index 0000000..3d71613 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/grafica_indices_ANG006.m b/grafica_indices_ANG006.m index 7f6e294..8d90013 100644 --- a/grafica_indices_ANG006.m +++ b/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); diff --git a/indices_vegetacion.m b/indices_vegetacion.m index 7be3a94..e404878 100644 --- a/indices_vegetacion.m +++ b/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); diff --git a/reflactancia.m b/reflactancia.m index 7d8c5ea..6f5d00f 100644 --- a/reflactancia.m +++ b/reflactancia.m @@ -154,6 +154,8 @@ end end + + %calculo indice de vegetacion NDVI: (B4-B3)/(B4+B3) end \ No newline at end of file diff --git a/test_indices_ANG006.m b/test_indices_ANG006.m index 4c263ba..3ddad1f 100644 --- a/test_indices_ANG006.m +++ b/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})]); diff --git a/test_multicorte_geotiff.m b/test_multicorte_geotiff.m new file mode 100644 index 0000000..bf6af16 --- /dev/null +++ b/test_multicorte_geotiff.m @@ -0,0 +1 @@ +%dar ruta no \ No newline at end of file diff --git a/valores_histograma.m b/valores_histograma.m new file mode 100644 index 0000000..7f11a70 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/vector2shape.m b/vector2shape.m index cc99d65..b2840b0 100644 --- a/vector2shape.m +++ b/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)