forTal_2.R

# clear working directory
rm(list=ls())

library(raster)

## Loading required package: sp

library(ggplot2)

# First detrended DEM

setwd("~/Dropbox/BGU/Tal_Cohen/p_01_raster_values_quantiles/data") # setwd('N:\\Tal\\gis\\roughness\\ASCII\\19.10.16')
files = list.files(pattern = "\\.txt$")

rastervalues = NULL
final = NULL
for(inFile in files) {
  myIn = raster(inFile)
  myIn = myIn * 1000
  realData <- myIn[]
  myECDF <- ecdf(realData)
  x = seq(min(realData, na.rm=TRUE), max(realData, na.rm=TRUE), length.out=100)
  y = myECDF(x)
  dat = data.frame(
    x, y, 
    name = inFile,
    stringsAsFactors = FALSE
  )
  final = rbind(final, dat)
  val = data.frame(
      elev = realData,
      name = inFile,
      stringsAsFactors = FALSE
    )
    rastervalues = rbind(rastervalues, val)
}
rastervalues = rastervalues[!is.na(rastervalues$elev), ]

# Examine result
head(final)

##           x            y      name
## 1 -46.20079 1.447381e-06 bar_1.txt
## 2 -44.97872 2.272388e-04 bar_1.txt
## 3 -43.75665 9.132975e-04 bar_1.txt
## 4 -42.53458 1.635541e-03 bar_1.txt
## 5 -41.31251 2.586470e-03 bar_1.txt
## 6 -40.09044 3.903587e-03 bar_1.txt

tail(final)

##            x         y          name
## 795 48.99024 0.9996643 flat4_new.txt
## 796 49.99522 0.9997830 flat4_new.txt
## 797 51.00021 0.9998690 flat4_new.txt
## 798 52.00519 0.9999277 flat4_new.txt
## 799 53.01017 0.9999672 flat4_new.txt
## 800 54.01516 1.0000000 flat4_new.txt

head(rastervalues)

##         elev      name
## 859 33.23616 bar_1.txt
## 860 35.28116 bar_1.txt
## 861 36.87033 bar_1.txt
## 862 37.46021 bar_1.txt
## 863 38.74188 bar_1.txt
## 864 40.04066 bar_1.txt

tail(rastervalues)

##              elev          name
## 25959171 53.85053 flat4_new.txt
## 25959172 53.57552 flat4_new.txt
## 25959173 52.95861 flat4_new.txt
## 25959174 53.04146 flat4_new.txt
## 25959175 52.93798 flat4_new.txt
## 25959176 53.02095 flat4_new.txt

final$type = substr(final$name, 1, 3)
final$num = gsub('\\D+','', final$name)
rastervalues$type = substr(rastervalues$name, 1, 3)
rastervalues$num = gsub('\\D+','', rastervalues$name)

# Density of raster values
ggplot(rastervalues, aes(x = elev, fill = num)) +
  geom_density(alpha = 0.25) +
  facet_grid(type ~ .)

# Plot
ggplot(final, aes(x = x, y = y, colour = num, linetype = type)) +
  geom_line() +
  scale_x_continuous(
    "Size",
    limits = c(-150, 150), 
    breaks = seq(from = -150, to = 150, by = 50),
    labels = round(seq(from = -150, to = 150, by = 50), 2)
  ) +
  scale_y_continuous("Cumulative (%)") +
  theme_bw()

########################################################################

# clear working directory
rm(list=ls())

library(plyr)
library(raster)
library(ggplot2)
library(e1071)

## 
## Attaching package: 'e1071'

## The following object is masked from 'package:raster':
## 
##     interpolate

# First detrended DEM

setwd("~/Dropbox/BGU/Tal_Cohen/p_01_raster_values_quantiles/data") # setwd('N:\\Tal\\gis\\roughness\\ASCII\\19.10.16')
files = list.files(pattern = "\\.txt$")

types = substr(files, 1, 3)

rastervalues = NULL
final = NULL
for(type in types) {
  typevalues = NULL
  for(inFile in files[types == type]) {
    myIn = raster(inFile)
    myIn = myIn * 1000
    realData <- myIn[]
    typevalues = c(typevalues, realData)
  }
  myECDF <- ecdf(typevalues)
    x = seq(min(typevalues, na.rm=TRUE), max(typevalues, na.rm=TRUE), length.out=100)
    y = myECDF(x)
    dat = data.frame(
      x, y, 
      name = type,
      stringsAsFactors = FALSE
    )
    final = rbind(final, dat)
    val = data.frame(
      elev = typevalues,
      name = type,
      stringsAsFactors = FALSE
    )
    rastervalues = rbind(rastervalues, val)
}
rastervalues = rastervalues[!is.na(rastervalues$elev), ]

# Examine result
head(final)

##           x            y name
## 1 -85.51580 1.760382e-07  bar
## 2 -83.26362 8.801912e-07  bar
## 3 -81.01145 2.622970e-05  bar
## 4 -78.75928 1.431191e-04  bar
## 5 -76.50711 3.538369e-04  bar
## 6 -74.25493 6.140214e-04  bar

tail(final)

##            x         y name
## 795 80.63376 0.9999793  fla
## 796 81.97617 0.9999865  fla
## 797 83.31857 0.9999903  fla
## 798 84.66097 0.9999928  fla
## 799 86.00338 0.9999980  fla
## 800 87.34578 1.0000000  fla

head(rastervalues)

##         elev name
## 859 33.23616  bar
## 860 35.28116  bar
## 861 36.87033  bar
## 862 37.46021  bar
## 863 38.74188  bar
## 864 40.04066  bar

tail(rastervalues)

##               elev name
## 103836954 53.85053  fla
## 103836955 53.57552  fla
## 103836956 52.95861  fla
## 103836957 53.04146  fla
## 103836958 52.93798  fla
## 103836959 53.02095  fla

# Density of raster values
ggplot(rastervalues, aes(x = elev, fill = name)) +
  geom_density(alpha = 0.25)

# Plot
ggplot(final, aes(x = x, y = y, colour = name, linetype = name)) +
  geom_line() +
  scale_x_continuous(
    "Size",
    limits = c(-150, 150), 
    breaks = seq(from = -150, to = 150, by = 50),
    labels = round(seq(from = -150, to = 150, by = 50), 2)
  ) +
  scale_y_continuous("Cumulative (%)") +
  theme_bw()

# Statistics - ECDF
tab = ddply(
  final,
  "name",
  summarize,
  q05 = quantile(y, probs = 0.05),
  q16 = quantile(y, probs = 0.16),
  q50 = quantile(y, probs = 0.5),
  q84 = quantile(y, probs = 0.84),
  q95 = quantile(y, probs = 0.95),
  sd = sd(y),
  kurtosis = kurtosis(y)
)
tab

##   name          q05        q16      q50       q84       q95        sd
## 1  bar 0.0006010121 0.01265176 0.860720 0.9989522 0.9999556 0.4242226
## 2  fla 0.0008569272 0.02049738 0.970043 0.9999116 0.9999797 0.4267965
##    kurtosis
## 1 -1.609663
## 2 -1.463008

# write.csv(tab, "tab.csv", row.names = FALSE)

# Statistics - ECDF
tab1 = ddply(
  rastervalues,
  "name",
  summarize,
  q05 = quantile(elev, probs = 0.05),
  q16 = quantile(elev, probs = 0.16),
  q50 = quantile(elev, probs = 0.5),
  q84 = quantile(elev, probs = 0.84),
  q95 = quantile(elev, probs = 0.95),
  sd = sd(elev),
  kurtosis = kurtosis(elev),
  mean = mean(elev)
)
tab1

##   name       q05        q16       q50      q84      q95       sd kurtosis
## 1  bar -36.35633 -21.047231 0.7351041 23.74136 40.65323 23.75135 0.955106
## 2  fla -17.48431  -9.442985 0.5496144 10.99813 17.95375 10.98315 1.069536
##        mean
## 1 1.2543416
## 2 0.4978649

# write.csv(tab1, "tab1.csv", row.names = FALSE)

########################################################################

# clear working directory
rm(list=ls())

library(plyr)
library(raster)
library(ggplot2)
library(e1071)

# First detrended DEM

setwd("~/Dropbox/BGU/Tal_Cohen/p_01_raster_values_quantiles/data") # setwd('N:\\Tal\\gis\\roughness\\ASCII\\19.10.16')
files = list.files(pattern = "\\.txt$")

rastervalues = NULL
final = NULL
for(file in files) {
  filevalues = NULL
  myIn = raster(file)
  myIn = myIn * 1000
  realData <- myIn[]
  filevalues = c(filevalues, realData)
  myECDF <- ecdf(filevalues)
  x = seq(min(filevalues, na.rm=TRUE), max(filevalues, na.rm=TRUE), length.out=100)
  y = myECDF(x)
  dat = data.frame(
    x, y, 
    name = file,
    stringsAsFactors = FALSE
  )
  final = rbind(final, dat)
  val = data.frame(
    elev = filevalues,
    name = file,
    stringsAsFactors = FALSE
  )
  rastervalues = rbind(rastervalues, val)
}
rastervalues = rastervalues[!is.na(rastervalues$elev), ]

# Examine result
head(final)

##           x            y      name
## 1 -46.20079 1.447381e-06 bar_1.txt
## 2 -44.97872 2.272388e-04 bar_1.txt
## 3 -43.75665 9.132975e-04 bar_1.txt
## 4 -42.53458 1.635541e-03 bar_1.txt
## 5 -41.31251 2.586470e-03 bar_1.txt
## 6 -40.09044 3.903587e-03 bar_1.txt

tail(final)

##            x         y          name
## 795 48.99024 0.9996643 flat4_new.txt
## 796 49.99522 0.9997830 flat4_new.txt
## 797 51.00021 0.9998690 flat4_new.txt
## 798 52.00519 0.9999277 flat4_new.txt
## 799 53.01017 0.9999672 flat4_new.txt
## 800 54.01516 1.0000000 flat4_new.txt

head(rastervalues)

##         elev      name
## 859 33.23616 bar_1.txt
## 860 35.28116 bar_1.txt
## 861 36.87033 bar_1.txt
## 862 37.46021 bar_1.txt
## 863 38.74188 bar_1.txt
## 864 40.04066 bar_1.txt

tail(rastervalues)

##              elev          name
## 25959171 53.85053 flat4_new.txt
## 25959172 53.57552 flat4_new.txt
## 25959173 52.95861 flat4_new.txt
## 25959174 53.04146 flat4_new.txt
## 25959175 52.93798 flat4_new.txt
## 25959176 53.02095 flat4_new.txt

# Density of raster values
ggplot(rastervalues, aes(x = elev, fill = name)) +
  geom_density(alpha = 0.25)

# Plot
ggplot(final, aes(x = x, y = y, colour = name, linetype = name)) +
  geom_line() +
  scale_x_continuous(
    "Size",
    limits = c(-150, 150), 
    breaks = seq(from = -150, to = 150, by = 50),
    labels = round(seq(from = -150, to = 150, by = 50), 2)
  ) +
  scale_y_continuous("Cumulative (%)") +
  theme_bw()

# Statistics - ECDF
tab = ddply(
  final,
  "name",
  summarize,
  q05 = quantile(y, probs = 0.05),
  q16 = quantile(y, probs = 0.16),
  q50 = quantile(y, probs = 0.5),
  q84 = quantile(y, probs = 0.84),
  q95 = quantile(y, probs = 0.95),
  sd = sd(y),
  kurtosis = kurtosis(y)
)
tab

##            name          q05         q16       q50       q84       q95
## 1     bar_1.txt 0.0038377312 0.049618022 0.7884312 0.9938062 0.9989426
## 2      bar2.txt 0.0011887113 0.016334777 0.7725735 0.9981279 0.9996272
## 3      bar3.txt 0.0005071380 0.015174912 0.9203551 0.9984146 0.9993352
## 4  bar4newb.txt 0.0097851005 0.102079744 0.9184718 0.9990118 0.9999051
## 5     flat1.txt 0.0001037852 0.009522049 0.7089112 0.9984467 0.9996702
## 6 flat2_new.txt 0.0002096575 0.006352209 0.7634972 0.9995992 0.9998568
## 7 flat3_new.txt 0.0010685214 0.012811855 0.9740051 0.9997972 0.9999534
## 8 flat4_new.txt 0.0007050874 0.025828822 0.5383820 0.9962271 0.9996702
##          sd  kurtosis
## 1 0.3938868 -1.486669
## 2 0.4174689 -1.647658
## 3 0.4266219 -1.556260
## 4 0.3863469 -1.212720
## 5 0.4406260 -1.804811
## 6 0.4375298 -1.767853
## 7 0.4306220 -1.483118
## 8 0.4050482 -1.705278

# write.csv(tab, "tab.csv", row.names = FALSE)

# Statistics - ECDF
tab1 = ddply(
  rastervalues,
  "name",
  summarize,
  q05 = quantile(elev, probs = 0.05),
  q16 = quantile(elev, probs = 0.16),
  q50 = quantile(elev, probs = 0.5),
  q84 = quantile(elev, probs = 0.84),
  q95 = quantile(elev, probs = 0.95),
  sd = sd(elev),
  kurtosis = kurtosis(elev),
  mean = mean(elev)
)
tab1

##            name       q05        q16        q50       q84       q95
## 1     bar_1.txt -26.77466 -16.704657 -1.4427240 18.796151 32.765107
## 2      bar2.txt -39.10840 -25.364218  1.5677210 27.549557 44.003407
## 3      bar3.txt -25.65387 -14.742271  1.4896690 17.207030 32.009986
## 4  bar4newb.txt -42.07639 -23.525521 -0.8345842 24.438555 46.639464
## 5     flat1.txt -15.60789  -8.006394  0.2566576 10.677671 15.546059
## 6 flat2_new.txt -12.05655  -7.544558  0.1411438  6.055854  9.301591
## 7 flat3_new.txt -16.11497  -9.532332  0.4805326 10.646084 16.561929
## 8 flat4_new.txt -25.81322 -15.152033  2.4253130 17.694229 22.788409
##          sd   kurtosis       mean
## 1 18.245296  0.3814870  0.3352742
## 2 26.472995  0.2980146  1.4620484
## 3 18.013250  2.6557628  1.8575508
## 4 26.115906  0.8692506  0.5746597
## 5  9.899647  1.1822108  0.5520368
## 6  6.642285  0.1531999 -0.4801201
## 7 10.524702  1.4840665  0.4565730
## 8 15.315681 -0.4806180  1.3862402

# write.csv(tab1, "tab1.csv", row.names = FALSE)