Este cuaderno es una guía para facilitar el desarrollo de la actividades en clase de los cursos: Hidroclimatológia (pregrado y posgrado) y Modelación de procesos hidrológicos de la Facultad de Estudios Ambientales y Rurales de la Pontificia Universidad Javeriana, y fue consolidado con la ayuda de Jhon Arévalo-Pulgarín, estudiante de la Maestría en Hidrosistemas en la PUJ.

Iniciamos Cargando las librerías

require("raster") #Para que se puedan leer archivos raster
## Loading required package: raster
## Loading required package: sp
require("RSAGA") #Para que R articule con SAGA-GIS
## Loading required package: RSAGA
## Loading required package: gstat
## Loading required package: shapefiles
## Loading required package: foreign
## 
## Attaching package: 'shapefiles'
## The following objects are masked from 'package:foreign':
## 
##     read.dbf, write.dbf
## Loading required package: plyr
require("terra")
## Loading required package: terra
## terra 1.7.23
require("gdalraster") #reemplazo de rgdal
## Loading required package: gdalraster
## GDAL 3.5.2, released 2022/09/02, PROJ 8.2.1
## 
## Attaching package: 'gdalraster'
## The following object is masked from 'package:terra':
## 
##     rasterize
## The following objects are masked from 'package:raster':
## 
##     calc, rasterize
require("lubridate")
## Loading required package: lubridate
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:terra':
## 
##     intersect, union
## The following objects are masked from 'package:raster':
## 
##     intersect, union
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
require("gstat")
require("hydroTSM")
## Loading required package: hydroTSM
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following object is masked from 'package:terra':
## 
##     time<-
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'hydroTSM'
## The following object is masked from 'package:terra':
## 
##     extract
## The following object is masked from 'package:raster':
## 
##     extract
require("sf")
## Loading required package: sf
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

Cargamos el DEM y generamos la gráfica para asegurarnos que este bien cargada

dem=raster("dem_3116.tif")
projection(dem)
## [1] "+proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs"
spplot(dem, col.regions=terrain.colors(255))

Ahora cargamos la límite de la cuenca y hacemos la transformación de las coordenadas

cuenca=shapefile("Cuenca_SAGA.shp")
projection(cuenca)
## [1] "+proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs"
class(cuenca)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
#projection(cuenca) = crs("+init=epsg:3116")

#cuenca= crs("EPSG:3116")
projection(dem) = crs("+init=epsg:3116")
cuenca <- spTransform(cuenca, crs("+init=epsg:3116"))
## Warning: PROJ support is provided by the sf and terra packages among others

1. Estimación de las Características Morfométricas Generales

Antes de iniciar los cálculos, asegurese de tener los archivos en coordenadas planas, pues de lo contrario no podrá calcular el valor del área y el perímetro.

A= area(cuenca)/1000000 # undiades en Km^2
y=vect(cuenca)
P=perim(y)/1000 # Km
cat("El Área de la zona hidrográfica es",A, "Km^2 y el perímetro es igual a: ",P, "Km")
## El Área de la zona hidrográfica es 4556.384 Km^2 y el perímetro es igual a:  478.33 Km
if (A<5) {
  C="Unidad"
}else if(A>=5 & A<20){
  C="Sector"
}else if(A>=20 & A<100){
  C="Microcuenca"
}else if(A>=100 & A<300){
  C="Sub-cuenca"
}else if(A>=300){
  C="Cuenca"
}

cat(" Conforme al Área (en Km^2) la zona hidrográfica se clasifica como: ", C)
##  Conforme al Área (en Km^2) la zona hidrográfica se clasifica como:  Cuenca

2. Estimación de las Características Morfométricas - Relativas a la Forma

2.1. Indice de compacidad o de Gravelius (Kc)

Kc=0.28*(P/sqrt(A))

if (Kc>=1 & Kc<1.25) {
  C_Kc="Forma redonda a oval redonda"
}else if(Kc>=1.25 & Kc<1.5){
  C_Kc="Forma oval redonda a oval oblonga"
}else if(Kc>=1.5 & Kc<1.75){
  C_Kc="Forma oval oblonga a rectangular oblonga"
}else if(Kc>=1.75){
  C_Kc="Forma rectangular oblonga"
}

cat("De acuerdo con el Indice de Compacidad o de Gravelius (Kc) la ",C," tiene ", C_Kc)
## De acuerdo con el Indice de Compacidad o de Gravelius (Kc) la  Cuenca  tiene  Forma rectangular oblonga

En este código no se contempla el cálculo de la longitud axial, Vertientes y Ancho mC!ximo, ya que estos deben ser tomados directamente del mapa. Con estos datos puede estiamr de manera adicional el factor forma, C-ndice de sinuosidad, alargamiento y asimétrico

2.2. Lados del rectangulo equivalente (RE)

Para poder calcular el rectangulo equivalente, el Kc debe ser mayor o igual a 1.12, por lo cual:

if (Kc>=1.12) {
  L1=((Kc*sqrt(A))/1.12)*(1-(sqrt(1-(1.12/Kc))^2))
  L2=((Kc*sqrt(A))/1.12)*(1+(sqrt(1-(1.12/Kc))^2))
}else{
  L1="No es posible el cálculo, Kc < 1.12"
  L2="No es posible el cálculo, Kc < 1.12"
}

3. Estimación de las Características Morfométricas - Relativos al relieve

Para esta etapa, iniciamos con la mascara del DEM con la cuenca.

dem = mask(dem,cuenca) ## Aquí no hacemos la máscara pues ya esta recortado el DEM con el límite de la Cuenca
ext = c(cuenca@bbox[1,1] ,cuenca@bbox[1,2],cuenca@bbox[2,1] ,cuenca@bbox[2,2])
dem = crop(dem,ext)
spplot(dem, col.regions=terrain.colors(255))

3.1. Elevación

cat("La elevación mínima es:",dem@data@min,"m.s.n.m.")
## La elevación mínima es: 0 m.s.n.m.
cat("La elevación media es:",mean(na.omit(dem@data@values)),"m.s.n.m.")
## La elevación media es: 2428.701 m.s.n.m.
cat("La elevación máxima es:",dem@data@max,"m.s.n.m.")
## La elevación máxima es: 5294.35 m.s.n.m.

3.1.1. Curva hipsométrica

sgdf <- as(dem, 'SpatialGridDataFrame')
jpeg("Curva_hypso.jpeg", quality = 75)
par(mar=c(5, 5, 4, 3)) # abajo, izquierda, arriba y derecha
#  histograma de altitudes
histograma <- hist(values(dem), breaks = 100, plot = FALSE)

#  curva hipsometrica acumulando los conteos de cada clase
curva_hipsometrica <- cumsum(histograma$counts) / sum(histograma$counts)
plot(histograma$mids, curva_hipsometrica, type = "l", xlab = "Altitud", ylab = "Porcentaje acumulado de area", main = "Curva hipsometrica")

dev.off()
## png 
##   2

3.2. Pendiente

pend=terrain(dem, 'slope', unit='degrees', neighbors=8) # da en %

pend_c=mean(na.omit(pend@data@values))

cat("La pendiente media es:",pend_c,"%")
## La pendiente media es: 21.20521 %
if (pend_c>=0 & pend_c<3) {
  pend_clas="Plano"
}else if(pend_c>=3 & pend_c<7){
  pend_clas="Suave"
}else if(pend_c>=7 & pend_c<12){
  pend_clas="Medianamente accidentado"
}else if(pend_c>=12 & pend_c<20){
  pend_clas="Accidentado"
}else if(pend_c>=20 & pend_c<35){
  pend_clas="Fuertemente accidentado"
}else if(pend_c>=35 & pend_c<50){
  pend_clas="Muy fuertemente accidentado"
}else if(pend_c>=50 & pend_c<75){
  pend_clas="Escarpado"
}else if(pend_c>=75){
  pend_clas="Muy escarpado"
}

cat("Conforme con la pendiente media (en %) la zona tiene un tipo de relieve:",pend_clas)
## Conforme con la pendiente media (en %) la zona tiene un tipo de relieve: Fuertemente accidentado

3.3. Curvas de Nivel

curvas=rasterToContour(dem,nlevels=100) # Así se generan cada 50 m

4. Estimación de las Características Morfométricas - Relativos a la red de drenaje

Aquí inciamos a trabajar con SAGA desde Rstudio. Para ello debemos cargar el ambiente de SAGA. Para esto, por favor revise las instrucciones dadas en el script de Delimtiación de Cuencas con SAGA.

Primero deberá descargar el programa SAGA de la página de internet Link y seguir las instrucciones de instalación. Se recomienda que se instale la herramienta en C://Program Files.

Para que R corra SAGA-GIS se debe direccionar el comando a que lea la ubicación del programa. Adicionalmente, es necesario habilitar la librería con el código que se muestra a continuación.Ejecute estas líneas ÚNICAMENTE cuando SAGA este instalado en el computador.

envi = rsaga.env(path = "C:/Program Files/SAGA") # Direccion del ejecutable de SAGA GIS
## Verify specified path to SAGA command line program...
## Found SAGA command line program. Search for not specified SAGA modules path... 
## Done
rsaga.get.version(env = envi)
## [1] "8.1.1"
# Exportar raster
writeRaster(dem,"dem_cuenca.tif",overwrite=TRUE) # Este, porque ya esta recortado a la forma de la cuenca, (lo lee de la carpeta de trabajo)

# Crea archivo .sgrd para trabajar con el dem en SAGA
rsaga.geoprocessor(lib = "io_gdal", 0,
                   param = list(GRIDS ="dem_sgrd",
                                FILES= "dem_cuenca.tif",
                                TRANSFORM= 1),env = envi)
## 
## 
## SAGA Version: 8.1.1 (64 bit)
## 
## library path: C:\PROGRA~1\SAGA\tools\
## library name: io_gdal
## library     : io_gdal
## tool        : Import Raster
## identifier  : 0
## author      : O.Conrad (c) 2007 (A.Ringeler)
## processors  : 8 [8]
## 
## [Import Raster] Execution started...
## 
## [Import Raster] Parameters:
## 
## Grids: No objects
## Files: "dem_cuenca.tif"
## Multiple Bands Output: automatic
## Select from Multiple Bands: 
## Transformation: true
## Resampling: Nearest Neighbour
## Extent: original
## 
## 
## loading: dem_cuenca.tif
## 
## 
## 
## Driver: GTiff
## 
## Bands: 1
## 
## Rows: 2232
## 
## Columns: 2962
## 
## 
## loading: dem_cuenca
## 
## total execution time: 0 milliseconds (less than 1 millisecond)
## 
## [Import Raster] Execution succeeded (less than 1 millisecond)
## Saving grid: dem_sgrd...
# Llenar datos vacC-os
rsaga.geoprocessor(lib = "ta_preprocessor", 2,
                   param = list(DEM ="dem_sgrd.sgrd",
                                DEM_PREPROC= "demsink_fill"),env = envi)
## 
## 
## SAGA Version: 8.1.1 (64 bit)
## 
## library path: C:\PROGRA~1\SAGA\tools\
## library name: ta_preprocessor
## library     : ta_preprocessor
## tool        : Sink Removal
## identifier  : 2
## author      : O. Conrad (c) 2001
## processors  : 8 [8]
## 
## Loading grid: dem_sgrd.sgrd...
## [Sink Removal] Execution started...
## 
## [Sink Removal] Parameters:
## 
## Grid System: 31; 2232x 2962y; 751970.713617x 723106.370313y
## DEM: dem_sgrd
## Sink Route: <not set>
## Preprocessed DEM: Preprocessed DEM
## Method: Fill Sinks
## Threshold: false
## 
## Create index: dem_sgrd [no sinks]
## Find Pits
## Find Outlets
## Routing
## Finalize
## 
## number of processed sinks: 15404
## Initializing direction matrix...
## I'm fillin'...
## 
## total execution time: 20000 milliseconds (20s)
## 
## [Sink Removal] Execution succeeded (20s)
## Saving grid: demsink_fill...
# Calcular red de drenajes, lo hace por Strahler
rsaga.get.modules("ta_channels",env = envi)
## $ta_channels
##   code                                      name interactive
## 1    0                           Channel Network       FALSE
## 2    1                          Watershed Basins       FALSE
## 3    2               Watershed Basins (Extended)       FALSE
## 4    3      Vertical Distance to Channel Network       FALSE
## 5    4 Overland Flow Distance to Channel Network       FALSE
## 6    5       Channel Network and Drainage Basins       FALSE
## 7    6                            Strahler Order       FALSE
## 8    7                              Valley Depth       FALSE
rsaga.get.usage('ta_channels', 5,env = envi)
## library path: C:\PROGRA~1\SAGA\tools\
## library name: ta_channels
## library     : ta_channels
## Usage: saga_cmd ta_channels 5 [-DEM <str>] [-DIRECTION <str>] [-CONNECTION <str>] [-ORDER <str>] [-BASIN <str>] [-SEGMENTS <str>] [-BASINS <str>] [-NODES <str>] [-THRESHOLD <num>]
##   -DEM:<str>         Elevation
##  Grid, input
##   -DIRECTION:<str>   Flow Direction
##  Grid, output, optional
##   -CONNECTION:<str>  Flow Connectivity
##  Grid, output, optional
##   -ORDER:<str>       Strahler Order
##  Grid, output, optional
##   -BASIN:<str>       Drainage Basins
##  Grid, output, optional
##   -SEGMENTS:<str>    Channels
##  Shapes, output
##   -BASINS:<str>      Drainage Basins
##  Shapes, output
##   -NODES:<str>       Junctions
##  Shapes, output, optional
##   -THRESHOLD:<num>   Threshold
##  Integer
##  Minimum: 1
rsaga.geoprocessor(lib = "ta_channels", 5,
                   param = list(DEM ="demsink_fill.sgrd",
                                ORDER="Ordenes_Strahler_raster",
                                SEGMENTS="Drenaje_ordenes",
                                THRESHOLD=5), # se deja el default que es 5
                   env = envi) 
## 
## 
## SAGA Version: 8.1.1 (64 bit)
## 
## library path: C:\PROGRA~1\SAGA\tools\
## library name: ta_channels
## library     : ta_channels
## tool        : Channel Network and Drainage Basins
## identifier  : 5
## author      : O.Conrad (c) 2003
## processors  : 8 [8]
## 
## Loading grid: demsink_fill.sgrd...
## [Channel Network and Drainage Basins] Execution started...
## 
## [Channel Network and Drainage Basins] Parameters:
## 
## Grid System: 31; 2232x 2962y; 751970.713617x 723106.370313y
## Elevation: demsink_fill
## Flow Direction: <not set>
## Flow Connectivity: <not set>
## Strahler Order: Strahler Order
## Drainage Basins: <not set>
## Channels: Channels
## Drainage Basins: Drainage Basins
## Junctions: <not set>
## Threshold: 5
## 
## Flow Direction
## Stream Order
## Junctions
## Drainage Basins
## Vectorising Grid Classes
## 
## [Vectorising Grid Classes] Parameters:
## 
## Grid System: 31; 2232x 2962y; 751970.713617x 723106.370313y
## Grid: Drainage Basins
## Polygons: Polygons
## Class Selection: all classes
## Vectorised class as...: one single (multi-)polygon object
## Keep Vertices on Straight Lines: false
## 
## class identification
## Create index: Drainage Basins
## edge detection
## edge collection
## 
## [Vectorising Grid Classes] execution time: 04s
## Channels
## 
## total execution time: 8000 milliseconds (08s)
## 
## [Channel Network and Drainage Basins] Execution succeeded (08s)
## Saving grid: Ordenes_Strahler_raster...
## Saving shapes: Drenaje_ordenes...

IMPORTANTE: al generarse el shp queda cargado con muchos atributos calculados, ahorra trabajo

Cargar red de drenajes generada

Ordenes=raster("Ordenes_Strahler_raster.sdat")
Ordenes=shapefile("Drenaje_ordenes.shp")
## Warning: [vect] Z coordinates ignored
## M coordinates ignored
Ordenes=spTransform(Ordenes,crs("+init=epsg:3116"))
## Warning: PROJ support is provided by the sf and terra packages among others

4.1. Análisis de la Red de Drenaje

Orden_cuen=max(Ordenes@data$ORDER)

cat("El orden máximo de la cuenca es: ",Orden_cuen)
## El orden máximo de la cuenca es:  7
Long_total=sum(Ordenes@data$LENGTH)/1000 # Km

cat("La longitud total de la corriente es: ",Long_total,"Km.")
## La longitud total de la corriente es:  4611.624 Km.
No_corrientes=max(Ordenes@data$SEGMENT_ID)
No_corrientes <- as.numeric(No_corrientes)

cat("El número total de corrientes es: ",No_corrientes)
## El número total de corrientes es:  4619
Den_corrientes= No_corrientes/A # Unid/Km2

cat("La densidad de las corrientes es: ",Den_corrientes," Unid/Km2")
## La densidad de las corrientes es:  1.013742  Unid/Km2
Den_red_dren=Long_total/A # Km/Km2

cat("La densidad de la red de drenajes es: ",Den_red_dren," Km/Km2")
## La densidad de la red de drenajes es:  1.012124  Km/Km2

4.2. Análisis de la Red de Drenaje Principal —-

Long_dren_prin=sum(Ordenes@data$LENGTH[which(Ordenes@data$ORDER==Orden_cuen)])/1000

cat("La longitud del cauce principal es igual a (por el método de Strahler): ",Long_dren_prin," Km") # Solo muestra la acumulación de las redes, un segmento, no como Horton
## La longitud del cauce principal es igual a (por el método de Strahler):  22.04349  Km
Orden_prin=subset(Ordenes,Ordenes@data$ORDER==Orden_cuen)

red_prin = mask(dem,Orden_prin)
exten = c(Orden_prin@bbox[1,1] ,Orden_prin@bbox[1,2],Orden_prin@bbox[2,1] ,Orden_prin@bbox[2,2])
red_prin = crop(red_prin,exten)

cat("La elevación mínima de la red principal es:",red_prin@data@min,"m.s.n.m.")
## La elevación mínima de la red principal es: 789.2186 m.s.n.m.
cat("La elevación media de la red principal es:",mean(na.omit(red_prin@data@values)),"m.s.n.m.")
## La elevación media de la red principal es: 856.5589 m.s.n.m.
cat("La elevación máxima de la red principal es:",red_prin@data@max,"m.s.n.m.")
## La elevación máxima de la red principal es: 926.726 m.s.n.m.
pend_red_prin=mask(pend,Orden_prin)
exten = c(Orden_prin@bbox[1,1] ,Orden_prin@bbox[1,2],Orden_prin@bbox[2,1] ,Orden_prin@bbox[2,2])
pend_red_prin = crop(pend_red_prin,exten)

s_red_prin=mean(na.omit(pend_red_prin@data@values))/100 # m/m

cat("La pendiente media de la red principal es:",mean(na.omit(pend_red_prin@data@values)),"%")
## La pendiente media de la red principal es: 4.992134 %

4.3. Perfil de Elevacion de la Red de Drenaje Principal

writeRaster(red_prin,"red_prin_dem.tif", overwrite=TRUE)

class(Orden_prin)
## [1] "SpatialLinesDataFrame"
## attr(,"package")
## [1] "sp"
raster::shapefile(Orden_prin, "Orden_prin.shp", overwrite=TRUE)

rsaga.geoprocessor(lib = "io_gdal", 0,
                   param = list(GRIDS ="red_prin_dem",
                                FILES= "red_prin_dem.tif",
                                TRANSFORM= 1),env = envi)
## 
## 
## SAGA Version: 8.1.1 (64 bit)
## 
## library path: C:\PROGRA~1\SAGA\tools\
## library name: io_gdal
## library     : io_gdal
## tool        : Import Raster
## identifier  : 0
## author      : O.Conrad (c) 2007 (A.Ringeler)
## processors  : 8 [8]
## 
## [Import Raster] Execution started...
## 
## [Import Raster] Parameters:
## 
## Grids: No objects
## Files: "red_prin_dem.tif"
## Multiple Bands Output: automatic
## Select from Multiple Bands: 
## Transformation: true
## Resampling: Nearest Neighbour
## Extent: original
## 
## 
## loading: red_prin_dem.tif
## 
## 
## 
## Driver: GTiff
## 
## Bands: 1
## 
## Rows: 350
## 
## Columns: 194
## 
## 
## loading: red_prin_dem
## 
## total execution time: 0 milliseconds (less than 1 millisecond)
## 
## [Import Raster] Execution succeeded (less than 1 millisecond)
## Saving grid: red_prin_dem...
rsaga.geoprocessor(lib = "ta_profiles", 4,
                   param = list(DEM ="red_prin_dem.sgrd",
                                LINES= "Orden_prin.shp",
                                PROFILE="Puntos_datos_perfil"
                                ),env = envi)
## 
## 
## SAGA Version: 8.1.1 (64 bit)
## 
## library path: C:\PROGRA~1\SAGA\tools\
## library name: ta_profiles
## library     : ta_profiles
## tool        : Profiles from Lines
## identifier  : 4
## author      : O.Conrad (c) 2006
## processors  : 8 [8]
## 
## Loading grid: red_prin_dem.sgrd...
## Loading shapes: Orden_prin.shp...
## [Profiles from Lines] Execution started...
## 
## [Profiles from Lines] Parameters:
## 
## Grid System: 31; 350x 194y; 801849.713617x 760616.370313y
## DEM: red_prin_dem
## Values: No objects
## Lines: Orden_prin
## Name: <not set>
## Profiles: Profiles
## Profiles: No objects
## Each Line as new Profile: false
## 
## 
## total execution time: 0 milliseconds (less than 1 millisecond)
## 
## [Profiles from Lines] Execution succeeded (less than 1 millisecond)
## Saving shapes: Puntos_datos_perfil...
Datos_perfil=shapefile("Puntos_datos_perfil.shp")
Datos_perfil=spTransform(Datos_perfil,crs("+init=epsg:9377"))
## Warning: PROJ support is provided by the sf and terra packages among others
Datos_per_orden=Datos_perfil@data[order(Datos_perfil@data$X),]

Datos_elev=as.data.frame(matrix(NA,nrow = length(Datos_per_orden$X),ncol = 3))
Datos_elev[,1]=Datos_per_orden$X
Datos_elev[,2]=Datos_per_orden$Z
Datos_elev[,3]=Datos_per_orden$DIST

jpeg("Perfil_elevacion.jpeg", quality = 75)
par(mar=c(5, 5, 4, 5))
plot(Datos_elev$V1,Datos_elev$V2,type="l", col="#FF8800", lwd = 3,lty = 1,
     xlab="",ylab="ElevaciCon (m.s.n.m.)",axes=F,cex.lab=0.8,
     main=paste("Perfil de elevaciC3n del cauce principal (Orden: ",Orden_cuen," - Met. Strahler)",sep = ""),
    cex.main=1)
axis(2,ylim=Datos_elev$V2,las=1,tcl=-0.2,cex.axis=0.8)
box()
par(new = T)
plot(Datos_elev$V3,Datos_elev$V2,type="n",xlab="Distancia (m.)",ylab="",axes=F,cex.lab=0.8)
axis(1,xlim=Datos_elev$V3,las=1,tcl=-0.2,cex.axis=0.8)
dev.off()
## png 
##   2

5. Tiempos de ConcentraciCón

5.1. Método de Kirpich —-

Tc_K=0.01947*((Long_dren_prin*1000)^0.77)/(s_red_prin^(0.385)) # min
Tc_K=Tc_K/60 # horas

5.2. MC)todo de Temez —-

Tc_T=0.3*((Long_dren_prin/(s_red_prin^0.25))^0.76) # horas

5.3. MC)todo de Giandotti —-

Tc_G=((4*sqrt(A))+(1.5*Long_dren_prin))/(25.3*sqrt(s_red_prin*Long_dren_prin)) # horas

5.4. MC)todo de Passini —-

Tc_P=(0.108*((A*Long_dren_prin)^(1/3)))/(s_red_prin^0.5) # horas

5.5. Tiempo de concentraciC3n promedio —-

Tc=(Tc_K+Tc_T+Tc_G+Tc_P)/4 # horas

6. Exportar todos los datos obtenidos —-

nombres=c("Area","Perimetro","Clas. A","Indice de compacidad","Clas. Kc","Lado 1 Rect. equ.",
          "Lado 2 Rect. equ.", "Elevacion min. de la cuenca","Elevacion med. de la cuenca","Elevacion max. de la cuenca",
          "Pend. media de la cuenca", "Clas. pend.","Orden de la cuenca", "Long. total de las corrientes",
          "No. de corrientes","Dens. corrientes","Dens. red drenaje","Long. cauce principal",
          "Elevacion min. de la RD","Elevacion med. de la RD","Elevacion max. de la RD","Pend. media de la RD",
          "Tiempo de concentracion Kirpich","Tiempo de concentracion Temez","Tiempo de concentracion Giandotti",
          "Tiempo de concentracion Passini","Tiempo de concentracion promedio")

Union_resul=as.data.frame(matrix(NA, nrow = length(nombres), ncol = 4))
names(Union_resul)=c("Relativo","Parametro","Valor","Unidad")
Union_resul[,2]=nombres

Union_resul[1:3,1]="Parametros generales"
Union_resul[4:7,1]="A la forma"
Union_resul[8:12,1]="Al relieve"
Union_resul[13:22,1]="A la red de drenaje"
Union_resul[23:27,1]="Tiempo de concentraciC3n"

Union_resul[1,3]=A
Union_resul[2,3]=P
Union_resul[3,3]=C
Union_resul[4,3]=Kc
Union_resul[5,3]=C_Kc
Union_resul[6,3]=L1
Union_resul[7,3]=L2
Union_resul[8,3]=dem@data@min
Union_resul[9,3]=mean(na.omit(dem@data@values))
Union_resul[10,3]=dem@data@max
Union_resul[11,3]=pend_c
Union_resul[12,3]=pend_clas
Union_resul[13,3]=Orden_cuen
Union_resul[14,3]=Long_total
Union_resul[15,3]=No_corrientes
Union_resul[16,3]=Den_corrientes
Union_resul[17,3]=Den_red_dren
Union_resul[18,3]=Long_dren_prin
Union_resul[19,3]=red_prin@data@min
Union_resul[20,3]=mean(na.omit(red_prin@data@values))
Union_resul[21,3]=red_prin@data@max
Union_resul[22,3]=mean(na.omit(pend_red_prin@data@values))
Union_resul[23,3]=Tc_K
Union_resul[24,3]=Tc_T
Union_resul[25,3]=Tc_G
Union_resul[26,3]=Tc_P
Union_resul[27,3]=Tc

Union_resul[1,4]="Km2"
Union_resul[2,4]="Km"
Union_resul[3,4]="-"
Union_resul[4,4]="Km"
Union_resul[5,4]="-"
Union_resul[6,4]="Km"
Union_resul[7,4]="Km"
Union_resul[8,4]="m.s.n.m."
Union_resul[9,4]="m.s.n.m."
Union_resul[10,4]="m.s.n.m."
Union_resul[11,4]="%"
Union_resul[12,4]="-"
Union_resul[13,4]="Unid."
Union_resul[14,4]="Km"
Union_resul[15,4]="Unid."
Union_resul[16,4]="Unid./Km"
Union_resul[17,4]="Km/Km2"
Union_resul[18,4]="Km"
Union_resul[19,4]="m.s.n.m."
Union_resul[20,4]="m.s.n.m."
Union_resul[21,4]="m.s.n.m."
Union_resul[22,4]="%"
Union_resul[23,4]="Horas"
Union_resul[24,4]="Horas"
Union_resul[25,4]="Horas"
Union_resul[26,4]="Horas"
Union_resul[27,4]="Horas"


write.csv(Union_resul,"Resultados_analisis_morfom2.csv")