library(raster)
library(sp)
library(rgdal)
#######################################
setwd('/Users/COLSON/Documents/Estaciones_CNA')
getwd()
list.files()
#load raster in an R object called 'DEM'
# este es el raster que cortamos usando el poligono de
# la cuenca.
DEM <- raster("cuenca_MED.tif")
DEM
# damos un vistazo
plot(DEM)
cellStats(DEM,'min') # 396 metros altitud minima
cellStats(DEM,'max') # 1522 metros altitud maxima
###############################
# extraemos todas las altitudes
###############################
e <- extent(DEM)
valores=extract(DEM, e)
valores<- na.omit(valores)
valores[1]
hist(valores)
percent=(quantile(valores, prob = seq(0, 1, length = 11), type = 5))
x=c(0,10,20,30,40,50,60,70,80,90,100)
plot(x,(percent), pch=19, xaxt="n", xlab="Deciles", ylab="Altitud (metros)",
col="red")
axis(1, at = seq(0, 100, by = 10), las=2)
lines(x,(percent))
text(x+2, (percent)+8,labels=(percent), cex= 0.8, offset = 10)
396486047/length(valores) # número de pixeles entre el área
(length(which(valores>556))*209.8426)/10000 #3959.625 hectáreas
(length(which(valores>656))*209.8426)/10000 #1454.356 hectáreas
(length(which(valores>756))*209.8426)/10000 #1177.951 hectáreas
(length(which(valores>856))*209.8426)/10000 #1006.006 hectáreas
(length(which(valores>956))*209.8426)/10000 #852.4226 hectáreas
(length(which(valores>1056))*209.8426)/10000 #733.0012 hectáreas
(length(which(valores>1156))*209.8426)/10000 #601.5977 hectáreas
(length(which(valores>1256))*209.8426)/10000 #430.8698 hectáreas
(length(which(valores>1356))*209.8426)/10000 #263.0167 hectáreas
(length(which(valores>1456))*209.8426)/10000 #45.70372 hectáreas
# pendiente media
((1.432-.398)/62)*100
install.packages("rasterVis")
install.packages("rgl")
library(rasterVis)
library(rgl)
plot3D(DEM) # note: 3D not 3d
