GeomorfometrÃa
Daniela Ballari
Octubre, 2015
MaestrÃa en EcohidrologÃa
Universidad de Cuenca
http://www.ucuenca.edu.ec/maestria-ecohidrologia/
Daniela Ballari (PhD)
daniela.ballari@ucuenca.edu.ec
http://www.researchgate.net/profile/Daniela_Ballari/publications
Objetivo de la práctica. Explorar el uso de la libreria RSAGA para analisis geomorfométrico. Este material fue creado utilizando: - https://cran.r-project.org/web/packages/RSAGA/vignettes/RSAGA-landslides.pdf
- Vern Cimmery (2010) User Guide for SAGA (version 2.0.5) Volume 2. Capitulos 4 y 5
- http://www.headwateranalytics.com/blog/delineate-a-watershed-with-saga-gis-and-r4
- http://www.headwateranalytics.com/blog/delineating-a-watershed-with-saga-gis-and-r-part-2
Temario
A- Cargar librerias y datos
B- Explorar RSAGA
C- Pendiente
D-Delimitación de cuencas hidrográficas
A- Cargar librerias y datos
library(raster)## Loading required package: splibrary(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: plyrlibrary(rgdal)## rgdal: version: 1.0-7, (SVN revision 559)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
## Path to GDAL shared files: C:/Users/dani/Documents/R/win-library/3.2/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: C:/Users/dani/Documents/R/win-library/3.2/rgdal/proj
## Linking to sp version: 1.2-0library(latticeExtra)## Loading required package: RColorBrewer
## Loading required package: latticesetwd("C:/R_ecohidrologia/geomorfometria")
dem <- raster("s03_w079_3arc_v2.tif")
projection(dem)## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"dem_utm <- projectRaster(dem,crs=("+init=epsg:32717"))
zona_estudio<-readOGR("C:/R_ecohidrologia/SRS","prov_zona_estudio")## OGR data source with driver: ESRI Shapefile
## Source: "C:/R_ecohidrologia/SRS", layer: "prov_zona_estudio"
## with 7 features
## It has 11 fieldsspplot(dem_utm, col.regions=terrain.colors(255))
spplot(zona_estudio[1], fill="transparent", col="black", under = FALSE)+as.layer(spplot(dem_utm, col.regions=terrain.colors(255)))
writeRaster(dem_utm, filename="dem_utm.sgrd", format="SAGA", overwrite=TRUE)## class : RasterLayer
## dimensions : 1213, 1214, 1472582 (nrow, ncol, ncell)
## resolution : 92.7, 92.7 (x, y)
## extent : 721782.3, 834320.1, 9667472, 9779917 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## data source : C:\R_ecohidrologia\geomorfometria\dem_utm.sgrd
## names : s03_w079_3arc_v2
## values : 311.812, 5284.223 (min, max)B- Explorar RSAGA
rsaga.env()## $workspace
## [1] "."
##
## $cmd
## [1] "saga_cmd.exe"
##
## $path
## [1] "C:/Progra~1/SAGA-GIS"
##
## $modules
## [1] "C:/Progra~1/SAGA-GIS/modules"
##
## $version
## [1] "2.1.0"
##
## $cores
## [1] NA
##
## $parallel
## [1] FALSE
##
## $lib.prefix
## [1] ""rsaga.get.version()## [1] "2.1.0"# rsaga.get.modules()
#funciones en R
rsaga.hillshade("dem_utm.sgrd","hillshade")##
##
## library path: C:\Progra~1\SAGA-GIS\modules\ta_lighting.dll
## library name: Terrain Analysis - Lighting, Visibility
## module name : Analytical Hillshading
## author : O.Conrad, V.Wichmann (c) 2003-2013
##
## Load grid: dem_utm.sgrd...
##
##
## Parameters
##
##
## Grid system: 92.7; 1214x 1213y; 721828.605297x 9667518.100606y
## Elevation: dem_utm
## Analytical Hillshading: Analytical Hillshading
## Shading Method: Standard
## Azimuth [Degree]: 315.000000
## Declination [Degree]: 45.000000
## Exaggeration: 4.000000
## Number of Directions: 8
## Search Radius: 100.000000
##
## Save grid: hillshade.sgrd...hillshade <- raster("hillshade.sgrd")
spplot(hillshade,col.regions=rev(gray(seq(0,1,.01))))
#funciones en saga
rsaga.get.usage("ta_lighting",0)## library path: C:\Progra~1\SAGA-GIS\modules\ta_lighting.dll
## library name: Terrain Analysis - Lighting, Visibility
## module name : Analytical Hillshadingrsaga.geoprocessor("ta_lighting",0,list(ELEVATION="dem_utm.sgrd",SHADE="hillshade2",EXAGGERATION=2))##
##
## library path: C:\Progra~1\SAGA-GIS\modules\ta_lighting.dll
## library name: Terrain Analysis - Lighting, Visibility
## module name : Analytical Hillshading
## author : O.Conrad, V.Wichmann (c) 2003-2013
##
## Load grid: dem_utm.sgrd...
##
##
## Parameters
##
##
## Grid system: 92.7; 1214x 1213y; 721828.605297x 9667518.100606y
## Elevation: dem_utm
## Analytical Hillshading: Analytical Hillshading
## Shading Method: Standard
## Azimuth [Degree]: 315.000000
## Declination [Degree]: 45.000000
## Exaggeration: 2.000000
## Number of Directions: 8
## Search Radius: 100.000000
##
## Save grid: hillshade2...hillshade2 <- raster("hillshade2.sgrd")
spplot(hillshade2,col.regions=rev(gray(seq(0,1,.01))))
C-Pendiente
#Preprocessing/Sink Removal
rsaga.sink.removal("dem_utm.sgrd",,"demsink",method=1)##
##
## library path: C:\Progra~1\SAGA-GIS\modules\ta_preprocessor.dll
## library name: Terrain Analysis - Preprocessing
## module name : Sink Removal
## author : O. Conrad (c) 2001
##
## Load grid: dem_utm.sgrd...
##
##
## Parameters
##
##
## Grid system: 92.7; 1214x 1213y; 721828.605297x 9667518.100606y
## DEM: dem_utm
## Sink Route: <not set>
## Preprocessed DEM: Preprocessed DEM
## Method: Fill Sinks
## Threshold: no
## Threshold Height: 100.000000
##
## Find Pits
## Create index: dem_utm [no sinks]
## Find Outlets
## Routing
## Finalize
## number of processed sinks: 8778
## Initializing direction matrix...
## I'm fillin'...
## Save grid: demsink...demsink <- raster("demsink.sdat")
spplot(demsink,col.regions=terrain.colors(255))
#Calculo de pendiente con función rsaga.slope. rsaga.slope(in.dem, out.slope, method = "maxslope", env = rsaga.env(), ...)
rsaga.slope("demsink.sgrd",out.slope="slope",0)##
##
## library path: C:\Progra~1\SAGA-GIS\modules\ta_morphometry.dll
## library name: Terrain Analysis - Morphometry
## module name : Slope, Aspect, Curvature
## author : O.Conrad (c) 2001
##
## Load grid: demsink.sgrd...
##
##
## Parameters
##
##
## Grid system: 92.7; 1214x 1213y; 721828.605297x 9667518.100606y
## Elevation: demsink
## Slope: Slope
## Aspect: Aspect
## Curvature: <not set>
## Plan Curvature: <not set>
## Profile Curvature: <not set>
## Method: Maximum Slope (Travis et al. 1975)
##
## Save grid: slope...
## Save grid: C:\Users\dani\AppData\Local\Temp\Rtmp0K73ov\file292049f512ef...slope <- raster("slope.sgrd")
#convertir de radianes a grados
spplot(slope, zlim=c(0,40))
rsaga.get.modules("ta_morphometry")## $ta_morphometry
## code name
## 1 0 Slope, Aspect, Curvature
## 2 1 Convergence Index
## 3 2 Convergence Index (Search Radius)
## 4 3 Surface Specific Points
## 5 4 Curvature Classification
## 6 5 Hypsometry
## 7 6 Real Surface Area
## 8 7 Morphometric Protection Index
## 9 8 Multiresolution Index of Valley Bottom Flatness (MRVBF)
## 10 9 Downslope Distance Gradient
## 11 10 Mass Balance Index
## 12 11 Effective Air Flow Heights
## 13 12 Diurnal Anisotropic Heating
## 14 13 Land Surface Temperature
## 15 14 Relative Heights and Slope Positions
## 16 15 Wind Effect (Windward / Leeward Index)
## 17 16 Terrain Ruggedness Index (TRI)
## 18 17 Vector Ruggedness Measure (VRM)
## 19 18 Topographic Position Index (TPI)
## 20 19 TPI Based Landform Classification
## 21 20 Terrain Surface Texture
## 22 21 Terrain Surface Convexity
## 23 22 Terrain Surface Classification
## 24 23 Morphometric Features
## 25 24 Valley and Ridge Detection (Top Hat Approach)
## 26 25 Fuzzy Landform Element Classification
## interactive
## 1 FALSE
## 2 FALSE
## 3 FALSE
## 4 FALSE
## 5 FALSE
## 6 FALSE
## 7 FALSE
## 8 FALSE
## 9 FALSE
## 10 FALSE
## 11 FALSE
## 12 FALSE
## 13 FALSE
## 14 FALSE
## 15 FALSE
## 16 FALSE
## 17 FALSE
## 18 FALSE
## 19 FALSE
## 20 FALSE
## 21 FALSE
## 22 FALSE
## 23 FALSE
## 24 FALSE
## 25 FALSE
## 26 FALSErsaga.get.usage("ta_morphometry",0)## library path: C:\Progra~1\SAGA-GIS\modules\ta_morphometry.dll
## library name: Terrain Analysis - Morphometry
## module name : Slope, Aspect, Curvaturersaga.geoprocessor("ta_morphometry",0,list(ELEVATION="demsink.sgrd",SLOPE="slope2",METHOD =0))##
##
## library path: C:\Progra~1\SAGA-GIS\modules\ta_morphometry.dll
## library name: Terrain Analysis - Morphometry
## module name : Slope, Aspect, Curvature
## author : O.Conrad (c) 2001
##
## Load grid: demsink.sgrd...
##
##
## Parameters
##
##
## Grid system: 92.7; 1214x 1213y; 721828.605297x 9667518.100606y
## Elevation: demsink
## Slope: Slope
## Aspect: Aspect
## Curvature: <not set>
## Plan Curvature: <not set>
## Profile Curvature: <not set>
## Method: Maximum Slope (Travis et al. 1975)
##
## Save grid: slope2...slope2 <- raster("slope2.sgrd")
spplot(slope2, zlim=c(0,40))
D-Delimitación de cuencas hidrográficas
# 1- Area de Captación - Catchment Area
# Para una celda, el area de captación es el area de todas las celdas aguas arriba que proveerán con flujo de agua a dicha celda. Algoritmos de ruteo de flujo: D8, Rho8,FD8, Dinf, KRA, DEMON
rsaga.get.modules("ta_hydrology")## $ta_hydrology
## code name interactive
## 1 0 Catchment Area (Parallel) FALSE
## 2 1 Catchment Area (Recursive) FALSE
## 3 2 Catchment Area (Flow Tracing) FALSE
## 4 4 Upslope Area FALSE
## 5 6 Flow Path Length FALSE
## 6 7 Slope Length FALSE
## 7 10 Cell Balance FALSE
## 8 13 Edge Contamination FALSE
## 9 15 SAGA Wetness Index FALSE
## 10 16 Lake Flood FALSE
## 11 18 Catchment Area (Mass-Flux Method) FALSE
## 12 19 Flow Width and Specific Catchment Area FALSE
## 13 20 Topographic Wetness Index (TWI) FALSE
## 14 21 Stream Power Index FALSE
## 15 22 LS Factor FALSE
## 16 23 Melton Ruggedness Number FALSE
## 17 24 TCI Low FALSErsaga.get.usage("ta_hydrology",0) ; rsaga.get.usage("ta_hydrology","Catchment Area (Parallel)")## library path: C:\Progra~1\SAGA-GIS\modules\ta_hydrology.dll
## library name: Terrain Analysis - Hydrology
## module name : Catchment Area (Parallel)## library path: C:\Progra~1\SAGA-GIS\modules\ta_hydrology.dll
## library name: Terrain Analysis - Hydrology
## module name : Catchment Area (Parallel)rsaga.parallel.processing(in.dem = "demsink", out.carea = "carea", out.cslope = "cslope", out.caspect="caspect", method=4)#Multiple Flow Direction##
##
## library path: C:\Progra~1\SAGA-GIS\modules\ta_hydrology.dll
## library name: Terrain Analysis - Hydrology
## module name : Catchment Area (Parallel)
## author : O.Conrad (c) 2001-2010, T.Grabs portions (c) 2010
##
## Load grid: demsink.sgrd...
##
##
## Parameters
##
##
## Grid system: 92.7; 1214x 1213y; 721828.605297x 9667518.100606y
## Elevation: demsink
## Sink Routes: <not set>
## Weight: <not set>
## Material: <not set>
## Target: <not set>
## Catchment Area: Catchment Area
## Catchment Height: <not set>
## Catchment Slope: Catchment Slope
## Total accumulated Material: <not set>
## Accumulated Material from _left_ side: <not set>
## Accumulated Material from _right_ side: <not set>
## Step: 1
## Catchment Aspect: Catchment Aspect
## Flow Path Length: <not set>
## Method: Multiple Flow Direction
## Linear Flow: no
## Linear Flow Threshold: 500.000000
## Linear Flow Threshold Grid: <not set>
## Channel Direction: <not set>
## Convergence: 1.100000
##
## Create index: demsink
## Save grid: carea...
## Save grid: cslope...
## Save grid: caspect...carea <- raster("carea.sdat")
spplot(carea) #0 sin contribución y 1 contribución total
# cslope está en radianes, debe convertirse a grados
rad_a_grado = round(180/pi, 4)
formula = paste(rad_a_grado, "*a", sep = "")
rsaga.grid.calculus("cslope", "cslopegrados", formula)##
##
## library path: C:\Progra~1\SAGA-GIS\modules\grid_calculus.dll
## library name: Grid - Calculus
## module name : Grid Calculator
## author : Copyrights (c) 2003 by Andre Ringeler
##
## Load grid: cslope.sgrd...
##
##
## Parameters
##
##
## Grid system: 92.7; 1214x 1213y; 721828.605297x 9667518.100606y
## Grids: 1 object (cslope)
## Grids from different Systems: No objects
## Result: Result
## Formula: 57.2958*a
## Name: Calculation
## Take Formula: no
## Use NoData: no
## Data Type: 4 byte floating point number
##
## Save grid: cslopegrados...slope <- raster("cslopegrados.sdat")
spplot(slope) 
aspect=raster("caspect.sdat")
spplot(aspect) 
#2- Redes de drenaje - Channel Networks
# Identifca cuales son cauces y cuales no. Cauces: flujo encauzado, no cause: flujo en laderas
# Area de Captación es comúnmente usado como initation grid
rsaga.get.modules("ta_channels")## $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 FALSErsaga.get.usage("ta_channels","Channel Network")## library path: C:\Progra~1\SAGA-GIS\modules\ta_channels.dll
## library name: Terrain Analysis - Channels
## module name : Channel Networkrsaga.geoprocessor("ta_channels",0,list(ELEVATION="demsink.sgrd",CHNLNTWRK="channel_net.sdat",CHNLROUTE="route.sdat",SHAPES='channel_net_shp.shp',INIT_GRID="carea.sgrd",INIT_METHOD= 2,INIT_VALUE=5000000))##
##
## library path: C:\Progra~1\SAGA-GIS\modules\ta_channels.dll
## library name: Terrain Analysis - Channels
## module name : Channel Network
## author : (c) 2001 by O.Conrad
##
## Load grid: demsink.sgrd...
## Load grid: carea.sgrd...
##
##
## Parameters
##
##
## Grid system: 92.7; 1214x 1213y; 721828.605297x 9667518.100606y
## Elevation: demsink
## Flow Direction: <not set>
## Channel Network: Channel Network
## Channel Direction: Channel Direction
## Channel Network: Channel Network
## Initiation Grid: carea
## Initiation Type: Greater than
## Initiation Threshold: 5000000.000000
## Divergence: <not set>
## Tracing: Max. Divergence: 5
## Tracing: Weight: <not set>
## Min. Segment Length: 10
##
## Channel Network: Pass 1
## Channel Network: Pass 2
## Channel Network: Pass 3
## Create index: demsink
## Channel Network: Pass 4
## Channel Network: Pass 5
## Channel Network: Pass 6
## Save grid: channel_net.sdat...
## Save grid: route.sdat...
## Save shapes: channel_net_shp.shp...channel <- readOGR(".","channel_net_shp")## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "channel_net_shp"
## with 1011 features
## It has 3 fieldsspplot(channel, "Order")
#3- Delimitación de cuencas -Watershed Basins
rsaga.get.usage("ta_channels","Watershed Basins")## library path: C:\Progra~1\SAGA-GIS\modules\ta_channels.dll
## library name: Terrain Analysis - Channels
## module name : Watershed Basinsrsaga.geoprocessor("ta_channels",1,list(ELEVATION="demsink.sgrd", CHANNELS="channel_net.sgrd", BASINS="basins.sgrd",MINSIZE=100))##
##
## library path: C:\Progra~1\SAGA-GIS\modules\ta_channels.dll
## library name: Terrain Analysis - Channels
## module name : Watershed Basins
## author : (c) 2001 by O.Conrad
##
## Load grid: demsink.sgrd...
## Load grid: channel_net.sgrd...
##
##
## Parameters
##
##
## Grid system: 92.7; 1214x 1213y; 721828.605297x 9667518.100606y
## Elevation: demsink
## Channel Network: channel_net
## Sink Route: <not set>
## Watershed Basins: Watershed Basins
## Min. Size: 100
##
## Create index: demsink
## Save grid: basins.sgrd...basin <- raster("basins.sdat")
spplot(basin)
#Upslope Area
rsaga.get.usage('ta_hydrology', 4)## library path: C:\Progra~1\SAGA-GIS\modules\ta_hydrology.dll
## library name: Terrain Analysis - Hydrology
## module name : Upslope Arearsaga.geoprocessor(lib = 'ta_hydrology', 4,
param = list(TARGET_PT_X =808160.095,
TARGET_PT_Y =9701160.897,
ELEVATION = 'demsink.sgrd',
AREA = 'upslope.sgrd',
METHOD = 2))##
##
## library path: C:\Progra~1\SAGA-GIS\modules\ta_hydrology.dll
## library name: Terrain Analysis - Hydrology
## module name : Upslope Area
## author : (c) 2001 by O.Conrad
##
## Load grid: demsink.sgrd...
##
##
## Parameters
##
##
## Grid system: 92.7; 1214x 1213y; 721828.605297x 9667518.100606y
## Target Area: <not set>
## Target X coordinate: 808160.095000
## Target Y coordinate: 9701160.897000
## Elevation: demsink
## Sink Routes: <not set>
## Upslope Area: Upslope Area
## Method: Multiple Flow Direction
## Convergence: 1.100000
##
## Create index: demsink
## Save grid: upslope.sgrd...upslope <- raster("upslope.sdat")
spplot(upslope)
rsaga.geoprocessor(lib = 'shapes_grid', 6,
param = list(GRID = 'upslope.sgrd',
POLYGONS = 'boundary.shp',
CLASS_ALL = 0,
CLASS_ID = 100,
SPLIT = 0))##
##
## library path: C:\Progra~1\SAGA-GIS\modules\shapes_grid.dll
## library name: Shapes - Grid
## module name : Vectorising Grid Classes
## author : (c) 2008 by O.Conrad
##
## Load grid: upslope.sgrd...
##
##
## Parameters
##
##
## Grid system: 92.7; 1214x 1213y; 721828.605297x 9667518.100606y
## Grid: upslope
## Polygons: Polygons
## Class Selection: one single class specified by class identifier
## Class Identifier: 100.000000
## Vectorised class as...: one single (multi-)polygon object
##
## class identification
## edge detection
## edge collection
## Save shapes: boundary.shp...basin <- readOGR(".","boundary")## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "boundary"
## with 1 features
## It has 3 fields# plot it onto the DEM
plot(dem_utm)
plot(basin, add = T) 
Has llegado al final del material!