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 Pulgarin, estudiante de la Maestría en Hidrosistemas en la PUJ.

Iniciamos Cargando las librerías

library("raster") #Para que se puedan leer archivos raster
## Loading required package: sp
library("RSAGA") #Para que R articule con SAGA-GIS
## 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
library("terra")
## terra 1.5.21
library("rgdal")
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.3.2, released 2021/09/01
## Path to GDAL shared files: C:/Users/maria.arenasb/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/maria.arenasb/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-7
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## 
## Attaching package: 'rgdal'
## The following object is masked from 'package:terra':
## 
##     project
library("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
library("gstat")
library("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
## Loading required package: xts
## 
## Attaching package: 'hydroTSM'
## The following object is masked from 'package:terra':
## 
##     extract
## The following object is masked from 'package:raster':
## 
##     extract

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

dem=raster("dem_on.tif")
projection(dem)
## [1] "+proj=tmerc +lat_0=4 +lon_0=-73 +k=0.9992 +x_0=5000000 +y_0=2000000 +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_on.shp")
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum Marco_Geocentrico_Nacional_de_Referencia in Proj4
## definition: +proj=tmerc +lat_0=4 +lon_0=-73 +k=0.9992 +x_0=5000000 +y_0=2000000
## +ellps=GRS80 +units=m +no_defs
projection(cuenca)
## [1] "+proj=tmerc +lat_0=4 +lon_0=-73 +k=0.9992 +x_0=5000000 +y_0=2000000 +ellps=GRS80 +units=m +no_defs"
class(cuenca)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
cuenca=spTransform(cuenca,crs("+init=epsg:9377"))
projection(dem) = crs("+init=epsg:9377")

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)
## Warning: [vect] argument 'crs' should be a character value
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 4892.359 Km^2 y el perímetro es igual a:  433.1413 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 oval oblonga a 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: 722 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: 2490.168 m.s.n.m.
cat("La elevación máxima es:",dem@data@max,"m.s.n.m.")
## La elevación máxima es: 5388 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
hypsometric(sgdf, band=1, main="Curva Hipsométrica",
            xlab="Crea relativa por encima de la elevación (a/A)",
            ylab="Elevación relativa (h/H)", col="blue")
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.06585 %
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)
## Warning in rsaga.geoprocessor(lib = NULL, prefix = "--version", show.output.on.console = FALSE, : This RSAGA version has been tested with SAGA GIS versions 2.3.1 - 6.3.0.
## You seem to be using SAGA GIS 8.1.1, which may cause problems due to
## changes in names and definitions of SAGA module arguments, etc.
## [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)
## Warning in rsaga.geoprocessor(lib = "io_gdal", 0, param = list(GRIDS = "dem_sgrd", : This RSAGA version has been tested with SAGA GIS versions 2.3.1 - 6.3.0.
## You seem to be using SAGA GIS 8.1.1, which may cause problems due to
## changes in names and definitions of SAGA module arguments, etc.
## Warning in rsaga.geoprocessor(lib, module = NULL, env = env, intern = TRUE, : This RSAGA version has been tested with SAGA GIS versions 2.3.1 - 6.3.0.
## You seem to be using SAGA GIS 8.1.1, which may cause problems due to
## changes in names and definitions of SAGA module arguments, etc.
## 
## 
## 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  : 6 [6]
## 
## [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: 2483
## 
## Columns: 3400
## 
## 
## 
## Transformation:
## 
##   x' = 4624217.175893 + x * 31.100171 + y * 0.000000
## 
##   y' = 1895116.516419 + x * 0.000000 + y * -31.100171
## 
## 
## loading: dem_cuenca
## translation: dem_cuenca
## 
## warning: top-to-bottom and left-to-right cell sizes differ.
##  Difference: 0.000000
## 
##  using cellsize: 31.100171
## 
## 
## total execution time: 1000 milliseconds (01s)
## 
## [Import Raster] Execution succeeded (01s)
## 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)
## Warning in rsaga.geoprocessor(lib = "ta_preprocessor", 2, param = list(DEM = "dem_sgrd.sgrd", : This RSAGA version has been tested with SAGA GIS versions 2.3.1 - 6.3.0.
## You seem to be using SAGA GIS 8.1.1, which may cause problems due to
## changes in names and definitions of SAGA module arguments, etc.

## Warning in rsaga.geoprocessor(lib = "ta_preprocessor", 2, param = list(DEM = "dem_sgrd.sgrd", : This RSAGA version has been tested with SAGA GIS versions 2.3.1 - 6.3.0.
## You seem to be using SAGA GIS 8.1.1, which may cause problems due to
## changes in names and definitions of SAGA module arguments, etc.
## 
## 
## 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  : 6 [6]
## 
## Loading grid: dem_sgrd.sgrd...
## [Sink Removal] Execution started...
## 
## [Sink Removal] Parameters:
## 
## Grid System: 31.100171; 2483x 3400y; 4624232.725979x 1789391.484738y
## 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: 30091
## Initializing direction matrix...
## I'm fillin'...
## 
## total execution time: 25000 milliseconds (25s)
## 
## [Sink Removal] Execution succeeded (25s)
## Saving grid: demsink_fill...
# Calcular red de drenajes, lo hace por Strahler
# rsaga.get.modules("ta_channels",env = envi)
# rsaga.get.usage('ta_channels', 5,env = envi)
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) 
## Warning in rsaga.geoprocessor(lib = "ta_channels", 5, param = list(DEM = "demsink_fill.sgrd", : This RSAGA version has been tested with SAGA GIS versions 2.3.1 - 6.3.0.
## You seem to be using SAGA GIS 8.1.1, which may cause problems due to
## changes in names and definitions of SAGA module arguments, etc.

## Warning in rsaga.geoprocessor(lib = "ta_channels", 5, param = list(DEM = "demsink_fill.sgrd", : This RSAGA version has been tested with SAGA GIS versions 2.3.1 - 6.3.0.
## You seem to be using SAGA GIS 8.1.1, which may cause problems due to
## changes in names and definitions of SAGA module arguments, etc.
## 
## 
## 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  : 6 [6]
## 
## Loading grid: demsink_fill.sgrd...
## [Channel Network and Drainage Basins] Execution started...
## 
## [Channel Network and Drainage Basins] Parameters:
## 
## Grid System: 31.100171; 2483x 3400y; 4624232.725979x 1789391.484738y
## Elevation: demsink_fill
## Flow Direction: <not set>
## Flow Connectivity: <not set>
## Strahler Order: <not set>
## 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.100171; 2483x 3400y; 4624232.725979x 1789391.484738y
## 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: 06s
## Channels
## 
## total execution time: 10000 milliseconds (10s)
## 
## [Channel Network and Drainage Basins] Execution succeeded (10s)
## 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 in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum Marco_Geocentrico_Nacional_de_Referencia in Proj4
## definition: +proj=tmerc +lat_0=4 +lon_0=-73 +k=0.9992 +x_0=5000000 +y_0=2000000
## +ellps=GRS80 +units=m +no_defs
## Warning in rgdal::readOGR(dirname(x), fn, stringsAsFactors = stringsAsFactors, :
## Z-dimension discarded
Ordenes=spTransform(Ordenes,crs("+init=epsg:9377"))

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:  5655.505 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:  999
Den_corrientes= No_corrientes/A # Unid/Km2

cat("La densidad de las corrientes es: ",Den_corrientes," Unid/Km2")
## La densidad de las corrientes es:  0.204196  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.155987  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):  35.51789  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: 722 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: 820.7112 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: 929 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: 3.065191 %

4.3. Perfil de Elevacion de la Red de Drenaje Principal

writeRaster(red_prin,"red_prin_dem.tif", overwrite=TRUE)
writeOGR(Orden_prin, ".","Orden_prin",driver="ESRI Shapefile",overwrite=TRUE)

rsaga.geoprocessor(lib = "io_gdal", 0,
                   param = list(GRIDS ="red_prin_dem",
                                FILES= "red_prin_dem.tif",
                                TRANSFORM= 1),env = envi)
## Warning in rsaga.geoprocessor(lib = "io_gdal", 0, param = list(GRIDS = "red_prin_dem", : This RSAGA version has been tested with SAGA GIS versions 2.3.1 - 6.3.0.
## You seem to be using SAGA GIS 8.1.1, which may cause problems due to
## changes in names and definitions of SAGA module arguments, etc.
## Warning in rsaga.geoprocessor(lib, module = NULL, env = env, intern = TRUE, : This RSAGA version has been tested with SAGA GIS versions 2.3.1 - 6.3.0.
## You seem to be using SAGA GIS 8.1.1, which may cause problems due to
## changes in names and definitions of SAGA module arguments, etc.
## 
## 
## 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  : 6 [6]
## 
## [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: 620
## 
## Columns: 190
## 
## 
## 
## Transformation:
## 
##   x' = 4682125.694496 + x * 31.100171 + y * 0.000000
## 
##   y' = 1832885.074033 + x * 0.000000 + y * -31.100171
## 
## 
## loading: red_prin_dem
## translation: red_prin_dem
## 
## warning: top-to-bottom and left-to-right cell sizes differ.
##  Difference: 0.000000
## 
##  using cellsize: 31.100171
## 
## 
## 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)
## Warning in rsaga.geoprocessor(lib = "ta_profiles", 4, param = list(DEM = "red_prin_dem.sgrd", : This RSAGA version has been tested with SAGA GIS versions 2.3.1 - 6.3.0.
## You seem to be using SAGA GIS 8.1.1, which may cause problems due to
## changes in names and definitions of SAGA module arguments, etc.

## Warning in rsaga.geoprocessor(lib = "ta_profiles", 4, param = list(DEM = "red_prin_dem.sgrd", : This RSAGA version has been tested with SAGA GIS versions 2.3.1 - 6.3.0.
## You seem to be using SAGA GIS 8.1.1, which may cause problems due to
## changes in names and definitions of SAGA module arguments, etc.
## 
## 
## 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  : 6 [6]
## 
## Loading grid: red_prin_dem.sgrd...
## Loading shapes: Orden_prin.shp...
## [Profiles from Lines] Execution started...
## 
## [Profiles from Lines] Parameters:
## 
## Grid System: 31.100171; 620x 190y; 4682141.244582x 1826991.591608y
## 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")
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum Marco_Geocentrico_Nacional_de_Referencia in Proj4
## definition: +proj=tmerc +lat_0=4 +lon_0=-73 +k=0.9992 +x_0=5000000 +y_0=2000000
## +ellps=GRS80 +units=m +no_defs
Datos_perfil=spTransform(Datos_perfil,crs("+init=epsg:9377"))

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="ElevaciC3n (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")