From c50a60fee5280205db1baf5dd6c2cb73e4e28362 Mon Sep 17 00:00:00 2001 From: David Pineiden Date: Wed, 8 Jan 2014 16:40:13 -0300 Subject: [PATCH] =?UTF-8?q?Se=20a=C3=B1ade=20proceso=20para=20indices=20de?= =?UTF-8?q?=20vegetacion=20INDICES=3D{'NDVI','SR','DVI','SAVI','TNDI','NDW?= =?UTF-8?q?I','EVI','TGDVI'}?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- AnalisisDirectorio.m | 15 ++--- buscar_tar.m | 18 +++--- corte_imagen.m | 11 ++-- cuadros_geograficos_zc.m | 74 ++++++++++++++++++++++++ gen_db.m | 0 indices_landsat.m | 41 +++++++++++++ indices_vegetacion.m | 121 +++++++++++++++++++++++++++++++++++++++ ndvi.m | 94 +++++++++++++++++++++++++----- reflactancia.m | 9 ++- set_MTL_to_db.m | 53 +++++++++++++++++ test_corte.m | 8 +-- test_indices.m | 11 ++++ 12 files changed, 412 insertions(+), 43 deletions(-) create mode 100644 cuadros_geograficos_zc.m create mode 100644 gen_db.m create mode 100644 indices_landsat.m create mode 100644 indices_vegetacion.m create mode 100644 set_MTL_to_db.m create mode 100644 test_indices.m diff --git a/AnalisisDirectorio.m b/AnalisisDirectorio.m index 87eb6cf..376dc88 100644 --- a/AnalisisDirectorio.m +++ b/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='img_proyecto'; -ProjectDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/img_proyecto'; +Directorio='DesAtacama'; +ProjectDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/233_076'; cd(ProjectDir); [A lista]=unix(['find -name *.txt']); RutaLog=[BaseDir,'/Logs']; @@ -27,24 +27,17 @@ %se descomprime aqui el archivo (dentro de carpeta) end fclose(fid); -archivos_MTL='archivos_MTL.txt'; +archivos_MTL=['archivos_MTL_' Directorio '.txt']; [a b]=size(MTLset); FID=fopen([RutaLog,'/',archivos_MTL],'w+'); for i=1:a fprintf(FID,MTLset{i,1}); end fclose(FID); - FID=fopen([RutaLog,'/',archivos_MTL],'r'); - savemat=0; Proceso='proyecto' - - - str='/'; - - - + str='/'; while ~feof(fid) MTL_lin=fgetl(fid); patron=MTL_lin; diff --git a/buscar_tar.m b/buscar_tar.m index 2a229ff..9040cee 100644 --- a/buscar_tar.m +++ b/buscar_tar.m @@ -1,8 +1,10 @@ BaseDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat'; -Directorio='img_ZonaControl'; -DirectorioFin='img_proyecto'; -Directorio_Padre=[BaseDir,'/',Directorio]; -[A lista]=unix(['find ',Directorio_Padre,' -name *.tar.gz']); +Directorio='DesAtacama'; +DirectorioFin='DesAtacama/Imagenes'; +Dir_tar='Imagenes/Comp'; +Directorio_Padre=[BaseDir,'/',Directorio,'/',Dir_tar]; +cd(Directorio_Padre); +[A lista]=unix(['find -name *.tar.gz']); %guardar lista en archivos archivos_targz='archivos_targz.txt'; RutaLog=[BaseDir,'/Logs']; @@ -15,7 +17,6 @@ unix(['mkdir ',BaseDir,'/',NewDir]); fid=fopen([RutaLog,'/',archivos_targz],'r'); while ~feof(fid) - cd(BaseDir); leer_linea=fgetl(FID); indice=1; nombre=desmembrar(leer_linea,'/'); @@ -25,7 +26,7 @@ indice=j; end end - cd(NewDir) + cd([BaseDir,'/',DirectorioFin]); for j=indice+1:a-1 unix(['mkdir ',nombre{j,1}]); cd(nombre{j,1}); @@ -34,8 +35,9 @@ unix(['mkdir ',nombre_toma(1:length(nombre_toma)-7)]); cd(nombre_toma(1:length(nombre_toma)-7)); %se han creado las carpetas + line_tar=[BaseDir,'/',Directorio,'/',Dir_tar,leer_linea(2:length(leer_linea))]; %se descomprime aqui el archivo (dentro de carpeta) - unix(['tar -xvf ',leer_linea]); - + unix(['tar -xf ',line_tar]); + unix(['rm ',line_tar]); end fclose(fid); \ No newline at end of file diff --git a/corte_imagen.m b/corte_imagen.m index 37b62d7..c4175ec 100644 --- a/corte_imagen.m +++ b/corte_imagen.m @@ -1,5 +1,5 @@ -function [IMc X R INFO]=corte_imagen(BaseDir,Ruta,File,nombre_corte,UTM_x,UTM_y,largo_pixel) +function [IMc, X ,R , INFO]=corte_imagen(BaseDir,File,UTM_x,UTM_y) [X, R] = geotiffread(File); %cortar seccion definida cd(BaseDir); @@ -13,6 +13,7 @@ %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 %entregar los dos puntos limite del corte %X_a__dx_1___ @@ -35,10 +36,10 @@ X_1=UTM_x; X_2=UTM_y; -left=abs(x(1)-X_1(1))/p -right=abs(x(1)-X_1(2))/p -upper=abs(y(1)-X_2(1))/p -lower=abs(y(1)-X_2(2))/p +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)) diff --git a/cuadros_geograficos_zc.m b/cuadros_geograficos_zc.m new file mode 100644 index 0000000..3a22fd9 --- /dev/null +++ b/cuadros_geograficos_zc.m @@ -0,0 +1,74 @@ +cuadro_zc.nombre={ + 'Río Loa Alto'; + 'Salar de Ascotan'; + 'Geyser del Tatio'; + 'Oasis de Atacama'; + 'Volcán Licancabur'; + 'Ayllus de San Pedro de Atacama'; + 'Salar de Atacama'; + 'Salar de Tara'; + 'Salar de aguas Calientes'; + 'Sistema de Lagunas Miscanti y Miñiques'; + 'Salar de Punta Negra'; + 'Salar de Pujsa'; + 'Laguna Lejía'; + 'Sistema Hidrológico Soncor' + }; + +%valores en wgs85 +cuadro_zc.x={ +[-68.85320,-68.09895] +[-68.26799,-67.85977] +[-68.08994,-67.81551] +[-68.30257,-68.15484] +[-67.78503,-67.60751] +[-68.65887,-68.02238] +[-69.06626,-68.88556] +[-68.97034,-68.86352] +[-70.60807,-70.36595] +[-70.58620,-70.44180] +[-70.63696,-70.48724] +[-70.09172,-69.93866] +[-69.56745,-69.49855] +[-68.40515,-68.18795] +[-68.67859,-68.53386] +}; + + +cuadro_zc.y={ +[-21.08443,-21.88897] +[-22.26893,-22.59543] +[-22.48025,-22.97415] +[-22.83610,-22.98443] +[-23.34728,-23.54745] +[-22.93933,-23.86970] +[-24.34012,-24.75150] +[-22.41437 ,-22.50059] +[-24.67875 ,-25.22364] +[-24.45511 ,-24.59138] +[-23.00353 ,-23.54435] +[-21.40138 ,-21.54488] +[-21.61380 ,-21.73473] +[-21.41117 ,-21.72944] +[-24.88029 ,-25.12285] +}; + +[n m]=size(cuadro_zc.x); +utmzone=19; +hemi='N'; +for i=1:n + Lat=cuadro_zc.y{i}(1); + Lon=cuadro_zc.x{i}(1); +[x,y,utmzone,utmhemi] = wgs2utm(Lat,Lon,utmzone,hemi); +cuadro_zc.UTMx{i}(1)=x; +cuadro_zc.UTMy{i}(1)=y; + Lat=cuadro_zc.y{i}(2); + Lon=cuadro_zc.x{i}(2); +[x,y,utmzone,utmhemi] = wgs2utm(Lat,Lon,utmzone,hemi); +cuadro_zc.UTMx{i}(2)=x; +cuadro_zc.UTMy{i}(2)=y; +end +%pasar información de recuadros de zona de control a database +%diccionario METADATA LANDSAT: +%http://glovis.usgs.gov/ImgViewer/landsat_dictionary.html + diff --git a/gen_db.m b/gen_db.m new file mode 100644 index 0000000..e69de29 diff --git a/indices_landsat.m b/indices_landsat.m new file mode 100644 index 0000000..a80efc7 --- /dev/null +++ b/indices_landsat.m @@ -0,0 +1,41 @@ +%buscar cada cuadro en la lista de imagenes de proyecto +%entregar cuatro ´indices en los que encuentra cada punto. +%esto indica bajo que imagenes landsat se encuentra el recuadro +%Puntos_x es un par de valores que indican minimo y maximo horizontal +%Puntos_y un par de valores que indican minimo maximo vertical +%Arreglo_x es una lista de valores de dos columnas que contiene los valores +%minimo y maximo horizontal de las imagenes landsat en UTM 19N +%Arreglo_y es una lista de valores de dos columnas que contiene los valores +%minimo y maximo vertical de las imagenes landsat en UTM 19N +function [x1 x2 y1 y2]=indices_landsat(Puntos_x,Puntos_y,Arreglo_x,Arreglo_y) +[n m]=size(Arreglo_x);%el mismo que Arreglo_y n filas y m=2 columnas +%tamaño de Puntos_- es de 1 fila, 2 columnas +j=1; +l=1; +h=1; +k=1; +x1=0; +x2=0; +y1=0; +y2=0; +for i=1:n + if and(Arreglo_y(i,1)<=Puntos_y(1),Puntos_y(1)<=Arreglo_y(i,2)) + y1(h)=i; + h=h+1; + end + if and(Arreglo_y(i,1)<=Puntos_y(2),Puntos_y(2)<=Arreglo_y(i,2)) + y2(k)=i; + k=1+k; + end + if and(Arreglo_x(i,1)<=Puntos_x(1),Puntos_x(1)<=Arreglo_x(i,2)) + x1(j)=i; + j=1+j; + end + if and(Arreglo_x(i,1)<=Puntos_x(2),Puntos_x(2)<=Arreglo_x(i,2)) + x2(l)=i; + l=1+l; + end + + +end +end \ No newline at end of file diff --git a/indices_vegetacion.m b/indices_vegetacion.m new file mode 100644 index 0000000..6ad78b5 --- /dev/null +++ b/indices_vegetacion.m @@ -0,0 +1,121 @@ + +%definir directorio de trabajo +%setear en este directorio +%Nombre del Proceso +%se desean guardar las variables en .mat?savemat 0:no, 1:si +%INDICES={"NDVI","SR","DVI","SAVI","TNDI","NDWI","EVI","TGDVI"} +function [NDVI,SR,DVI,SAVI,TNDI,NDWI,EVI,TGDVI]=indices_vegetacion(BaseDir,MTLDIR,savemat,Nombre_Proceso,CorteX,CorteY,INDICES) +[n m]=size(INDICES); +cd(BaseDir); +var_mat=savemat; +Carpeta_output='Output';%salida, guardar resultados +i=1; +Landsat=5; +Carpeta_img=strcat(BaseDir,'/img_',Nombre_Proceso); +%Nombres bandas:[blue,green,red,NIR,SWIR,-] +hemisferio='S'; +%reflactancia NIR 1 +Banda=4; +[R_l(:,:,1),L_l(:,:,1),DN7(:,:,1),UTM,COD_SAT_IMG,INFO.nir]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); +%reflactancia ROJA 2 +Banda=3; +[R_l(:,:,2),L_l(:,:,2),DN7(:,:,2),UTM,COD_SAT_IMG,INFO.red]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); +%reflactancia AZUL 3 +Banda=1; +[R_l(:,:,3),L_l(:,:,3),DN7(:,:,3),UTM,COD_SAT_IMG,INFO.blue]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); +%reflactancia SWIR 4 +Banda=5; +[R_l(:,:,4),L_l(:,:,4),DN7(:,:,4),UTM,COD_SAT_IMG,INFO.swir]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); +NDVI=0; +SR=0; +DVI=0; +SAVI=0; +TNDI=0; +NDWI=0; +EVI=0; +TGDVI=0; +for i=1:m +INDICES{1,i} + switch INDICES{1,i} + case 'NDVI' + Numerador=R_l(:,:,1)-R_l(:,:,2); + Denominador=R_l(:,:,1)+R_l(:,:,2); + NDVI=Numerador./Denominador; + INFO_x=INFO.nir; + Indice=NDVI; + case 'SR' + SR=R_l(:,:,1)./R_l(:,:,2); + INFO_x=INFO.nir; + Indice=SR; + case 'DVI' + DVI=R_l(:,:,1)-R_l(:,:,2); + INFO_x=INFO.nir; + Indice=DVI; + case 'SAVI' + L=0.5; + Numerador=R_l(:,:,1)-R_l(:,:,2); + Denominador=R_l(:,:,1)+R_l(:,:,2)+L; + SAVI=(Numerador./Denominador)*(1+L); + INFO_x=INFO.nir; + Indice=SAVI; + case 'TNDI' + Numerador=R_l(:,:,1)-R_l(:,:,2); + Denominador=R_l(:,:,1)+R_l(:,:,2); + pre_TNDI=(Numerador./Denominador) +1; + TNDI=sqrt(pre_TNDI); + INFO_x=INFO.nir; + Indice=TNDI; + case 'NDWI' + Numerador=R_l(:,:,1)-R_l(:,:,4); + Denominador=R_l(:,:,1)+R_l(:,:,4); + NDWI=Numerador./Denominador; + INFO_x=INFO.swir; + Indice=NDWI; + case 'EVI' + G=2.5; + L=1; + C_1=6; + C_2=7.5; + Numerador=R_l(:,:,1)-R_l(:,:,2); + Denominador=R_l(:,:,1)+C_1*R_l(:,:,2)-C_2*R_l(:,:,3)+L; + EVI=G*Numerador./Denominador; + INFO_x=INFO.nir; + Indice=EVI; + case 'TGDVI' + La_1=0.56e-6;%en metros + La_2=0.66e-6; + La_3=0.83e-6; + Numerador_1=R_l(:,:,1)-R_l(:,:,2); + Denominador_1=La_3-La_2; + Numerador_2=R_l(:,:,2)-R_l(:,:,3); + Denominador_2=La_2-La_1; + pre_TGDVI=Numerador_1./Denominador_1-Numerador_2./Denominador_2; + %pasar a 0 valores menores que 0 + pre_TGDVI(pre_TGDVI<0)=0; + TGDVI=pre_TGDVI; + INFO_x=INFO.blue; + Indice=TGDVI; + end + + +% I=mat2gray(Indice,[min(min(Indice)) max(max(Indice))]); +% [I_16b,Map_16b]=gray2ind(I,65536); +% Histograma=imhist(I_16b); +% Bit_depth=16;%'BitDepth'= 16 bits +% NDVI_16=uint16(NDVI*65535); +% 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); +end + +if var_mat==1 +savefile=strcat(Carpeta_output,'/',Nombre_Proceso,'/',COD_SAT_IMG,'.mat'); +save(savefile,'R_l','L_l','DN7','UTM'); +end + + +end \ No newline at end of file diff --git a/ndvi.m b/ndvi.m index 13c6d88..ba20688 100644 --- a/ndvi.m +++ b/ndvi.m @@ -3,37 +3,105 @@ %setear en este directorio %Nombre del Proceso %se desean guardar las variables en .mat?savemat 0:no, 1:si -function NDVI=ndvi(BaseDir,MTLDIR,savemat,Nombre_Proceso) +%INDICES={"NDVI","SR","DVI","SAVI","TNDI","NDWI","EVI","TGDVI"} +function [NDVI,SR,DVI,SAVI,TNDI,NDWI,EVI,TGDVI]=indices_vegetacion(BaseDir,MTLDIR,savemat,Nombre_Proceso,CorteX,CorteY,INDICES) +[n m]=size(INDICES) +%for --> case cd(BaseDir); var_mat=savemat; Carpeta_output='Output';%salida, guardar resultados i=1; Landsat=5; Carpeta_img=strcat(BaseDir,'/img_',Nombre_Proceso); +%Nombres bandas:[blue,green,red,NIR,SWIR,-] hemisferio='S'; +%reflactancia NIR 1 Banda=4; -[R_l(:,:,1),L_l(:,:,1),DN7(:,:,1),UTM,COD_SAT_IMG,INFO]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso); +[R_l(:,:,1),L_l(:,:,1),DN7(:,:,1),UTM,COD_SAT_IMG,INFO]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); +%reflactancia ROJA 2 Banda=3; -[R_l(:,:,2),L_l(:,:,2),DN7(:,:,2),UTM,COD_SAT_IMG,INFO]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso); -Numerador=R_l(:,:,2)-R_l(:,:,1); -Denominador=R_l(:,:,2)+R_l(:,:,1); -valor_NDVI=Numerador./Denominador; -if var_mat==1 -savefile=strcat(Carpeta_output,'/',Nombre_Proceso,'/',COD_SAT_IMG,'.mat'); -save(savefile,'R_l','L_l','DN7','UTM'); -end +[R_l(:,:,2),L_l(:,:,2),DN7(:,:,2),UTM,COD_SAT_IMG,INFO]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); +%reflactancia AZUL 3 +Banda=1; +[R_l(:,:,3),L_l(:,:,3),DN7(:,:,3),UTM,COD_SAT_IMG,INFO]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); +%reflactancia SWIR 4 +Banda=5; +[R_l(:,:,4),L_l(:,:,4),DN7(:,:,4),UTM,COD_SAT_IMG,INFO]=reflactancia(Banda,Landsat,Carpeta_img,MTLDIR,hemisferio,Nombre_Proceso,CorteX,CorteY); +NDVI=0; +SR=0; +DVI=0; +SAVI=0; +TNDI=0; +NDWI=0; +EVI=0; +TGDVI=0; +for i=1:n + switch INDICES{i} + case 'NDVI' + Numerador=R_l(:,:,1)-R_l(:,:,2); + Denominador=R_l(:,:,1)+R_l(:,:,1); + NDVI=Numerador./Denominador; + case 'SR' + SR=R_l(:,:,1)./R_l(:,:,2) + case 'DVI' + DVI=R_l(:,:,1).-R_l(:,:,2) + case 'SAVI' + L=0.5; + Numerador=R_l(:,:,1)-R_l(:,:,2); + Denominador=R_l(:,:,1)+R_l(:,:,2)+L; + SAVI=(Numerador./Denominador)*(1+L); + case 'TNDI' + Numerador=R_l(:,:,2)-R_l(:,:,1); + Denominador=R_l(:,:,2)+R_l(:,:,1); + pre_TNDI=(Numerador./Denominador) +0.5; + TNDI=sqrt(pre_TNDI) + case 'NDWI' + Numerador=R_l(:,:,2)-R_l(:,:,4); + Denominador=R_l(:,:,2)+R_l(:,:,4); + NDVI=Numerador./Denominador; + case 'EVI' + G=2.5; + L=1; + C_1=6; + C_2=7.5; + Numerador=R_l(:,:,1)-R_l(:,:,2); + Denominador=R_l(:,:,1)+C_1*R_l(:,:,4)-C_2*R_l(:,:,3)+L; + NDVI=G*Numerador./Denominador; + case 'TGDVI' + La_1=0.56e-6;%en metros + La_2=0.66e-6; + La_3=0.83e-6; + Numerador_1=R_l(:,:,1)-R_l(:,:,2); + Denominador_1=La_3-La_2; + Numerador_2=R_l(:,:,2)-R_l(:,:,3); + Denominador_2=La_2-La_1; + pre_TGDVI=Numerador_1./Denominador_1-Numerador_2./Denominador_2; + %pasar a 0 valores menores que 0 + pre_TGDVI(pre_TGDVI<0)=0; + TGDVI=pre_TGDVI; + + end + + I=mat2gray(valor_NDVI,[-1 1]); [I_16b,Map_16b]=gray2ind(I,65536); Histograma=imhist(I_16b); Bit_depth=16;%'BitDepth'= 16 bits NDVI_16=uint16(valor_NDVI*65535); -savefile_TIF=strcat(Carpeta_output,'/',Nombre_Proceso,'/',COD_SAT_IMG,'.TIF'); +savefile_TIF=strcat(Carpeta_output,'/',Nombre_Proceso,'/',COD_SAT_IMG,'_',INDICES{i},'.TIF'); %geotiffwrite(savefile_TIF,valor_NDVI,INFO,'GeoKeyDirectoryTag',UTM.geokey); geotiffwrite(savefile_TIF,I_16b,INFO,'GeoKeyDirectoryTag',UTM.geokey); %imwrite(I_16b,savefile_TIF,'tif'); %Info_tif=imfinfo(savefile_TIF); -savefile_PNG=strcat(Carpeta_output,'/',Nombre_Proceso,'/',COD_SAT_IMG,'.PNG'); +savefile_PNG=strcat(Carpeta_output,'/',Nombre_Proceso,'/',COD_SAT_IMG,'_',INDICES{i},'.PNG'); imwrite(I_16b,savefile_PNG,'png'); Info_png=imfinfo(savefile_PNG); -NDVI='hola'; end + +if var_mat==1 +savefile=strcat(Carpeta_output,'/',Nombre_Proceso,'/',COD_SAT_IMG'.mat'); +save(savefile,'R_l','L_l','DN7','UTM'); +end + + +end \ No newline at end of file diff --git a/reflactancia.m b/reflactancia.m index fde4d02..456c194 100644 --- a/reflactancia.m +++ b/reflactancia.m @@ -7,7 +7,7 @@ %entregar Path %Entregar hemisferio -function [R_l,L_l,DN7,UTM,COD_SAT_IMG,R]=reflactancia(Banda,Landsat,Path,MTLDIR,Hemisferio,Nombre_Proceso) +function [R_l,L_l,DN7,UTM,COD_SAT_IMG,R]=reflactancia(Banda,Landsat,Path,MTLDIR,Hemisferio,Nombre_Proceso,CorteX,CorteY) cd(Path); cd('..'); BaseDir=pwd; @@ -97,11 +97,16 @@ Busqueda=strcat('*_B',num2str(Banda),'.TIF'); archivo=get_list_files(Carpeta_img,Busqueda); FilePath=char(strcat(Carpeta_img,'/',archivo)); +INFO = geotiffinfo(FilePath); + +if ~strcmp(CorteX,'') | ~strcmp(CorteY,'') +[IMc X R INFO]=corte_imagen(BaseDir,File,CorteX,CorteY); +else [X, R] = geotiffread(FilePath); +end %X es la matriz con valores 0 a 255 8bit, grayscale= 2 dim %CMAP es el mapa de colores %R es la referencia geoespacial -INFO = geotiffinfo(FilePath); UTM.geokey=INFO.GeoTIFFTags.GeoKeyDirectoryTag; UTM.tipo=INFO.Projection; %pasa coordenadas a deg utm 19 N a deg diff --git a/set_MTL_to_db.m b/set_MTL_to_db.m new file mode 100644 index 0000000..3beaf07 --- /dev/null +++ b/set_MTL_to_db.m @@ -0,0 +1,53 @@ + +cd('/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076'); +BaseProject=pwd; +%buscar UL y LR en cada MTL de las imagenes de la zona +MTL_name='LANDSAT_SCENE_ID'; +%se declaran los string a buscar +UL_x='CORNER_UL_PROJECTION_X_PRODUCT'; +UL_y='CORNER_UL_PROJECTION_Y_PRODUCT'; +LR_x='CORNER_LR_PROJECTION_X_PRODUCT'; +LR_y='CORNER_LR_PROJECTION_Y_PRODUCT'; +%tomar archivo MTL de proyecto +Directorio='DesAtacama'; +archivos_MTL=['archivos_MTL_' Directorio '.txt']; +%pasar a basedatos NOMBRE imagen, buscar los UL LR y guardarlos en db +BaseLog='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/Logs'; +FID_MTL=fopen([BaseLog,'/',archivos_MTL],'r'); +str='='; +i=0; +UTM_MTL_x=[]; +UTM_MTL_y=[]; +FILE_name={}; +while ~feof(FID_MTL) +i=1+i; +MTL_lin=fgetl(FID_MTL); +FILE_name{i}=[BaseProject MTL_lin(2:length(MTL_lin))]; +fid=fopen(MTL_lin,'r'); +while ~feof(fid) +lin=fgetl(fid); +%leer archivo de cada línea MTL +%buscar '=' +patron=lin; +k_str=strfind(patron,str); +n= length(k_str); +if n>=1 +%nombre y valor: +UTM_name=lin(5:k_str(n)-2); +UTM_val=lin(k_str(length(k_str))+2:length(lin)); +%comparar por caso +switch UTM_name + case UL_x + UTM_MTL_x(i,1)=str2double(UTM_val); + case UL_y + UTM_MTL_y(i,1)=str2double(UTM_val); + case LR_x + UTM_MTL_x(i,2)=str2double(UTM_val); + case LR_y + UTM_MTL_y(i,2)=str2double(UTM_val); +end%switch +end +end%while +fclose(fid); +end +fclose(FID_MTL); \ No newline at end of file diff --git a/test_corte.m b/test_corte.m index e49bbfd..5391e10 100644 --- a/test_corte.m +++ b/test_corte.m @@ -2,15 +2,15 @@ %cargar imagen %definir ruta de imagen BaseDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/' -File='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/img_test/L5Ssat1/LT50010751985064XXX04_B1.TIF'; -Ruta='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/img_test/L5Ssat1'; +File='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/1984/12_diciembre/LM52330761984359AAA03_B1.TIF'; +Ruta='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/1984/12_diciembre'; nombre_corte='Prueba'; UTM_x=[420000,500000]; UTM_y=[-2300000,-2400000]; largo_pixel=30; -[IMc X R INFO]=corte_imagen(BaseDir,Ruta,File,nombre_corte,UTM_x,UTM_y,largo_pixel); +[IMc, X, R, INFO]=corte_imagen(BaseDir,File,UTM_x,UTM_y); figure imshow(IMc) figure imshow(X) -save_geotif=guardar_img_tif(IMc,R,INFO,File,nombre_corte,UTM_x,UTM_y) \ No newline at end of file +save_geotif=guardar_img_tif(IMc,R,INFO,File,nombre_corte,UTM_x,UTM_y); \ No newline at end of file diff --git a/test_indices.m b/test_indices.m new file mode 100644 index 0000000..07550f4 --- /dev/null +++ b/test_indices.m @@ -0,0 +1,11 @@ +INDICES={'NDVI','SR','DVI','SAVI','TNDI','NDWI','EVI','TGDVI'}; +%INDICES={'EVI','TGDVI'}; +BaseDir='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat' +File='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2001/11_noviembre/LM52330761984359AAA03_B1.TIF'; +Ruta='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2001/11_noviembre'; +Nombre_Proceso='Prueba_Indices'; +MTLDIR='/home/david/Documents/Proyectos_CEA/CNM008/Codigo_Mat/DesAtacama/Imagenes/233_076/2001/11_noviembre'; +savemat=0; +tic +[NDVI,SR,DVI,SAVI,TNDI,NDWI,EVI,TGDVI]=indices_vegetacion(BaseDir,MTLDIR,savemat,Nombre_Proceso,'','',INDICES); +toc \ No newline at end of file