# NOT RUN {
m <- datasets::volcano
m <- m[nrow(m):1, ncol(m):1]
x <- seq(from = 2667405, length.out = ncol(m), by = 10)
y <- seq(from = 6478705, length.out = nrow(m), by = 10)
r1 <- raster::raster(m, xmn = min(x), xmx = max(x), ymn = min(y), ymx = max(y),
crs = "+init=epsg:27200")
r2 <- min(r1[]) - r1 / 10
r3 <- r1 - r2
rs <- raster::stack(r1, r2, r3)
names(rs) <- c("r1", "r2", "r3")
xy <- rbind(c(2667508, 6479501), c(2667803, 6479214), c(2667508, 6478749))
transect <- sp::SpatialLines(list(sp::Lines(list(sp::Line(xy)), ID = "Transect")),
proj4string = raster::crs(rs))
PlotMap(r1, pal = terrain.colors, scale.loc = "top", arrow.loc = "topright",
shade = list(alpha = 0.3), contour.lines = list(col = "#1F1F1FA6"))
lines(transect)
raster::text(as(transect, "SpatialPoints"), labels = c("A", "BEND", "A'"),
cex = 0.7, pos = c(3, 4, 1), offset = 0.1, font = 4)
dev.new()
explanation <- "Vertical thickness between layers, in meters."
PlotCrossSection(transect, rs, geo.lays = c("r1", "r2"), val.lays = "r3",
ylab = "Elevation", asp = 5, unit = "METERS",
explanation = explanation, scale.loc = "bottomright")
val <- PlotCrossSection(transect, rs, geo.lays = c("r1", "r2"), val.lays = "r3",
ylab = "Elevation", asp = 5, unit = "METERS",
explanation = explanation, file = "Rplots.png")
print(val)
file.remove("Rplots.png")
graphics.off()
# }
Run the code above in your browser using DataLab