From d6ead70b4601e0d9e8ed00a458f85e14167cacc2 Mon Sep 17 00:00:00 2001 From: dpineiden Date: Tue, 27 May 2014 11:13:08 -0400 Subject: [PATCH] =?UTF-8?q?Se=20a=C3=B1aden=20algunas=20mejoas=20en=20proc?= =?UTF-8?q?esamiento=20por=20paquetes=20de=20imagenes?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- corte_imagen.m | 8 ++-- desmembrar.m | 4 +- geotiff2shape.m | 15 ++++++ get_list_files.m | 1 + img_lands_one.m | 104 ++++++++++++++++++++++++++++++++++++++++++ indices2txt_ascii.m | 47 +++++++++++++++++++ indices_vegetacion.m | 2 +- map_vector2ascii.m | 8 ++++ test_indices_ANG006.m | 14 +++--- vector2shape.m | 18 ++++++++ 10 files changed, 208 insertions(+), 13 deletions(-) create mode 100644 geotiff2shape.m create mode 100644 img_lands_one.m create mode 100644 indices2txt_ascii.m create mode 100644 map_vector2ascii.m create mode 100644 vector2shape.m diff --git a/corte_imagen.m b/corte_imagen.m index cd4eec7..e6110b9 100644 --- a/corte_imagen.m +++ b/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)) diff --git a/desmembrar.m b/desmembrar.m index 2fa3d0a..59d82e4 100644 --- a/desmembrar.m +++ b/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 \ No newline at end of file + end diff --git a/geotiff2shape.m b/geotiff2shape.m new file mode 100644 index 0000000..eb494ab --- /dev/null +++ b/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 \ No newline at end of file diff --git a/get_list_files.m b/get_list_files.m index d44ee2d..77924c3 100644 --- a/get_list_files.m +++ b/get_list_files.m @@ -8,3 +8,4 @@ list_dir=dir(fullfile(path,type)); list_dir={list_dir.name}; out=list_dir; +end \ No newline at end of file diff --git a/img_lands_one.m b/img_lands_one.m new file mode 100644 index 0000000..086756b --- /dev/null +++ b/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 \ No newline at end of file diff --git a/indices2txt_ascii.m b/indices2txt_ascii.m new file mode 100644 index 0000000..dbea2b3 --- /dev/null +++ b/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); \ No newline at end of file diff --git a/indices_vegetacion.m b/indices_vegetacion.m index 7670189..7be3a94 100644 --- a/indices_vegetacion.m +++ b/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 diff --git a/map_vector2ascii.m b/map_vector2ascii.m new file mode 100644 index 0000000..4298aaa --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test_indices_ANG006.m b/test_indices_ANG006.m index be68737..4c263ba 100644 --- a/test_indices_ANG006.m +++ b/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]); diff --git a/vector2shape.m b/vector2shape.m new file mode 100644 index 0000000..cc99d65 --- /dev/null +++ b/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