Spectral Analysis

Spectral analysis

Leemos los datos…

library(tidyverse)

── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

✓ ggplot2 3.3.5     ✓ purrr   0.3.4
✓ tibble  3.1.6     ✓ dplyr   1.0.8
✓ tidyr   1.2.0     ✓ stringr 1.4.0
✓ readr   2.1.2     ✓ forcats 0.5.1

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

pixels <- read_csv("~/Downloads/new_data.csv")

Rows: 1094522 Columns: 14

── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): clase
dbl (13): B2, B3, B4, B5, B6, B7, B8, B11, B12, DEM, Slopes, lat, lon

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Registramos algunas funciones para obtener la desviacion estandar por debajo y por encima del promedio y para extraer los registros para cada clase…

##
getSDRange <- function(data, the_class, sd_factor){
  sample <- data |> 
    filter(clase == the_class) |> 
    mutate(BDEM = DEM, BSlope = Slopes) |>
    select(B2, B3, B4, B5, B6, B7, B8, B11, B12, BDEM, BSlope)
  print(summary(sample))
  print(sample |> tally())
  
  sds   = sapply(sample, sd)
  means = sapply(sample, mean)
  stats = bind_rows(means + (sd_factor * sds), means, means - (sd_factor * sds)) |> 
    mutate(id = row_number()) |> 
    pivot_longer(cols = starts_with("B"), 
                 names_to = "band", values_to = "value") |>
    mutate(band = factor(
      band, 
      levels = c("B2","B3","B4","B5","B6","B7","B8","B9","B11","B12", "BDEM", "BSlope"))
    )
  
  return(stats)
}

##
plotClass <- function(the_class){
  bands <- pixels |> 
    filter(clase == the_class) |> 
    mutate(BDEM = DEM, BSlope = Slopes) |>
    select(B2, B3, B4, B5, B6, B7, B8, B11, B12, BDEM, BSlope) |>
    mutate(id = row_number()) |>
    pivot_longer(cols = starts_with("B"), names_to = "band", values_to = "value") |>
    mutate(clase = the_class, 
           band = factor(
             band, 
             levels = c("B2","B3","B4","B5","B6","B7","B8","B9","B11","B12", "BDEM", "BSlope")
           )
    )
  
  return(bands)
}

Graficamos los cinco numeros para las clases principales…

palma <- plotClass("palma")
cafe1 <- plotClass("cafe_expuesto")
cafe2 <- plotClass("cafe_semisombra")
cafe3 <- plotClass("cafe_sombra")
ggplot(data = palma, aes(x=band, y= value, fill=clase)) + geom_boxplot()

ggplot(data = cafe1, aes(x=band, y= value, fill=clase)) + geom_boxplot()

ggplot(data = cafe2, aes(x=band, y= value, fill=clase)) + geom_boxplot()

ggplot(data = cafe3, aes(x=band, y= value, fill=clase)) + geom_boxplot()

y para las demas…

vege <- plotClass("vegetacion")
ciud <- plotClass("ciudad")
agua <- plotClass("agua")
ggplot(data = vege, aes(x=band, y= value, fill=clase)) + geom_boxplot()

ggplot(data = agua, aes(x=band, y= value, fill=clase)) + geom_boxplot()

ggplot(data = ciud, aes(x=band, y= value, fill=clase)) + geom_boxplot()

Vemos como se distribuye los 3 tipos de cafe…

##
cafe <- cafe1 |> bind_rows(cafe2) |> bind_rows(cafe3)
ggplot(data = cafe, aes(x=band, y= value, fill=clase)) + geom_boxplot()

Registramos una funcion para graficar las firmas espectrales para las clases con su promedio y +1.5 y -1.5 desviaciones estandar…

##
runGgplot <- function(data, the_class, factor){
  range = getSDRange(pixels, the_class, factor)
  g = ggplot(data = data, aes(x=band, y= value, group=id)) + 
    geom_line(size=0.2, alpha=0.01) +
    stat_summary(aes(y = value, group = 1), 
                 fun = mean, colour = "blue", geom = "line", group = 1) +
    geom_line(data = range |> filter(id == 1), 
              aes(x = band, y = value), colour= "red", linetype = "dashed") +
    geom_line(data = range |> filter(id == 3), 
              aes(x = band, y = value), colour= "red", linetype = "dashed") +
    theme_bw() 
  
  return(g)
}

Graficamos las firmas espectrales para las clases principales…

factor_default = 1.5

plot(runGgplot(palma, "palma", factor_default))

       B2               B3               B4               B5        
 Min.   :  59.0   Min.   : 167.0   Min.   :  66.0   Min.   :-100.0  
 1st Qu.: 170.0   1st Qu.: 376.0   1st Qu.: 218.0   1st Qu.: 668.0  
 Median : 203.0   Median : 421.0   Median : 269.0   Median : 742.0  
 Mean   : 229.5   Mean   : 451.8   Mean   : 377.2   Mean   : 797.9  
 3rd Qu.: 261.0   3rd Qu.: 496.0   3rd Qu.: 410.0   3rd Qu.: 875.0  
 Max.   :1594.0   Max.   :1910.0   Max.   :2276.0   Max.   :2434.0  
       B6             B7             B8            B11            B12        
 Min.   :-100   Min.   :  71   Min.   :  52   Min.   :  85   Min.   :  51.0  
 1st Qu.:2261   1st Qu.:2650   1st Qu.:2732   1st Qu.:1589   1st Qu.: 675.0  
 Median :2462   Median :3067   Median :3156   Median :1730   Median : 785.0  
 Mean   :2450   Mean   :2995   Mean   :3087   Mean   :1895   Mean   : 960.9  
 3rd Qu.:2650   3rd Qu.:3380   3rd Qu.:3486   3rd Qu.:2057   3rd Qu.:1064.0  
 Max.   :4146   Max.   :5358   Max.   :5492   Max.   :4945   Max.   :4208.0  
      BDEM           BSlope      
 Min.   : 35.0   Min.   : 0.000  
 1st Qu.: 78.0   1st Qu.: 2.502  
 Median : 93.0   Median : 3.701  
 Mean   : 93.5   Mean   : 4.177  
 3rd Qu.:106.0   3rd Qu.: 5.319  
 Max.   :239.0   Max.   :30.112  
# A tibble: 1 × 1
      n
  <int>
1 19019

plot(runGgplot(cafe1, "cafe_expuesto", factor_default))

       B2             B3               B4               B5        
 Min.   :   1   Min.   :  36.0   Min.   :   3.0   Min.   :-100.0  
 1st Qu.: 156   1st Qu.: 349.0   1st Qu.: 199.0   1st Qu.: 707.5  
 Median : 209   Median : 440.0   Median : 273.0   Median : 854.0  
 Mean   : 234   Mean   : 472.8   Mean   : 325.5   Mean   : 902.9  
 3rd Qu.: 277   3rd Qu.: 552.0   3rd Qu.: 393.0   3rd Qu.:1040.0  
 Max.   :3122   Max.   :3402.0   Max.   :2858.0   Max.   :4254.0  
       B6              B7              B8             B11      
 Min.   : -100   Min.   :  748   Min.   :  555   Min.   : 461  
 1st Qu.: 2078   1st Qu.: 2554   1st Qu.: 2654   1st Qu.:1610  
 Median : 2471   Median : 3020   Median : 3186   Median :1858  
 Mean   : 2620   Mean   : 3189   Mean   : 3384   Mean   :1928  
 3rd Qu.: 2998   3rd Qu.: 3620   3rd Qu.: 3868   3rd Qu.:2141  
 Max.   :10928   Max.   :12494   Max.   :12792   Max.   :8339  
      B12              BDEM          BSlope     
 Min.   : 195.0   Min.   :-100   Min.   : 1.68  
 1st Qu.: 733.0   1st Qu.:1365   1st Qu.:21.00  
 Median : 892.0   Median :1592   Median :27.39  
 Mean   : 961.8   Mean   :1560   Mean   :27.54  
 3rd Qu.:1105.0   3rd Qu.:1782   3rd Qu.:33.97  
 Max.   :5257.0   Max.   :2714   Max.   :65.51  
# A tibble: 1 × 1
      n
  <int>
1 16467

plot(runGgplot(cafe2, "cafe_semisombra", factor_default))

       B2               B3               B4               B5        
 Min.   :   1.0   Min.   :   1.0   Min.   :   1.0   Min.   :-100.0  
 1st Qu.: 128.0   1st Qu.: 316.0   1st Qu.: 184.0   1st Qu.: 657.0  
 Median : 174.0   Median : 397.0   Median : 250.0   Median : 791.0  
 Mean   : 192.5   Mean   : 426.3   Mean   : 299.5   Mean   : 830.3  
 3rd Qu.: 232.0   3rd Qu.: 504.0   3rd Qu.: 357.0   3rd Qu.: 956.0  
 Max.   :3299.0   Max.   :3300.0   Max.   :5316.0   Max.   :3794.0  
       B6              B7              B8             B11       
 Min.   : -100   Min.   :  738   Min.   :  515   Min.   :  601  
 1st Qu.: 2124   1st Qu.: 2623   1st Qu.: 2712   1st Qu.: 1574  
 Median : 2459   Median : 3084   Median : 3242   Median : 1815  
 Mean   : 2618   Mean   : 3260   Mean   : 3453   Mean   : 1891  
 3rd Qu.: 2943   3rd Qu.: 3703   3rd Qu.: 3958   3rd Qu.: 2116  
 Max.   :10474   Max.   :13264   Max.   :15440   Max.   :10586  
      B12              BDEM          BSlope      
 Min.   : 239.0   Min.   :-100   Min.   : 1.989  
 1st Qu.: 695.0   1st Qu.:1200   1st Qu.:20.968  
 Median : 837.0   Median :1427   Median :27.535  
 Mean   : 907.4   Mean   :1410   Mean   :27.585  
 3rd Qu.:1036.0   3rd Qu.:1621   3rd Qu.:34.383  
 Max.   :7736.0   Max.   :2277   Max.   :64.963  
# A tibble: 1 × 1
      n
  <int>
1 17690

plot(runGgplot(cafe3, "cafe_sombra", factor_default))

       B2               B3               B4               B5        
 Min.   :   1.0   Min.   :   8.0   Min.   :   2.0   Min.   : 116.0  
 1st Qu.: 116.0   1st Qu.: 296.0   1st Qu.: 172.0   1st Qu.: 625.0  
 Median : 158.0   Median : 376.0   Median : 228.0   Median : 745.0  
 Mean   : 172.5   Mean   : 400.3   Mean   : 264.7   Mean   : 784.5  
 3rd Qu.: 210.0   3rd Qu.: 474.0   3rd Qu.: 311.0   3rd Qu.: 896.0  
 Max.   :1644.0   Max.   :2124.0   Max.   :1910.0   Max.   :3940.0  
       B6              B7              B8             B11      
 Min.   :  485   Min.   : 1106   Min.   :  611   Min.   : 626  
 1st Qu.: 2184   1st Qu.: 2717   1st Qu.: 2768   1st Qu.:1568  
 Median : 2563   Median : 3198   Median : 3356   Median :1783  
 Mean   : 2707   Mean   : 3360   Mean   : 3545   Mean   :1858  
 3rd Qu.: 3062   3rd Qu.: 3826   3rd Qu.: 4090   3rd Qu.:2065  
 Max.   :11659   Max.   :13811   Max.   :15440   Max.   :7743  
      B12              BDEM          BSlope      
 Min.   : 232.0   Min.   : 666   Min.   : 1.585  
 1st Qu.: 681.0   1st Qu.:1124   1st Qu.:21.405  
 Median : 803.0   Median :1264   Median :27.375  
 Mean   : 856.8   Mean   :1326   Mean   :27.417  
 3rd Qu.: 964.0   3rd Qu.:1488   3rd Qu.:33.501  
 Max.   :3853.0   Max.   :2390   Max.   :61.055  
# A tibble: 1 × 1
      n
  <int>
1 18641

y para las secundarias…

plot(runGgplot(vege, "vegetacion", factor_default))

       B2               B3               B4               B5        
 Min.   :   1.0   Min.   :   1.0   Min.   :   1.0   Min.   :  43.0  
 1st Qu.: 138.0   1st Qu.: 345.0   1st Qu.: 227.0   1st Qu.: 670.0  
 Median : 200.0   Median : 405.0   Median : 316.0   Median : 765.0  
 Mean   : 198.7   Mean   : 416.9   Mean   : 420.5   Mean   : 811.1  
 3rd Qu.: 253.0   3rd Qu.: 477.0   3rd Qu.: 519.0   3rd Qu.: 915.0  
 Max.   :1415.0   Max.   :2480.0   Max.   :3458.0   Max.   :3595.0  
       B6              B7              B8             B11            B12      
 Min.   :  172   Min.   :  284   Min.   :  321   Min.   : 168   Min.   :  72  
 1st Qu.: 1701   1st Qu.: 2061   1st Qu.: 2142   1st Qu.:1689   1st Qu.: 780  
 Median : 1963   Median : 2420   Median : 2534   Median :1917   Median : 958  
 Mean   : 2164   Mean   : 2552   Mean   : 2695   Mean   :2013   Mean   :1079  
 3rd Qu.: 2401   3rd Qu.: 2821   3rd Qu.: 2976   3rd Qu.:2282   3rd Qu.:1284  
 Max.   :10589   Max.   :11119   Max.   :13864   Max.   :5520   Max.   :4352  
      BDEM            BSlope      
 Min.   :  40.0   Min.   : 0.000  
 1st Qu.:  75.0   1st Qu.: 3.356  
 Median : 177.0   Median : 8.788  
 Mean   : 566.6   Mean   :16.201  
 3rd Qu.: 856.0   3rd Qu.:28.604  
 Max.   :2472.0   Max.   :63.070  
# A tibble: 1 × 1
      n
  <int>
1 59523

plot(runGgplot(ciud, "ciudad", factor_default))

       B2              B3              B4              B5             B6      
 Min.   :    1   Min.   :  257   Min.   :   83   Min.   :-100   Min.   :-100  
 1st Qu.:  736   1st Qu.: 1068   1st Qu.: 1088   1st Qu.:1649   1st Qu.:2242  
 Median : 1036   Median : 1382   Median : 1476   Median :1912   Median :2443  
 Mean   : 1087   Mean   : 1432   Mean   : 1514   Mean   :1936   Mean   :2455  
 3rd Qu.: 1368   3rd Qu.: 1732   3rd Qu.: 1888   3rd Qu.:2184   3rd Qu.:2656  
 Max.   :12872   Max.   :12288   Max.   :12152   Max.   :6748   Max.   :6691  
       B7             B8             B11            B12            BDEM       
 Min.   : 494   Min.   :  405   Min.   : 815   Min.   : 413   Min.   :-100.0  
 1st Qu.:2442   1st Qu.: 2418   1st Qu.:2483   1st Qu.:2080   1st Qu.: 148.0  
 Median :2664   Median : 2700   Median :2748   Median :2416   Median : 158.0  
 Mean   :2694   Mean   : 2734   Mean   :2761   Mean   :2399   Mean   : 158.1  
 3rd Qu.:2920   3rd Qu.: 3017   3rd Qu.:3022   3rd Qu.:2742   3rd Qu.: 167.0  
 Max.   :6235   Max.   :11376   Max.   :7704   Max.   :8313   Max.   : 210.0  
     BSlope      
 Min.   : 0.000  
 1st Qu.: 2.104  
 Median : 3.066  
 Mean   : 3.515  
 3rd Qu.: 4.482  
 Max.   :28.471  
# A tibble: 1 × 1
       n
   <int>
1 139487

plot(runGgplot(agua, "agua", factor_default))

       B2            B3             B4               B5        
 Min.   : 56   Min.   : 103   Min.   :  88.0   Min.   :  94.0  
 1st Qu.:160   1st Qu.: 222   1st Qu.: 190.0   1st Qu.: 213.0  
 Median :181   Median : 253   Median : 346.0   Median : 253.0  
 Mean   :212   Mean   : 340   Mean   : 342.9   Mean   : 415.4  
 3rd Qu.:232   3rd Qu.: 405   3rd Qu.: 435.0   3rd Qu.: 456.0  
 Max.   :744   Max.   :1060   Max.   :1854.0   Max.   :1708.0  
       B6               B7             B8            B11            B12        
 Min.   :  33.0   Min.   :   1   Min.   :  31   Min.   :  60   Min.   :  38.0  
 1st Qu.: 120.0   1st Qu.: 175   1st Qu.: 147   1st Qu.: 126   1st Qu.:  88.0  
 Median : 149.0   Median :1925   Median :1956   Median :1473   Median : 736.0  
 Mean   : 558.8   Mean   :1661   Mean   :1680   Mean   :1187   Mean   : 649.3  
 3rd Qu.: 358.0   3rd Qu.:2421   3rd Qu.:2532   3rd Qu.:1979   3rd Qu.:1039.0  
 Max.   :3959.0   Max.   :5601   Max.   :5758   Max.   :3347   Max.   :2576.0  
      BDEM           BSlope      
 Min.   : 33.0   Min.   : 0.000  
 1st Qu.: 52.0   1st Qu.: 1.673  
 Median : 57.0   Median : 3.234  
 Mean   : 59.5   Mean   : 3.602  
 3rd Qu.: 57.0   3rd Qu.: 4.978  
 Max.   :138.0   Max.   :24.419  
# A tibble: 1 × 1
      n
  <int>
1 39341

Comparamos los promedios entre clases principales…

ggplot(data = cafe, aes(x=band, y= value, group = id)) + 
  stat_summary(data = palma, aes(y = value, group = 1), 
               fun = mean, colour = "orange", geom = "line", group = 1) +
  stat_summary(data = cafe1, aes(y = value, group = 1), 
               fun = mean, colour = "green", geom = "line", group = 1) +
  stat_summary(data = cafe2, aes(y = value, group = 1), 
               fun = mean, colour = "blue", geom = "line", group = 1) +
  stat_summary(data = cafe3, aes(y = value, group = 1), 
               fun = mean, colour = "red", geom = "line", group = 1) +
  theme_bw()

Comparamos cafe con las clases secundarias…

ggplot(data = cafe, aes(x=band, y= value, group = id)) + 
  stat_summary(data = cafe1, aes(y = value, group = 1), 
               fun = mean, colour = "green", geom = "line", group = 1) +
  stat_summary(data = cafe2, aes(y = value, group = 1), 
               fun = mean, colour = "blue", geom = "line", group = 1) +
  stat_summary(data = cafe3, aes(y = value, group = 1), 
               fun = mean, colour = "red", geom = "line", group = 1) +
  stat_summary(data = vege, aes(y = value, group = 1), 
               fun = mean, colour = "darkgreen", geom = "line", group = 1) +
  stat_summary(data = ciud, aes(y = value, group = 1), 
               fun = mean, colour = "darkgray", geom = "line", group = 1) +
  stat_summary(data = agua, aes(y = value, group = 1), 
               fun = mean, colour = "darkblue", geom = "line", group = 1) +
  theme_bw()

Comparamos palma con las clases secundarias…

ggplot(data = palma, aes(x=band, y= value, group = id)) + 
  stat_summary(data = palma, aes(y = value, group = 1), 
               fun = mean, colour = "orange", geom = "line", group = 1) +
  stat_summary(data = vege, aes(y = value, group = 1), 
               fun = mean, colour = "darkgreen", geom = "line", group = 1) +
  stat_summary(data = ciud, aes(y = value, group = 1), 
               fun = mean, colour = "darkgray", geom = "line", group = 1) +
  stat_summary(data = agua, aes(y = value, group = 1), 
               fun = mean, colour = "darkblue", geom = "line", group = 1) +
  theme_bw()