code_01_slope_calculation.R

library(raster)

## Loading required package: sp

library(sp)
library(rgdal)

## rgdal: version: 1.2-8, (SVN revision 663)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15
##  Path to GDAL shared files: /usr/share/gdal/2.2
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.2-5

library(rgeos)

## rgeos version: 0.3-27, (SVN revision 560)
##  GEOS runtime version: 3.5.1-CAPI-1.9.1 r4246 
##  Linking to sp version: 1.2-7 
##  Polygon checking: TRUE

library(sf)

## Linking to GEOS 3.5.1, GDAL 2.2.2, proj.4 4.9.2

library(plyr)

setwd("/home/michael/Dropbox/BGU/Chaya_Sud/p_02_slope_between_two_points_on_DEM/")

# Read data
r = raster("DEM_zinZ_apr18.tif")
line = st_read("New_Shapefile.shp")

## Reading layer `New_Shapefile' from data source `/home/michael/Dropbox/BGU/Chaya_Sud/p_02_slope_between_two_points_on_DEM/New_Shapefile.shp' using driver `ESRI Shapefile'
## Simple feature collection with 73 features and 2 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -86.2232 ymin: 0.5826478 xmax: 51.47417 ymax: 134.2049
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=31.73439361111111 +lon_0=35.20451694444445 +k=1.0000067 +x_0=219529.584 +y_0=626907.39 +ellps=GRS80 +units=m +no_defs

# Set line ID
line$id = 1:nrow(line)

# Plot 1
plot(r)
plot(st_geometry(line), add = TRUE)

# Find highest point
pnt = st_cast(line, "POINT")

## Warning in st_cast.sf(line, "POINT"): repeating attributes for all sub-
## geometries for which they may not be constant

pnt$seq = ave(pnt$id, pnt$id, FUN = seq_along)
flag = ddply(
  st_set_geometry(pnt, NULL), 
  "id",
  summarize,
  flag = seq == min(seq) | seq == max(seq)
)$flag
pnt = pnt[flag, ]
type = ddply(
  st_set_geometry(pnt, NULL), 
  "id",
  summarize,
  type = ifelse(seq == min(seq), "first", "last")
)
pnt$type = type$type

# Plot 2
plot(st_geometry(line))
plot(st_geometry(pnt[pnt$type == "first", ]), col = "red", add = TRUE)
plot(st_geometry(pnt[pnt$type == "last", ]), col = "blue", add = TRUE)

# Find second point
line_sp = as(line, "Spatial")
second = list()
for(i in unique(line_sp$id)) {
  second[[i]] = gInterpolate(line_sp[i, ], 1)
}
second = do.call(rbind, second)
second = st_as_sf(second)
second$id = line$id
first = pnt[pnt$type == "first", ]

# Plot 2
plot(st_geometry(line))
plot(st_geometry(first), col = "red", add = TRUE)
plot(st_geometry(second), col = "blue", add = TRUE)

# Extract data
dat = data.frame(id = first$id, elev1 = NA, elev2 = NA, dist = NA, stringsAsFactors = FALSE)
for(i in unique(line_sp$id)) {
  dat$elev1[i] = extract(r, first[i, ])
  dat$elev2[i] = extract(r, second[i, ])
  dat$dist[i] = st_distance(first[i, ], second[i, ]) %>% as.numeric
}

# Calculate angle
dat$slope = atan((dat$elev1 - dat$elev2) / dat$dist)
dat$slope = (dat$slope * 180) / (pi)

# Save layer
dat = dplyr::left_join(first, dat, "id")
st_write(dat, "slope.shp", delete_dsn = TRUE)

## Deleting source `/home/michael/Dropbox/BGU/Chaya_Sud/p_02_slope_between_two_points_on_DEM/slope.shp' using driver `ESRI Shapefile'
## Writing layer `slope' to data source `/home/michael/Dropbox/BGU/Chaya_Sud/p_02_slope_between_two_points_on_DEM/slope.shp' using driver `ESRI Shapefile'

## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 6: Normalized/laundered field
## name: 'id' to 'id_1'

## features:       73
## fields:         9
## geometry type:  Point

# Plot 4
plot(st_geometry(line))
text(as(dat, "Spatial"), round(dat$slope, 2))