KD - 1_Databehandling

1. Databehandling av Klimadata (Kjøres via Rstudio serveren)

1.1 Nedbør - Mengde

Total mengde nedbør (mm) i perioden Mai - Oktober.

# Set up for parallel run 
UseCores = 12
cl = makeCluster(UseCores)
registerDoParallel(cl)

# List files in R:/ directory (daily data; rr = precipitation)
setwd("/home/NINA.NO/markus.israelsen/Mounts/")
AllFilesInRR = list.files("/data/R/GeoSpatialData/Meteorology/Norway_SeNorge/Original/rr",full.names = TRUE, pattern = ".bil$")

# Define years  (1957 - 2020)
uniq_y = c(1957:2020) 

# Define months
uniq_m = c("05", "06", "07", "08", "09", "10") # May - October

# Subsetting the dataset by selecting only the months we are interested in (all years)
SampleMonths = foreach(i = 1:length(AllFilesInRR), .combine = c) %dopar% {
  if((substring(AllFilesInRR[i], first = nchar(AllFilesInRR[i])-8, last = nchar(AllFilesInRR[i])-7)) %in% uniq_m){ # if the month in line "i" is one of the            months we are interested in, append that line to SampleMonths
    AllFilesInRR[i] 
  }
}

# Summing the amount of precipitation for each year
sumPrecip_all_years = list()
sumPrecip_all_years = foreach(i = 1:length(uniq_y), .packages = c("doParallel", "raster")) %dopar% {
  sumPrecip_all_years = stack(SampleMonths[grepl(paste("rr", uniq_y[i], sep = "_"), SampleMonths, fixed = TRUE)])
  sumPrecip_year = list() # New code
  sumPrecip_year = foreach(j = 1:nlayers(sumPrecip_all_years), .packages = c("doParallel", "raster")) %dopar% {
    rs = raster(sumPrecip_all_years, layer = j)
    rs[rs > 0] = rs[rs > 0] / 10
    sumPrecip_year[[j]] = rs
  }
  sumPrecip_all_years[[i]] = stackApply(stack(sumPrecip_year), c(1), fun = sum)
}

# Stacking the years to one raster stack
sumPrecip_all_years = stack(sumPrecip_all_years) 

names(sumPrecip_all_years) = uniq_y # Name by years

# Rasterstack directory;
RasterLocation <- paste("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/", paste(paste("rr_sum_", uniq_y[1], uniq_y[length(uniq_y)], sep = "_"),".grd", sep = ""), sep = "")

# Write raster
writeRaster(sumPrecip_all_years, RasterLocation, format="raster", overwrite = TRUE)
    
stopImplicitCluster()

1.2 Nedbør - Dager

Antall dager med nedbør i perioden Mai - Oktober.

# Set up for parallel run 
UseCores = 12
cl = makeCluster(UseCores)
registerDoParallel(cl)

# List files in R:/ directory (daily data; rr = precipitation)
setwd("/home/NINA.NO/markus.israelsen/Mounts/")
AllFilesInRR = list.files("/data/R/GeoSpatialData/Meteorology/Norway_SeNorge/Original/rr",full.names = TRUE, pattern = ".bil$")

# Define years  (1957 - 2020)
uniq_y = c(1957:2020) 

# Define months
uniq_m = c("05", "06", "07", "08", "09", "10") # May - October

# Subsetting the dataset by selecting only the months we are interested in (all years)
SampleMonths = foreach(i = 1:length(AllFilesInRR), .combine = c) %dopar% {
  if((substring(AllFilesInRR[i], first = nchar(AllFilesInRR[i])-8, last = nchar(AllFilesInRR[i])-7)) %in% uniq_m){ # if the month in line "i" is one of the            months we are interested in, append that line to SampleMonths
  AllFilesInRR[i] 
  }
}

# Summing the days with precipitation for each year
daysPrecip_all_years = list()
daysPrecip_all_years = foreach(i = 1:length(uniq_y), .packages = c("doParallel", "raster")) %dopar% {
  
  daysPrecip_year = stack(SampleMonths[grepl(paste("rr", uniq_y[i], sep = "_"), SampleMonths, fixed = TRUE)])
  daysPrecip = list()
  daysPrecip = foreach(j = 1:nlayers(daysPrecip_year), .packages = c("doParallel", "raster")) %dopar% {
    r = raster(daysPrecip_year, layer = j)
    r[r > 0] = 1 
    daysPrecip[[j]] = r
  }
  
  daysPrecip_all_years[[i]] = stackApply(stack(daysPrecip), c(1), fun = sum)
}

# Stacking the years to one raster stack
daysPrecip_all_years = stack(daysPrecip_all_years) 

names(daysPrecip_all_years) = uniq_y # Name by years

# Rasterstack directory; change
RasterLocation <- paste("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/", paste(paste("rr_days_", uniq_y[1], uniq_y[length(uniq_y)], sep = "_"),".grd", sep = ""), sep = "")

# Write raster
writeRaster(daysPrecip_all_years, RasterLocation, format="raster", overwrite = TRUE)

stopImplicitCluster()

1.3 Temperatur - Sommer

Gjennomsnittlig sommertemperatur, \(^\circ\)C (Juni - Aug.).

# Set up for parallel run 
UseCores = 12 
cl = makeCluster(UseCores)
registerDoParallel(cl)

# List files in R:/ directory (daily data; tm = temperature)
setwd("/home/NINA.NO/markus.israelsen/Mounts/")
AllFilesInTM <- list.files("/data/R/GeoSpatialData/Meteorology/Norway_SeNorge/Original/tm",full.names = TRUE, pattern = ".bil$")

# Define years
uniq_y <- c(1957:2020) 

# Define months
uniq_m <- c("06", "07", "08") # July - August

SampleMonths = foreach(i = 1:length(AllFilesInTM), .combine = c, .packages = c("doParallel", "raster")) %dopar% {
  if((substring(AllFilesInTM[i], first = nchar(AllFilesInTM[i])-8, last = nchar(AllFilesInTM[i])-7)) %in% uniq_m){ # if the month in line "i" is one of the months we are interested in, append that line to SampleMonths
  AllFilesInTM[i] 
  }
}

meanSummer_all_years = list()

meanSummer_all_years = foreach(i = 1:length(uniq_y), .packages = c("doParallel", "raster")) %dopar% {
  meanSummer_years = stack(SampleMonths[grepl(paste("tm", uniq_y[i], sep = "_"), SampleMonths, fixed = TRUE)])
  
  meanSummer_year = list()
  meanSummer_year = foreach(j = 1:nlayers(meanSummer_years), .packages = c("doParallel", "raster")) %dopar% {
    ms = raster(meanSummer_years, layer = j)
    ms[ms > 0] = (ms[ms > 0] - 2731) / 10 # convert from 1/10 kelvin to degrees centigrade
    meanSummer_year[[j]] = ms
  }
  meanSummer_all_years[[i]] = stackApply(stack(meanSummer_year), c(1), fun = mean) 
}

# Stacking the years to one raster stack
meanSummer_all_years = stack(meanSummer_all_years) 

names(meanSummer_all_years) = uniq_y # Name by years

# Rasterstack directory; change
RasterLocation <- paste("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/",
                    paste(paste("tm_summer_", uniq_y[1], uniq_y[length(uniq_y)], sep = "_"),".grd", sep = ""), sep = "")

# Write raster
writeRaster(meanSummer_all_years, RasterLocation, format="raster", overwrite = TRUE)

stopImplicitCluster()

1.4 Temperatur - Vinter

Gjennomsnittlig vintertemperatur, \(^\circ\)C (Des. - Feb.).

# Set up for parallel run 
UseCores = 12 
cl = makeCluster(UseCores)
registerDoParallel(cl)

# List files in R:/ directory (daily data; tm = temperature)
setwd("/home/NINA.NO/markus.israelsen/Mounts/")
AllFilesInTM <- list.files("/data/R/GeoSpatialData/Meteorology/Norway_SeNorge/Original/tm",full.names = TRUE, pattern = ".bil$")

# Define years
uniq_y <- c(1957:2020) 

# Calculating the mean winter temperature for each year
uniq_mw = c("12", "01", "02") # Define months for winter

# Subsetting the dataset by selecting only the months we are interested in (all years)
SampleMonths = foreach(i = 1:length(AllFilesInTM), .combine = c, .packages = c("doParallel", "raster")) %dopar% {
  if((substring(AllFilesInTM[i], first = nchar(AllFilesInTM[i])-8, last = nchar(AllFilesInTM[i])-7)) %in% uniq_mw){ # if the month in line "i" is one of the months we are interested in, append that line to SampleMonths
  AllFilesInTM[i] 
  }
}

meanWinter_all_years = list()

meanWinter_all_years = foreach(i = 2:length(uniq_y), .packages = c("doParallel", "raster", "rlist")) %dopar% {
  
  meanWinter_years = SampleMonths[grep(paste("tm", uniq_y[i]-1, 12, sep = "_"), SampleMonths)]
  meanWinter_years = list.append(meanWinter_years, SampleMonths[grep(paste("tm", uniq_y[i], uniq_mw[2], sep = "_"), SampleMonths)], SampleMonths[grep(paste("tm", uniq_y[i], uniq_mw[3], sep = "_"), SampleMonths)])
  meanWinter_years = stack(meanWinter_years)
  
  meanWinter_year = list()
  meanWinter_year = foreach(j = 1:nlayers(meanWinter_years), .packages = c("doParallel", "raster")) %dopar% {
    mw = raster(meanWinter_years, layer = j)
    mw[mw > 0] = (mw[mw > 0] - 2731) / 10 # convert from 1/10 kelvin to degrees centigrade
    meanWinter_year[[j]] = mw
  }
  meanWinter_all_years[[i-1]] = stackApply(stack(meanWinter_year), c(1), fun = mean)
}

  
# Stacking the years to one raster stack
meanWinter_all_years = stack(meanWinter_all_years) 

layerNames = c()
for(i in 1:length(uniq_y)){
  layerNames = append(layerNames, 
                      paste(substring(uniq_y[i], first = nchar(uniq_y[i])-1, last = nchar(uniq_y[i])), substring(uniq_y[i+1], first = nchar(uniq_y[i+1])-1, last = nchar(uniq_y[i+1])), sep = "-"))
}
layerNames = layerNames[1:(length(uniq_y)-1)]
names(meanWinter_all_years) = layerNames # Name by years

# Rasterstack directory; change
RasterLocation <- paste("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/",
                    paste(paste("tm_winter_", uniq_y[1], uniq_y[length(uniq_y)], sep = "_"),".grd", sep = ""), sep = "")

# Write raster
writeRaster(meanWinter_all_years, RasterLocation, format="raster", overwrite = TRUE)

stopImplicitCluster()

1.5 Snødekke

Antall dager med snødekke i perioden Oktober - Juni.

# Set up for parallel run 
UseCores = 12 
cl = makeCluster(UseCores)
registerDoParallel(cl)

# List files in R:/ directory (daily data; rr = precipitation)
setwd("/home/NINA.NO/markus.israelsen/Mounts/")
AllFilesInSD = list.files("/data/R/GeoSpatialData/Meteorology/Norway_SeNorge/Original/sd",full.names = TRUE, pattern = ".bil$")

# Define years
uniq_y = c(1957:2020) 

# Define months
uniq_m = c("10", "11", "12", "01", "02", "03", "04", "05", "06") # Oct - June (Data starts in Sept 1957, should therefore start at i = 2 with mean of Oct 1957+Jan 58+Feb 58+Mar 58+Apr 58+May 58+June 58)

# Subsetting the dataset by selecting only the months we are interested in (all years)
SampleMonths = foreach(i = 1:length(AllFilesInSD), .combine = c, .packages = c("doParallel", "raster")) %dopar% {
  if((substring(AllFilesInSD[i], first = nchar(AllFilesInSD[i])-8, last = nchar(AllFilesInSD[i])-7)) %in% uniq_m){ # if the month in line "i" is one of the months we are interested in, append that line to SampleMonths
    AllFilesInSD[i] 
  }
}


# Calculating the number of days with snow cover for each year
daysSnowCover_all_years = list()

daysSnowCover_all_years = foreach(i = 2:length(uniq_y), .packages = c("doParallel","raster", "rlist")) %dopar% {
# if (i == 1){
#     daysSnowCover_year = SampleMonths[grepl(paste("sd", uniq_y[i], sep = "_"), SampleMonths, fixed = TRUE)]
#     daysSnowCover_year = list.append(daysSnowCover_year, SampleMonths[grepl(paste("sd", uniq_y[i+1], uniq_m[4], sep = "_"), SampleMonths, fixed = TRUE)], SampleMonths[grepl(paste("sd", uniq_y[i+1], uniq_m[5], sep = "_"), SampleMonths, fixed = TRUE)], SampleMonths[grepl(paste("sd", uniq_y[i+1], uniq_m[6], sep = "_"), SampleMonths, fixed = TRUE)], SampleMonths[grepl(paste("sd", uniq_y[i+1], uniq_m[7], sep = "_"), SampleMonths, fixed = TRUE)], SampleMonths[grepl(paste("sd", uniq_y[i+1], uniq_m[8], sep = "_"), SampleMonths, fixed = TRUE)], SampleMonths[grepl(paste("sd", uniq_y[i+1], uniq_m[9], sep = "_"), SampleMonths, fixed = TRUE)])
#     daysSnowCover_year = stack(daysSnowCover_year)
#     daysSnowCover = list()
#     daysSnowCover = foreach(j = 1:nlayers(daysSnowCover_year), .packages = c("doParallel","raster")) %dopar% {
#       r = raster(daysSnowCover_year, layer = j)
#       r[r > 0] = 1
#       daysSnowCover[[j]] = r
#     }
#   daysSnowCover_all_years[[i]] = stackApply(stack(daysSnowCover), c(1), fun = sum)
# }


  daysSnowCover_year = SampleMonths[grepl(paste("sd", uniq_y[i]-1, 12, sep = "_"), SampleMonths, fixed = TRUE) 
                                    | grepl(paste("sd", uniq_y[i]-1, 11, sep = "_"), SampleMonths, fixed = TRUE) 
                                    | grepl(paste("sd", uniq_y[i]-1, 10, sep = "_"), SampleMonths, fixed = TRUE)]
  daysSnowCover_year = list.append(daysSnowCover_year, SampleMonths[grepl(paste("sd", uniq_y[i], uniq_m[4], sep = "_"), SampleMonths, fixed = TRUE)], SampleMonths[grepl(paste("sd", uniq_y[i], uniq_m[5], sep = "_"), SampleMonths, fixed = TRUE)], SampleMonths[grepl(paste("sd", uniq_y[i], uniq_m[6], sep = "_"), SampleMonths, fixed = TRUE)], SampleMonths[grepl(paste("sd", uniq_y[i], uniq_m[7], sep = "_"), SampleMonths, fixed = TRUE)], SampleMonths[grepl(paste("sd", uniq_y[i], uniq_m[8], sep = "_"), SampleMonths, fixed = TRUE)], SampleMonths[grepl(paste("sd", uniq_y[i], uniq_m[9], sep = "_"), SampleMonths, fixed = TRUE)])
  daysSnowCover_year = stack(daysSnowCover_year)
  daysSnowCover = list()
  daysSnowCover = foreach(j = 1:nlayers(daysSnowCover_year), .packages = c("doParallel","raster")) %dopar% {
    r = raster(daysSnowCover_year, layer = j)
    r[r > 0] = 1
    daysSnowCover[[j]] = r
  }
  daysSnowCover_all_years[[i-1]] = stackApply(stack(daysSnowCover), c(1), fun = sum)
}


# Stacking the years to one raster stack
daysSnowCover_all_years = stack(daysSnowCover_all_years) 

layerNames = c()
for(i in 1:length(uniq_y)){
  layerNames = append(layerNames, paste(substring(uniq_y[i], first = nchar(uniq_y[i])-1, last = nchar(uniq_y[i])), substring(uniq_y[i+1], first = nchar(uniq_y[i+1])-1, last = nchar(uniq_y[i+1])), sep="-"))
}
layerNames = layerNames[1:(length(uniq_y)-1)]
names(daysSnowCover_all_years) = layerNames # Name by years

# Rasterstack directory; 
RasterLocation <- paste("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", paste(paste("sd_s", uniq_y[1], uniq_y[length(uniq_y)], sep = "_"),".grd", sep = ""), sep = "")

# Write raster
writeRaster(daysSnowCover_all_years, RasterLocation, format="raster", overwrite = TRUE)

stopImplicitCluster()

# if(SimpleOrComplex == "Complex"){
#   
#   # Temperature
#   AllFilesInTM <- list.files("/data/R/GeoSpatialData/Meteorology/Norway_SeNorge/Original/tm",full.names = TRUE, pattern = ".bil$")
# 
#   SampleMonthsTemp = foreach(i = 1:length(AllFilesInTM), .combine = c, .packages = c("doParallel", "raster")) %dopar% {
#     if((substring(AllFilesInTM[i], first = nchar(AllFilesInTM[i])-8, last = nchar(AllFilesInTM[i])-7)) %in% uniq_m){
#       AllFilesInTM[i]
#     }
#   }
# 
#   SnowTempCover_all_years = list()
# 
#   SnowTempCover_all_years = foreach(i = 1:length(uniq_y), .packages = c("doParallel", "raster", "rlist")) %dopar% {
#     if (i == 1){
#       SnowRaster = SampleMonths[grepl(paste("sd", uniq_y[i], sep = "_"), SampleMonths, fixed = TRUE)]
#       TempRaster = SampleMonthsTemp[grepl(paste("tm", uniq_y[i], sep = "_"), SampleMonthsTemp, fixed = TRUE)]
#       SnowRaster = stack(SnowRaster)
#       TempRaster = stack(TempRaster)
# 
#       SnowTempCover = list()
# 
#       SnowTempCover = foreach(j = 1:nlayers(SnowRaster), .packages = c("doParallel", "raster")) %dopar% {
#         r = raster(SnowRaster, layer = j)
# 
#         k = raster(TempRaster, layer = j)
# 
#         TempRows = which(k[] > 2731)
#         SnowRows = which(r[] > 0)
# 
#         if(length(TempRows) >= length(SnowRows)){
#           TempSnowRows = which(SnowRows %in% TempRows)
#           r[][SnowRows[TempSnowRows]] = -1337
#           r[r != -1337] = 0
#           r[r == -1337] = 1
#         }
#         else {
#           TempSnowRows = which(TempRows %in% SnowRows)
#           r[][TempRows[TempSnowRows]] = -1337
#           r[r != -1337] = 0
#           r[r == -1337] = 1
#         }
# 
#         SnowTempCover[[j]] = r
#       }
#       SnowTempCover_all_years[[i]] = stackApply(stack(SnowTempCover), c(1), fun = sum)
#     }
# 
#     else if(i > 1) {
#       SnowRaster = SampleMonths[grepl(paste("sd", uniq_y[i]-1, 12, sep = "_"), SampleMonths, fixed = TRUE)]
#       SnowRaster = list.append(SnowRaster, SampleMonths[grepl(paste("sd", uniq_y[i], uniq_m[2], sep = "_"), SampleMonths, fixed = TRUE)], .data = SampleMonths[grepl(paste("sd", uniq_y[i], uniq_m[4], sep = "_"), SampleMonths, fixed = TRUE)], SampleMonths[grepl(paste("sd", uniq_y[i], uniq_m[3], sep = "_"), SampleMonths, fixed = TRUE)])
#       TempRaster = SampleMonthsTemp[grepl(paste("tm", uniq_y[i]-1, 12, sep = "_"), SampleMonthsTemp, fixed = TRUE)]
#       TempRaster = list.append(TempRaster, SampleMonthsTemp[grepl(paste("tm", uniq_y[i], uniq_m[2], sep = "_"), SampleMonthsTemp, fixed = TRUE)], SampleMonthsTemp[grepl(paste("tm", uniq_y[i], uniq_m[3], sep = "_"), SampleMonthsTemp, fixed = TRUE)], SampleMonthsTemp[grepl(paste("tm", uniq_y[i], uniq_m[4], sep = "_"), SampleMonthsTemp, fixed = TRUE)])
#       SnowRaster = stack(SnowRaster)
#       TempRaster = stack(TempRaster)
#   
#       SnowTempCover = list()
#   
#       SnowTempCover = foreach(j = 1:nlayers(SnowRaster), .packages = c("doParallel", "raster")) %dopar% {
#         r = raster(SnowRaster, layer = j)
#     
#         k = raster(TempRaster, layer = j)
#     
#         TempRows = which(k[] > 2731)
#         SnowRows = which(r[] > 0)
#     
#         if(length(TempRows) >= length(SnowRows)){
#           TempSnowRows = which(SnowRows %in% TempRows)
#           r[][SnowRows[TempSnowRows]] = -1337 # random variable value
#           r[r != -1337] = 0
#           r[r == -1337] = 1
#         }
#         else {
#           TempSnowRows = which(TempRows %in% SnowRows)
#           r[][TempRows[TempSnowRows]] = -1337
#           r[r != -1337] = 0
#           r[r == -1337] = 1
#         }
#     
#         SnowTempCover[[j]] = r
#       }
#       SnowTempCover_all_years[[i]] = stackApply(stack(SnowTempCover), c(1), fun = sum)
#     }
#   }
# 
#   SnowTempCover_all_years = stack(SnowTempCover_all_years)
# 
#   names(SnowTempCover_all_years) = uniq_y # Name by years
# 
#   # Rasterstack directory
#   RasterLocation = paste("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/M/Klima/All years 1957 - 2020/Days Snow Cover Condition on Temperature/", paste(paste("sd_c", uniq_y[1], uniq_y[length(uniq_y)], sep = "_"), ".grd", sep = ""), sep = "")
# 
#   # Write Raster
#   writeRaster(SnowTempCover_all_years, RasterLocation, format = "raster", overwrite = TRUE)
#   }

1.6 Vekstsesong

Vekstsesongens lengde i antall dager (uten snødekke og døgntemp. > 5\(^\circ\)C).

# Set up for parallel run
UseCores = 12
cl = makeCluster(UseCores)
registerDoParallel(cl)

# Snow cover
AllFilesInSD = list.files("/data/R/GeoSpatialData/Meteorology/Norway_SeNorge/Original/sd", full.names = TRUE, pattern = ".bil$")

# Define years 
uniq_y = c(1957:2020)

SampleDataSnow = AllFilesInSD[grepl(paste("sd_", uniq_y[1], sep=""), AllFilesInSD, fixed = TRUE)]
for(i in 2:length(uniq_y)){
  SampleDataSnow = list.append(SampleDataSnow, AllFilesInSD[grepl(paste("sd_", uniq_y[i], sep=""), AllFilesInSD, fixed = TRUE)])
}

# Temperature
AllFilesInTM <- list.files("/data/R/GeoSpatialData/Meteorology/Norway_SeNorge/Original/tm", full.names = TRUE, pattern = ".bil$")

SampleDataTemp = AllFilesInTM[grepl(paste("tm_", uniq_y[1], sep = ""), AllFilesInTM, fixed = TRUE)]
for(i in 2:length(uniq_y)){
  SampleDataTemp = list.append(SampleDataTemp, AllFilesInTM[grepl(paste("tm_", uniq_y[i], sep = ""), AllFilesInTM, fixed = TRUE)])
}

rm(AllFilesInTM)

GrowthSeason_all_years = list()

GrowthSeason_all_years = foreach(i = 1:length(uniq_y), .packages = c("doParallel", "raster", "rlist")) %dopar% {
  
  SnowRaster = list()
  TempRaster = list()
  
  SnowRaster[[i]] = stack(SampleDataSnow[grep(paste("sd_", uniq_y[i], sep=""), SampleDataSnow)])
  TempRaster[[i]] = stack(SampleDataTemp[grep(paste("tm_", uniq_y[i], sep=""), SampleDataTemp)])
    
  SnowTempCover = list()
  
  SnowTempCover = foreach(j = 1:nlayers(SnowRaster[[i]]), .packages = c("doParallel", "raster")) %dopar% {
    r = raster(SnowRaster[[i]], layer = j)
    
    k = raster(TempRaster[[i]], layer = j)
    
    r_values = values(r)
    
    k_values = values(k)
    k_values_rows = which(k_values >= 2781)
    
    test_val = r_values[k_values_rows]
    test_val = which(test_val != 0)
    if(length(test_val) > 0){
      k_values_rows = k_values_rows[-test_val]  
    }
    
    r_values[k_values_rows] = -1337
    r_values[r_values != -1337] = 0
    r_values[r_values == -1337] = 1
    values(r) = r_values
    
    SnowTempCover[[j]] = r
    
    # TempRows = which(k[] >= 2781)
    # SnowRows = which(r[] == 0)
    # 
    # if(length(TempRows) >= length(SnowRows)){
    #   TempSnowRows = which(SnowRows %in% TempRows)
    #   r[][SnowRows[TempSnowRows]] = -1337
    #   r[r != -1337] = 0
    #   r[r == -1337] = 1
    # }
    # 
    # else {
    #   TempSnowRows = which(TempRows %in% SnowRows)
    #   r[][TempRows[TempSnowRows]] = -1337
    #   r[r != -1337] = 0
    #   r[r == -1337] = 1
    # }
    
    SnowTempCover[[j]] = r
  }
  GrowthSeason_all_years[[i]] = stackApply(stack(SnowTempCover), c(1), fun = sum)
}

GrowthSeason_all_years = stack(GrowthSeason_all_years)

names(GrowthSeason_all_years) = uniq_y # Name by years

# Rasterstack directory
RasterLocation = paste("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", paste(paste("sd_gw", uniq_y[1], uniq_y[length(uniq_y)], sep = "_"), ".grd", sep = ""), sep = "")

# Write Raster
writeRaster(GrowthSeason_all_years, RasterLocation, format = "raster", overwrite = TRUE)

stopImplicitCluster()

1.7 Vinterregn

Total mengde nedbør (mm) for dager med temperatur > 2\(^\circ\)C i perioden Januar - Mars.

  rm(list = ls())
  # Set up for parallel run, more efficient 
  UseCores = 12 
  cl = makeCluster(UseCores)
  registerDoParallel(cl)
  
  # List files in R:/ directory (daily data; rr = precipitation)
  WD = getwd() 
  setwd("/home/NINA.NO/markus.israelsen/Mounts/")
  AllFilesInRR = list.files("/data/R/GeoSpatialData/Meteorology/Norway_SeNorge/Original/rr",full.names = TRUE, pattern = ".bil$")
  AllFilesInRR = AllFilesInRR[-grep("arome", AllFilesInRR)] # remove the two obs that are written as "arome06_rr_2015_11_20.bil" (duplicated with the two that are written normally "rr_2015_11_20.bil")
  
  # Define years
  uniq_y = c(1957:2020) 

  # Define months
  uniq_m = c("01", "02", "03")
  
  # Subsetting the precipitation dataset by selecting only the months we are interested in (all years)
  SampleMonths = foreach(i = 1:length(AllFilesInRR), .combine = c, .packages = c("doParallel", "raster")) %dopar% {
    if((substring(AllFilesInRR[i], first = nchar(AllFilesInRR[i])-8, last = nchar(AllFilesInRR[i])-7)) %in% uniq_m){ # if the month in line "i" is one of the months we are interested in, append that line to SampleMonths
      AllFilesInRR[i] 
    }
  }
  
  # Get the temperature data
  AllFilesInTM <- list.files("/data/R/GeoSpatialData/Meteorology/Norway_SeNorge/Original/tm",full.names = TRUE, pattern = ".bil$")
  
  # Subsetting the temperature dataset by selecting only the months we are interested in (all years)
  SampleMonthsTemp = foreach(i = 1:length(AllFilesInTM), .combine = c, .packages = c("doParallel", "raster")) %dopar% {
      if((substring(AllFilesInTM[i], first = nchar(AllFilesInTM[i])-8, last = nchar(AllFilesInTM[i])-7)) %in% uniq_m){
        AllFilesInTM[i]
      }
    }
  
  rm(AllFilesInTM) # Remove the full temperature dataset from memory as it is quite big and is not needed anymore
  
  # 
  RainTempCover_all_years = list()

  RainTempCover_all_years = foreach(i = 1:length(uniq_y), .packages = c("doParallel", "raster", "rlist")) %dopar% {
    
    RainRaster = list()
    TempRaster = list()
    
    RainRaster[[i]] = stack(SampleMonths[grep(paste("rr_",uniq_y[i], sep=""), SampleMonths)])
    TempRaster[[i]] = stack(SampleMonthsTemp[grep(paste("tm_",uniq_y[i],sep=""), SampleMonthsTemp)])
    
    RainTempCover = list()
    RainTempCover = foreach(j = 1:nlayers(RainRaster[[i]]), .packages = c("doParallel", "raster")) %dopar% {
      r = raster(RainRaster[[i]], layer = j)
      r_values = values(r)
      
      k = raster(TempRaster[[i]], layer = j)
      k_values = values(k)
      k_values_rows = which(k_values <= 2751.5)
      
      r_values[k_values_rows] = NA
      values(r) = r_values
      
      RainTempCover[[j]] = r 
    }
    RainTempCover_all_years[[i]] = stackApply(stack(RainTempCover), c(1), fun = sum) # Sum the amount of precipitation per year (Jan. - March)
  }

  WinterRain_all_years = stack(RainTempCover_all_years)
  
  names(WinterRain_all_years) = uniq_y # Name by years

  # Rasterstack directory
  RasterLocation = paste("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", paste(paste("sd_vr", uniq_y[1], uniq_y[length(uniq_y)], sep = "_"), ".grd", sep = ""), sep = "")

  # Write Raster
  writeRaster(WinterRain_all_years, RasterLocation, format = "raster", overwrite = TRUE)
  
  stopImplicitCluster()

1.8 Snødybde

Gjennomsnittlig snødybde (mm) i perioden Desember - Mai.

# Set up for parallel run 
rm(list = ls())
UseCores = 12 
cl = makeCluster(UseCores)
registerDoParallel(cl)

# List files in R:/ directory
setwd("/home/NINA.NO/markus.israelsen/Mounts/")
AllFilesInSD = list.files("/data/R/GeoSpatialData/Meteorology/Norway_SeNorge/Original/sd", full.names = TRUE, pattern = ".bil$")

# Define years
uniq_y = c(1957:2020) 

# Define months
uniq_m = c("12", "01", "02", "03", "04", "05") # Des - May

# Subsetting the dataset by selecting only the months we are interested in (all years)
SampleMonths = foreach(i = 1:length(AllFilesInSD), .combine = c, .packages = c("doParallel", "raster")) %dopar% {
  if((substring(AllFilesInSD[i], first = nchar(AllFilesInSD[i])-8, last = nchar(AllFilesInSD[i])-7)) %in% uniq_m){ # if the month in line "i" is one of the months we are interested in, append that line to SampleMonths
    AllFilesInSD[i] 
  }
}

# Calculating the number of days with snow cover for each year
snowDepth_all_years = list()

snowDepth_all_years = foreach(i = 2:length(uniq_y), .packages = c("doParallel","raster", "rlist")) %dopar% {

  snowDepth_year = SampleMonths[grepl(paste("sd", uniq_y[i]-1, 12, sep = "_"), SampleMonths, fixed = TRUE)]
  
  snowDepth_year = list.append(snowDepth_year, SampleMonths[grepl(paste("sd", uniq_y[i], uniq_m[2], sep = "_"), SampleMonths, fixed = TRUE)], SampleMonths[grepl(paste("sd", uniq_y[i], uniq_m[3], sep = "_"), SampleMonths, fixed = TRUE)], SampleMonths[grepl(paste("sd", uniq_y[i], uniq_m[4], sep = "_"), SampleMonths, fixed = TRUE)], SampleMonths[grepl(paste("sd", uniq_y[i], uniq_m[5], sep = "_"), SampleMonths, fixed = TRUE)], SampleMonths[grepl(paste("sd", uniq_y[i], uniq_m[6], sep = "_"), SampleMonths, fixed = TRUE)])
  
  snowDepth_year = stack(snowDepth_year)
  snowDepth_all_years[[i-1]] = stackApply(stack(snowDepth_year), c(1), fun = mean)
  
}

# Stacking the years to one raster stack
snowDepth_all_years = stack(snowDepth_all_years) 

layerNames = c()
for(i in 1:length(uniq_y)){
  layerNames = append(layerNames, paste(substring(uniq_y[i], first = nchar(uniq_y[i])-1, last = nchar(uniq_y[i])), substring(uniq_y[i+1], first = nchar(uniq_y[i+1])-1, last = nchar(uniq_y[i+1])), sep="-"))
}
layerNames = layerNames[1:(length(uniq_y)-1)]
names(snowDepth_all_years) = layerNames # Name by years

# Rasterstack directory; 
RasterLocation <- paste("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/", paste(paste("sd_snowdepth", uniq_y[1], uniq_y[length(uniq_y)], sep = "_"),".grd", sep = ""), sep = "")

# Write raster
writeRaster(snowDepth_all_years, RasterLocation, format="raster", overwrite = TRUE)

stopImplicitCluster()

2. Google Earth Engine (GEE)

Beregn antall “fjellpixler” og legg layer til infrastruktur raster (Må kjøres i GEE).

var infra = ee.Image("users/tsimonjakobsson/ecocond_2020-2021/NY_INFRA_IND"),
    norway_wgs84 = ee.FeatureCollection("users/tsimonjakobsson/ecocond_2020-2021/Norway_wgs84"),
    regnorway_wgs84 = ee.FeatureCollection("users/tsimonjakobsson/ecocond_2020-2021/regNorway_wgs84"),
    ecosystem_map = ee.Image("users/zandersamuel/NINA/Raster/Norway_ecosystem_types_Simon_5m"),
    infra_wgs84 = ee.Image("users/markusfisraelsen/Infra_wgs84"),
    infra_epsg = ee.Image("users/markusfisraelsen/Infra_ESPG25833"),
    ecoViz = {"min":101,"max":952,"palette":["#00911d","#bcbcbc","#f2e341","#eb56ff","#c2efff","#75b3ff","#2163ff","#3252a8","#ff0000"]};
    

// From Zander
function reduceImgResolution(image, reducer, projection){
  return image.reduceResolution({
    reducer: reducer, 
    bestEffort: false,  // for best accuracy, this should be left as false
    maxPixels: 420
  }).reproject(projection)
}


// To show the different infrastructure index values as a colour gradient
var infraViz = {min: 0, max: 15.38416862487793, palette: ['00bbbb', '0000bb']};

// Map infrastructure
var infra_viz = infra_epsg.visualize(infraViz);
//Map.addLayer(infra_viz,{}, 'infra_all');

// Map ecosystem
var ecosystem_viz = ecosystem_map.visualize(ecoViz);
//Map.addLayer(ecosystem_viz,{}, 'ecosystem_viz');

// Get the projection for infrastructure and ecosystem
var infra_projection = infra_viz.projection();
var ecosystem_projection = ecosystem_map.projection();


// Filter to get regions
var region1 = regnorway_wgs84.filter(ee.Filter.eq("Region_1", 1));
var region2 = regnorway_wgs84.filter(ee.Filter.eq("Region_1", 2));
var region3 = regnorway_wgs84.filter(ee.Filter.eq("Region_1", 3));
var region4 = regnorway_wgs84.filter(ee.Filter.eq("Region_1", 4));
var region5 = regnorway_wgs84.filter(ee.Filter.eq("Region_1", 5));

var region1_geo = region1.geometry();
var region2_geo = region2.geometry();
var region3_geo = region3.geometry();
var region4_geo = region4.geometry();
var region5_geo = region5.geometry();

// Clip infra for each region
var infra_region1 = infra_epsg.clip(region1_geo);
var infra_region2 = infra_epsg.clip(region2_geo);
var infra_region3 = infra_epsg.clip(region3_geo);
var infra_region4 = infra_epsg.clip(region4_geo);
var infra_region5 = infra_epsg.clip(region5_geo);
//Map.addLayer(infra_region1.visualize(infraViz));

// Print the native scale of the infra layer
var infra_scale = infra_region1.projection().nominalScale();
//print('Infra scale in meters', infra_scale);

// Clip ecosystem for each region
var eco_region1 = ecosystem_map.clip(region1_geo);
var eco_region2 = ecosystem_map.clip(region2_geo);
var eco_region3 = ecosystem_map.clip(region3_geo);
var eco_region4 = ecosystem_map.clip(region4_geo);
var eco_region5 = ecosystem_map.clip(region5_geo);
//Map.addLayer(eco_region1.visualize(ecoViz));

// Print the native scale of the ecosystem layer
var eco_scale = eco_region1.projection().nominalScale();
//print('Eco scale in meters', eco_scale);

// Mask out the mountain ecosystem
var mount_eco_region1 = eco_region1.eq(201).or(eco_region1.eq(202)).selfMask();
var mount_mask_region1 = ecosystem_map.updateMask(mount_eco_region1);

var mount_eco_region2 = eco_region2.eq(201).or(eco_region2.eq(202)).selfMask();
var mount_mask_region2 = ecosystem_map.updateMask(mount_eco_region2);

var mount_eco_region3 = eco_region3.eq(201).or(eco_region3.eq(202)).selfMask();
var mount_mask_region3 = ecosystem_map.updateMask(mount_eco_region3);

var mount_eco_region4 = eco_region4.eq(201).or(eco_region4.eq(202)).selfMask();
var mount_mask_region4 = ecosystem_map.updateMask(mount_eco_region4);

var mount_eco_region5 = eco_region5.eq(201).or(eco_region5.eq(202)).selfMask();
var mount_mask_region5 = ecosystem_map.updateMask(mount_eco_region5);
//Map.addLayer(mount_mask_region2.visualize({bands: "b1", palette: ["005500", "00ff00"]}));

// Get the count of mountain pixels for each region
var mountCount_region1 = reduceImgResolution(mount_mask_region1, ee.Reducer.count(), infra_projection);
mountCount_region1 = mountCount_region1.unmask(0);

var mountCount_region2 = reduceImgResolution(mount_mask_region2, ee.Reducer.count(), infra_projection);
mountCount_region2 = mountCount_region2.unmask(0);

var mountCount_region3 = reduceImgResolution(mount_mask_region3, ee.Reducer.count(), infra_projection);
mountCount_region3 = mountCount_region3.unmask(0);

var mountCount_region4 = reduceImgResolution(mount_mask_region4, ee.Reducer.count(), infra_projection);
mountCount_region4 = mountCount_region4.unmask(0);

var mountCount_region5 = reduceImgResolution(mount_mask_region5, ee.Reducer.count(), infra_projection);
mountCount_region5 = mountCount_region5.unmask(0);

// Append the mountain pixel count to the infra index raster
var mount_float_region1 = mountCount_region1.float(); // Convert the mountain count layer to float to be identical to the infra raster 
var mount_infra_region1 = infra_region1.rename('natureIndex')
  .addBands(mount_float_region1.rename('mountCount'));

var mount_float_region2 = mountCount_region2.float();
var mount_infra_region2 = infra_region2.rename('natureIndex')
  .addBands(mount_float_region2.rename('mountCount'));
  
var mount_float_region3 = mountCount_region3.float();
var mount_infra_region3 = infra_region3.rename('natureIndex')
  .addBands(mount_float_region3.rename('mountCount'));

var mount_float_region4 = mountCount_region4.float();
var mount_infra_region4 = infra_region4.rename('natureIndex')
  .addBands(mount_float_region4.rename('mountCount'));

var mount_float_region5 = mountCount_region5.float();
var mount_infra_region5 = infra_region5.rename('natureIndex')
  .addBands(mount_float_region5.rename('mountCount'));

//print(infra_projection);

// Exporting the image for each region
Export.image.toDrive({
  image: mount_infra_region1,
  description: "mount_infra_region1",
  region: mount_infra_region1.geometry(),
  scale: 100,
  crs: infra_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: mount_infra_region2,
  description: "mount_infra_region2",
  region: mount_infra_region2.geometry(),
  scale: 100,
  crs: infra_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: mount_infra_region3,
  description: "mount_infra_region3",
  region: mount_infra_region3.geometry(),
  scale: 100,
  crs: infra_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: mount_infra_region4,
  description: "mount_infra_region4",
  region: mount_infra_region4.geometry(),
  scale: 100,
  crs: infra_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: mount_infra_region5,
  description: "mount_infra_region5",
  region: mount_infra_region5.geometry(),
  scale: 100,
  crs: infra_projection,
  maxPixels: 1e10
});

3. Konvertere klimavariablene til en filtype støttet av Google Earth Engine

Konvertering av raster filene fra “.grd” filtype til “.tif”.

# Sum nedbør
rDir_snCon = "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/" 
snCon = stack(paste(rDir_snCon, "rr_sum_1957_2020.grd", sep=""))
crs(snCon) = CRS("+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs")
writeRaster(snCon, paste(rDir_snCon, "rr_sum_utm_1957_2020.grd", sep = ""), format = "GTiff")

# Dager med nedbør
rDir_dnCon = "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/" 
dnCon = stack(paste(rDir_dnCon, "rr_days_1957_2020.grd", sep=""))
crs(dnCon) = CRS("+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs")
writeRaster(dnCon, paste(rDir_dnCon, "rr_days_utm_1957_2020.grd", sep = ""), format = "GTiff")

# Temperatur - sommer
rDir_tsCon = "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/" 
tsCon = stack(paste(rDir_tsCon, "tm_summer_1957_2020.grd", sep=""))
crs(tsCon) = CRS("+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs")
writeRaster(tsCon, paste(rDir_tsCon, "tm_summer_utm_1957_2020.grd", sep = ""), format = "GTiff")

# Temperatur - vinter
rDir_tvCon = "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/" 
tvCon = stack(paste(rDir_tvCon, "tm_winter_1957_2020.grd", sep=""))
crs(tvCon) = CRS("+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs")
writeRaster(tvCon, paste(rDir_tvCon, "tm_winter_utm_1957_2020.grd", sep = ""), format = "GTiff")

# Snødekke
rDir_scCon = "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/" 
scCon = stack(paste(rDir_scCon, "sd_s_1957_2020.grd", sep=""))
crs(scCon) = CRS("+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs")
writeRaster(scCon, paste(rDir_scCon, "sd_s_utm_1957_2020.grd", sep = ""), format = "GTiff")

# Vekstsesong
rDir_vCon = "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/" 
vCon = stack(paste(rDir_vCon, "sd_gw_1957_2020.grd", sep=""))
crs(vCon) = CRS("+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs")
writeRaster(vCon, paste(rDir_vCon, "sd_gw_utm_1957_2020.grd", sep = ""), format = "GTiff")

# Vinterregn
rDir_vrCon = "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/" 
vrCon = stack(paste(rDir_vrCon, "sd_vr_1957_2020.grd", sep=""))
crs(vrCon) = CRS("+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs")
writeRaster(vrCon, paste(rDir_vrCon, "sd_vr_utm_1957_2020.grd", sep = ""), format = "GTiff", overwrite = TRUE)

# Snødybde
rDir_sdCon = "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/"
sdCon = stack(paste(rDir_sdCon, "sd_snowdepth_1957_2020.grd", sep = ""))
crs(sdCon) = CRS("+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs")
writeRaster(sdCon, paste(rDir_sdCon, "snowDepth_utm_1957_2020.grd", sep = ""), format = "GTiff", overwrite = TRUE)

4. Begrense klimavariablene til de ulike regionene i Norge (Må kjøres i GEE).

var infra_epsg = ee.Image("users/markusfisraelsen/Infra_ESPG25833"),
    ecosystem_map = ee.Image("users/zandersamuel/NINA/Raster/Norway_ecosystem_types_Simon_5m"),
    regnorway_wgs84 = ee.FeatureCollection("users/tsimonjakobsson/ecocond_2020-2021/regNorway_wgs84"),
    Mount_daysPrecip = ee.Image("users/markusfisraelsen/Mount_daysPrecip"),
    Mount_daysSnowCover = ee.Image("users/markusfisraelsen/Mount_daysSnowCover"),
    Mount_growthSeason = ee.Image("users/markusfisraelsen/Mount_growthSeason"),
    Mount_meanSummer = ee.Image("users/markusfisraelsen/Mount_meanSummer"),
    Mount_meanWinter = ee.Image("users/markusfisraelsen/Mount_meanWinter"),
    Mount_winterRain = ee.Image("users/markusfisraelsen/Mount_winterRain"),
    Mount_sumPrecip = ee.Image("users/markusfisraelsen/Mount_sumPrecip");



// To show the different infrastructure index values as a colour gradient
var infraViz = {min: 0, max: 15.38416862487793, palette: ['00bbbb', '0000bb']};

// Map infrastructure
var infra_viz = infra_epsg.visualize(infraViz);
//Map.addLayer(infra_viz,{}, 'infra_all');

// Get the projection for infrastructure, ecosystem and days precipitation
var infra_projection = infra_viz.projection();
var ecosystem_projection = ecosystem_map.projection();
var daysPrecip_projection = Mount_daysPrecip.projection(); // This projection is used for all climate variables

// Filter to get regions
var region1 = regnorway_wgs84.filter(ee.Filter.eq("Region_1", 1));
var region2 = regnorway_wgs84.filter(ee.Filter.eq("Region_1", 2));
var region3 = regnorway_wgs84.filter(ee.Filter.eq("Region_1", 3));
var region4 = regnorway_wgs84.filter(ee.Filter.eq("Region_1", 4));
var region5 = regnorway_wgs84.filter(ee.Filter.eq("Region_1", 5));

var region1_geo = region1.geometry();
var region2_geo = region2.geometry();
var region3_geo = region3.geometry();
var region4_geo = region4.geometry();
var region5_geo = region5.geometry();


// Clip sum precipitation for each region
var sumPrecip_region1 = Mount_sumPrecip.clip(region1_geo);
var sumPrecip_region2 = Mount_sumPrecip.clip(region2_geo);
var sumPrecip_region3 = Mount_sumPrecip.clip(region3_geo);
var sumPrecip_region4 = Mount_sumPrecip.clip(region4_geo);
var sumPrecip_region5 = Mount_sumPrecip.clip(region5_geo);

// Clip daysPrecip for each region
var daysPrecip_region1 = Mount_daysPrecip.clip(region1_geo);
var daysPrecip_region2 = Mount_daysPrecip.clip(region2_geo);
var daysPrecip_region3 = Mount_daysPrecip.clip(region3_geo);
var daysPrecip_region4 = Mount_daysPrecip.clip(region4_geo);
var daysPrecip_region5 = Mount_daysPrecip.clip(region5_geo);

// Clip mean summer temperature for each region
var meanSummer_region1 = Mount_meanSummer.clip(region1_geo);
var meanSummer_region2 = Mount_meanSummer.clip(region2_geo);
var meanSummer_region3 = Mount_meanSummer.clip(region3_geo);
var meanSummer_region4 = Mount_meanSummer.clip(region4_geo);
var meanSummer_region5 = Mount_meanSummer.clip(region5_geo);

// Clip mean winter temperature for each region
var meanWinter_region1 = Mount_meanWinter.clip(region1_geo);
var meanWinter_region2 = Mount_meanWinter.clip(region2_geo);
var meanWinter_region3 = Mount_meanWinter.clip(region3_geo);
var meanWinter_region4 = Mount_meanWinter.clip(region4_geo);
var meanWinter_region5 = Mount_meanWinter.clip(region5_geo);

// Clip daysSnowCover for each region
var daysSnowCover_region1 = Mount_daysSnowCover.clip(region1_geo);
var daysSnowCover_region2 = Mount_daysSnowCover.clip(region2_geo);
var daysSnowCover_region3 = Mount_daysSnowCover.clip(region3_geo);
var daysSnowCover_region4 = Mount_daysSnowCover.clip(region4_geo);
var daysSnowCover_region5 = Mount_daysSnowCover.clip(region5_geo);

// Clip growthSeason for each region
var growthSeason_region1 = Mount_growthSeason.clip(region1_geo);
var growthSeason_region2 = Mount_growthSeason.clip(region2_geo);
var growthSeason_region3 = Mount_growthSeason.clip(region3_geo);
var growthSeason_region4 = Mount_growthSeason.clip(region4_geo);
var growthSeason_region5 = Mount_growthSeason.clip(region5_geo);

// Clip winterRain for each region
var winterRain_region1 = Mount_winterRain.clip(region1_geo);
var winterRain_region2 = Mount_winterRain.clip(region2_geo);
var winterRain_region3 = Mount_winterRain.clip(region3_geo);
var winterRain_region4 = Mount_winterRain.clip(region4_geo);
var winterRain_region5 = Mount_winterRain.clip(region5_geo);

print(Mount_sumPrecip.projection());
print(daysPrecip_projection);


// Export the sumPrecip climate variables for each region
Export.image.toDrive({
  image: sumPrecip_region1,
  description: "Mount_sumPrecip_region1",
  scale: 1000,
  region: sumPrecip_region1.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: sumPrecip_region2,
  description: "Mount_sumPrecip_region2",
  scale: 1000,
  region: sumPrecip_region2.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: sumPrecip_region3,
  description: "Mount_sumPrecip_region3",
  scale: 1000,
  region: sumPrecip_region3.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: sumPrecip_region4,
  description: "Mount_sumPrecip_region4",
  scale: 1000,
  region: sumPrecip_region4.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: sumPrecip_region5,
  description: "Mount_sumPrecip_region5",
  scale: 1000,
  region: sumPrecip_region5.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

// Export the daysPrecip climate variables for each region
Export.image.toDrive({
  image: daysPrecip_region1,
  description: "Mount_daysPrecip_region1",
  scale: 1000,
  region: daysPrecip_region1.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: daysPrecip_region2,
  description: "Mount_daysPrecip_region2",
  scale: 1000,
  region: daysPrecip_region2.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: daysPrecip_region3,
  description: "Mount_daysPrecip_region3",
  scale: 1000,
  region: daysPrecip_region3.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: daysPrecip_region4,
  description: "Mount_daysPrecip_region4",
  scale: 1000,
  region: daysPrecip_region4.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: daysPrecip_region5,
  description: "Mount_daysPrecip_region5",
  scale: 1000,
  region: daysPrecip_region5.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

// Export the meanSummer climate variables for each region
Export.image.toDrive({
  image: meanSummer_region1,
  description: "Mount_meanSummer_region1",
  scale: 1000,
  region: meanSummer_region1.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: meanSummer_region2,
  description: "Mount_meanSummer_region2",
  scale: 1000,
  region: meanSummer_region2.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: meanSummer_region3,
  description: "Mount_meanSummer_region3",
  scale: 1000,
  region: meanSummer_region3.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: meanSummer_region4,
  description: "Mount_meanSummer_region4",
  scale: 1000,
  region: meanSummer_region4.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: meanSummer_region5,
  description: "Mount_meanSummer_region5",
  scale: 1000,
  region: meanSummer_region5.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

// Export the meanWinter climate variables for each region
Export.image.toDrive({
  image: meanWinter_region1,
  description: "Mount_meanWinter_region1",
  scale: 1000,
  region: meanWinter_region1.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: meanWinter_region2,
  description: "Mount_meanWinter_region2",
  scale: 1000,
  region: meanWinter_region2.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: meanWinter_region3,
  description: "Mount_meanWinter_region3",
  scale: 1000,
  region: meanWinter_region3.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: meanWinter_region4,
  description: "Mount_meanWinter_region4",
  scale: 1000,
  region: meanWinter_region4.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: meanWinter_region5,
  description: "Mount_meanWinter_region5",
  scale: 1000,
  region: meanWinter_region5.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

// Export the snowCover climate variables for each region
Export.image.toDrive({
  image: daysSnowCover_region1,
  description: "Mount_daysSnowCover_region1",
  scale: 1000,
  region: daysSnowCover_region1.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: daysSnowCover_region2,
  description: "Mount_daysSnowCover_region2",
  scale: 1000,
  region: daysSnowCover_region2.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: daysSnowCover_region3,
  description: "Mount_daysSnowCover_region3",
  scale: 1000,
  region: daysSnowCover_region3.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: daysSnowCover_region4,
  description: "Mount_daysSnowCover_region4",
  scale: 1000,
  region: daysSnowCover_region4.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: daysSnowCover_region5,
  description: "Mount_daysSnowCover_region5",
  scale: 1000,
  region: daysSnowCover_region5.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

// Export the growthSeason climate variables for each region
Export.image.toDrive({
  image: growthSeason_region1,
  description: "Mount_growthSeason_region1",
  scale: 1000,
  region: growthSeason_region1.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: growthSeason_region2,
  description: "Mount_growthSeason_region2",
  scale: 1000,
  region: growthSeason_region2.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: growthSeason_region3,
  description: "Mount_growthSeason_region3",
  scale: 1000,
  region: growthSeason_region3.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: growthSeason_region4,
  description: "Mount_growthSeason_region4",
  scale: 1000,
  region: growthSeason_region4.geometry(),
  //crs: Mount_growthSeason_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: growthSeason_region5,
  description: "Mount_growthSeason_region5",
  scale: 1000,
  region: growthSeason_region5.geometry(),
  //crs: Mount_growthSeason_projection,
  maxPixels: 1e10
});

// Export the winterRain climate variables for each region
Export.image.toDrive({
  image: winterRain_region1,
  description: "Mount_winterRain_region1",
  scale: 1000,
  region: winterRain_region1.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: winterRain_region2,
  description: "Mount_winterRain_region2",
  scale: 1000,
  region: winterRain_region2.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: winterRain_region3,
  description: "Mount_winterRain_region3",
  scale: 1000,
  region: winterRain_region3.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: winterRain_region4,
  description: "Mount_winterRain_region4",
  scale: 1000,
  region: winterRain_region4.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

Export.image.toDrive({
  image: winterRain_region5,
  description: "Mount_winterRain_region5",
  scale: 1000,
  region: winterRain_region5.geometry(),
  //crs: daysPrecip_projection,
  maxPixels: 1e10
});

5. Slå sammen behandlede klimadata og rasterne av fjellområdet (Kjøres via Rstudio lokalt)

Loading the mountain count raster for each region

# Mountain Count
mountCount_region1_full = raster("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/mountInfra/mount_infra_region1.tif", band = 2)
mountCount_region2_full = raster("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/mountInfra/mount_infra_region2.tif", band = 2)
mountCount_region3_full = raster("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/mountInfra/mount_infra_region3.tif", band = 2)
mountCount_region4_full = raster("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/mountInfra/mount_infra_region4.tif", band = 2)
mountCount_region5_full = raster("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/mountInfra/mount_infra_region5.tif", band = 2)

# Sum precipitation region 1 1957 - 2020
sumPrecip_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/", full.names = TRUE, pattern = "region1.tif$")
sumPrecip_region1 = stack(sumPrecip_region1) 
sumExtent_region1 = extent(sumPrecip_region1) 
sumShp_region1 = as(sumExtent_region1, "SpatialPolygons") 
crs(sumShp_region1) = crs(sumPrecip_region1)

# Sum precipitation region 2 1957 - 2020
sumPrecip_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/", full.names = TRUE, pattern = "region2.tif$")
sumPrecip_region2 = stack(sumPrecip_region2) 
sumExtent_region2 = extent(sumPrecip_region2) 
sumShp_region2 = as(sumExtent_region2, "SpatialPolygons") 
crs(sumShp_region2) = crs(sumPrecip_region2)

# Sum precipitation region 3 1957 - 2020
sumPrecip_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/", full.names = TRUE, pattern = "region3.tif$")
sumPrecip_region3 = stack(sumPrecip_region3) 
sumExtent_region3 = extent(sumPrecip_region3) 
sumShp_region3 = as(sumExtent_region3, "SpatialPolygons") 
crs(sumShp_region3) = crs(sumPrecip_region3)

# Sum precipitation region 4 1957 - 2020
sumPrecip_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/", full.names = TRUE, pattern = "region4.tif$")
sumPrecip_region4 = stack(sumPrecip_region4) 
sumExtent_region4 = extent(sumPrecip_region4) 
sumShp_region4 = as(sumExtent_region4, "SpatialPolygons") 
crs(sumShp_region4) = crs(sumPrecip_region4)

# Sum precipitation region 5 1957 - 2020
sumPrecip_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/", full.names = TRUE, pattern = "region5.tif$")
sumPrecip_region5 = stack(sumPrecip_region5) 
sumExtent_region5 = extent(sumPrecip_region5) 
sumShp_region5 = as(sumExtent_region5, "SpatialPolygons") 
crs(sumShp_region5) = crs(sumPrecip_region5)

# Clip mountain count to the sum precip shapefile
mountCount_region1 = crop(mountCount_region1_full, sumShp_region1)
mountCount_region2 = crop(mountCount_region2_full, sumShp_region2)
mountCount_region3 = crop(mountCount_region3_full, sumShp_region3)
mountCount_region4 = crop(mountCount_region4_full, sumShp_region4)
mountCount_region5 = crop(mountCount_region5_full, sumShp_region5)

5.1 Beregne antall fjellpixler i hver “sumPrecip” rasterpixel for hver region og slå de to lagene sammen til en raster

# Reduce mountain count region 1 (100 x 100m) to sumPrecip region 1 (1000 x 1000m). 
mountCount_region1_red = raster::aggregate(mountCount_region1, fact = 10, fun = sum)
extent(mountCount_region1_red) = extent(sumPrecip_region1)
nSumPrecip_region1 = list()
for(i in 1:nlayers(sumPrecip_region1)){
  nSumPrecip_region1[[i]] = overlay(sumPrecip_region1[[i]], mountCount_region1_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
nSumPrecip_region1 = stack(nSumPrecip_region1)
sumPrecipMountCount_region1 = addLayer(nSumPrecip_region1, mountCount_region1_red) 
names(sumPrecipMountCount_region1) = c("1957":"2020", "mountCount") 

# Reduce mountain count region 2 (100 x 100m) to sumPrecip region 2 (1000 x 1000m). 
mountCount_region2_red = raster::aggregate(mountCount_region2, fact = 10, fun = sum)
extent(mountCount_region2_red) = extent(sumPrecip_region2)
MCount_region2_red = resample(mountCount_region2_red, sumPrecip_region2, method = "ngb") # Had to resample because mountain count had one column too few
nSumPrecip_region2 = list()
for(i in 1:nlayers(sumPrecip_region2)){
  nSumPrecip_region2[[i]] = overlay(sumPrecip_region2[[i]], MCount_region2_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
nSumPrecip_region2 = stack(nSumPrecip_region2)
sumPrecipMountCount_region2 = addLayer(nSumPrecip_region2, MCount_region2_red) 
names(sumPrecipMountCount_region2) = c("1957":"2020", "mountCount") 

# Reduce mountain count region 3 (100 x 100m) to sumPrecip region 3 (1000 x 1000m). 
mountCount_region3_red = raster::aggregate(mountCount_region3, fact = 10, fun = sum)
extent(mountCount_region3_red) = extent(sumPrecip_region3)
MCount_region3_red = resample(mountCount_region3_red, sumPrecip_region3, method = "ngb") # Had to resample because mountain count had one column too few
nSumPrecip_region3 = list()
for(i in 1:nlayers(sumPrecip_region3)){
  nSumPrecip_region3[[i]] = overlay(sumPrecip_region3[[i]], MCount_region3_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
nSumPrecip_region3 = stack(nSumPrecip_region3)
sumPrecipMountCount_region3 = addLayer(nSumPrecip_region3, MCount_region3_red) 
names(sumPrecipMountCount_region3) = c("1957":"2020", "mountCount") 

# Reduce mountain count region 4 (100 x 100m) to sumPrecip region 4 (1000 x 1000m). 
mountCount_region4_red = raster::aggregate(mountCount_region4, fact = 10, fun = sum)
extent(mountCount_region4_red) = extent(sumPrecip_region4)
nSumPrecip_region4 = list()
for(i in 1:nlayers(sumPrecip_region4)){
  nSumPrecip_region4[[i]] = overlay(sumPrecip_region4[[i]], mountCount_region4_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
nSumPrecip_region4 = stack(nSumPrecip_region4)
sumPrecipMountCount_region4 = addLayer(nSumPrecip_region4, mountCount_region4_red) 
names(sumPrecipMountCount_region4) = c("1957":"2020", "mountCount") 

# Reduce mountain count region 5 (100 x 100m) to sumPrecip region 5 (1000 x 1000m). 
mountCount_region5_red = raster::aggregate(mountCount_region5, fact = 10, fun = sum)
extent(mountCount_region5_red) = extent(sumPrecip_region5)
nSumPrecip_region5 = list()
for(i in 1:nlayers(sumPrecip_region5)){
  nSumPrecip_region5[[i]] = overlay(sumPrecip_region5[[i]], mountCount_region5_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
nSumPrecip_region5 = stack(nSumPrecip_region5)
sumPrecipMountCount_region5 = addLayer(nSumPrecip_region5, mountCount_region5_red) 
names(sumPrecipMountCount_region5) = c("1957":"2020", "mountCount") 

# Set raster export directory
RDir_sumPrecip_r1 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/", paste(paste("sumPrecipMountCount_region1", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")
RDir_sumPrecip_r2 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/", paste(paste("sumPrecipMountCount_region2", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")
RDir_sumPrecip_r3 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/", paste(paste("sumPrecipMountCount_region3", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")
RDir_sumPrecip_r4 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/", paste(paste("sumPrecipMountCount_region4", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")
RDir_sumPrecip_r5 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/", paste(paste("sumPrecipMountCount_region5", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

# Export the rasters
writeRaster(sumPrecipMountCount_region1, RDir_sumPrecip_r1, format = "raster", overwrite = TRUE)
writeRaster(sumPrecipMountCount_region2, RDir_sumPrecip_r2, format = "raster", overwrite = TRUE)
writeRaster(sumPrecipMountCount_region3, RDir_sumPrecip_r3, format = "raster", overwrite = TRUE)
writeRaster(sumPrecipMountCount_region4, RDir_sumPrecip_r4, format = "raster", overwrite = TRUE)
writeRaster(sumPrecipMountCount_region5, RDir_sumPrecip_r5, format = "raster", overwrite = TRUE)

5.2 Beregne antall fjellpixler i hver “daysPrecip” rasterpixel for hver region og slå de to lagene sammen til en raster

# Days precipitation region 1 1957 - 2020
daysPrecipMountCount_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/", full.names = TRUE, pattern = "region1.tif$")
daysPrecipMountCount_region1 = stack(daysPrecipMountCount_region1)
extent(mountCount_region1_red) = extent(daysPrecipMountCount_region1)
dDaysPrecip_region1 = list()
for(i in 1:nlayers(daysPrecipMountCount_region1)){
  dDaysPrecip_region1[[i]] = overlay(daysPrecipMountCount_region1[[i]], mountCount_region1_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
dDaysPrecip_region1 = stack(dDaysPrecip_region1)
daysPrecipMountCount_region1 = addLayer(dDaysPrecip_region1, mountCount_region1_red) # use the reduced mountain count created above (5.1)
names(daysPrecipMountCount_region1) = c("1957":"2020", "mountCount") 

# Days precipitation region 2 1957 - 2020
daysPrecipMountCount_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/", full.names = TRUE, pattern = "region2.tif$")
daysPrecipMountCount_region2 = stack(daysPrecipMountCount_region2)
extent(mountCount_region2_red) = extent(daysPrecipMountCount_region2)
dDaysPrecip_region2 = list()
for(i in 1:nlayers(daysPrecipMountCount_region2)){
  dDaysPrecip_region2[[i]] = overlay(daysPrecipMountCount_region2[[i]], MCount_region2_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
dDaysPrecip_region2 = stack(dDaysPrecip_region2)
daysPrecipMountCount_region2 = addLayer(dDaysPrecip_region2, MCount_region2_red) # use the reduced mountain count created above (5.1)
names(daysPrecipMountCount_region2) = c("1957":"2020", "mountCount") 

# Days precipitation region 3 1957 - 2020
daysPrecipMountCount_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/", full.names = TRUE, pattern = "region3.tif$")
daysPrecipMountCount_region3 = stack(daysPrecipMountCount_region3)
extent(mountCount_region3_red) = extent(daysPrecipMountCount_region3)
dDaysPrecip_region3 = list()
for(i in 1:nlayers(daysPrecipMountCount_region3)){
  dDaysPrecip_region3[[i]] = overlay(daysPrecipMountCount_region3[[i]], MCount_region3_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
dDaysPrecip_region3 = stack(dDaysPrecip_region3)
daysPrecipMountCount_region3 = addLayer(dDaysPrecip_region3, MCount_region3_red) # use the reduced mountain count created above (5.1)
names(daysPrecipMountCount_region3) = c("1957":"2020", "mountCount") 

# Days precipitation region 4 1957 - 2020
daysPrecipMountCount_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/", full.names = TRUE, pattern = "region4.tif$")
daysPrecipMountCount_region4 = stack(daysPrecipMountCount_region4)
extent(mountCount_region4_red) = extent(daysPrecipMountCount_region4)
dDaysPrecip_region4 = list()
for(i in 1:nlayers(daysPrecipMountCount_region4)){
  dDaysPrecip_region4[[i]] = overlay(daysPrecipMountCount_region4[[i]], mountCount_region4_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
dDaysPrecip_region4 = stack(dDaysPrecip_region4)
daysPrecipMountCount_region4 = addLayer(dDaysPrecip_region4, mountCount_region4_red) # use the reduced mountain count created above (5.1)
names(daysPrecipMountCount_region4) = c("1957":"2020", "mountCount") 

# Days precipitation region 5 1957 - 2020
daysPrecipMountCount_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/", full.names = TRUE, pattern = "region5.tif$")
daysPrecipMountCount_region5 = stack(daysPrecipMountCount_region5)
extent(mountCount_region5_red) = extent(daysPrecipMountCount_region5)
dDaysPrecip_region5 = list()
for(i in 1:nlayers(daysPrecipMountCount_region5)){
  dDaysPrecip_region5[[i]] = overlay(daysPrecipMountCount_region5[[i]], mountCount_region5_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
dDaysPrecip_region5 = stack(dDaysPrecip_region5)
daysPrecipMountCount_region5 = addLayer(dDaysPrecip_region5, mountCount_region5_red) # use the reduced mountain count created above (5.1)
names(daysPrecipMountCount_region5) = c("1957":"2020", "mountCount") 


# Export the raster
RasDir_r1 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/", paste(paste("daysPrecipMountCount_region1", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasDir_r2 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/", paste(paste("daysPrecipMountCount_region2", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasDir_r3 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/", paste(paste("daysPrecipMountCount_region3", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasDir_r4 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/", paste(paste("daysPrecipMountCount_region4", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasDir_r5 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/", paste(paste("daysPrecipMountCount_region5", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

# Write Raster
writeRaster(daysPrecipMountCount_region1, RasDir_r1, format = "raster", overwrite = TRUE)
writeRaster(daysPrecipMountCount_region2, RasDir_r2, format = "raster", overwrite = TRUE)
writeRaster(daysPrecipMountCount_region3, RasDir_r3, format = "raster", overwrite = TRUE)
writeRaster(daysPrecipMountCount_region4, RasDir_r4, format = "raster", overwrite = TRUE)
writeRaster(daysPrecipMountCount_region5, RasDir_r5, format = "raster", overwrite = TRUE)

5.3 Beregne antall fjellpixler i hver “meanSummer” rasterpixel for hver region og slå de to lagene sammen til en raster

# Mean summer region 1 1957 - 2020
meanSummerMountCount_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/", full.names = TRUE, pattern = "region1.tif$")
meanSummerMountCount_region1 = stack(meanSummerMountCount_region1)
extent(mountCount_region1_red) = extent(meanSummerMountCount_region1)
mSummer_region1 = list()
for(i in 1:nlayers(meanSummerMountCount_region1)){
  mSummer_region1[[i]] = overlay(meanSummerMountCount_region1[[i]], mountCount_region1_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
mSummer_region1 = stack(mSummer_region1)
meanSummerMountCount_region1 = addLayer(mSummer_region1, mountCount_region1_red) # use the reduced mountain count created above (5.1)
names(meanSummerMountCount_region1) = c("1957":"2020", "mountCount") 

# Mean summer region 2 1957 - 2020
meanSummerMountCount_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/", full.names = TRUE, pattern = "region2.tif$")
meanSummerMountCount_region2 = stack(meanSummerMountCount_region2)
extent(mountCount_region2_red) = extent(meanSummerMountCount_region2)
mSummer_region2 = list()
for(i in 1:nlayers(meanSummerMountCount_region2)){
  mSummer_region2[[i]] = overlay(meanSummerMountCount_region2[[i]], MCount_region2_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
mSummer_region2 = stack(mSummer_region2)
meanSummerMountCount_region2 = addLayer(mSummer_region2, MCount_region2_red) # use the reduced mountain count created above (5.1)
names(meanSummerMountCount_region2) = c("1957":"2020", "mountCount") 

# Mean summer region 3 1957 - 2020
meanSummerMountCount_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/", full.names = TRUE, pattern = "region3.tif$")
meanSummerMountCount_region3 = stack(meanSummerMountCount_region3)
extent(mountCount_region3_red) = extent(meanSummerMountCount_region3)
mSummer_region3 = list()
for(i in 1:nlayers(meanSummerMountCount_region3)){
  mSummer_region3[[i]] = overlay(meanSummerMountCount_region3[[i]], MCount_region3_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
mSummer_region3 = stack(mSummer_region3)
meanSummerMountCount_region3 = addLayer(mSummer_region3, MCount_region3_red) # use the reduced mountain count created above (5.1)
names(meanSummerMountCount_region3) = c("1957":"2020", "mountCount") 

# Mean summer region 4 1957 - 2020
meanSummerMountCount_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/", full.names = TRUE, pattern = "region4.tif$")
meanSummerMountCount_region4 = stack(meanSummerMountCount_region4)
extent(mountCount_region4_red) = extent(meanSummerMountCount_region4)
mSummer_region4 = list()
for(i in 1:nlayers(meanSummerMountCount_region4)){
  mSummer_region4[[i]] = overlay(meanSummerMountCount_region4[[i]], mountCount_region4_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
mSummer_region4 = stack(mSummer_region4)
meanSummerMountCount_region4 = addLayer(mSummer_region4, mountCount_region4_red) # use the reduced mountain count created above (5.1)
names(meanSummerMountCount_region4) = c("1957":"2020", "mountCount") 

# Mean summer region 5 1957 - 2020
meanSummerMountCount_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/", full.names = TRUE, pattern = "region5.tif$")
meanSummerMountCount_region5 = stack(meanSummerMountCount_region5)
extent(mountCount_region5_red) = extent(meanSummerMountCount_region5)
mSummer_region5 = list()
for(i in 1:nlayers(meanSummerMountCount_region5)){
  mSummer_region5[[i]] = overlay(meanSummerMountCount_region5[[i]], mountCount_region5_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
mSummer_region5 = stack(mSummer_region5)
meanSummerMountCount_region5 = addLayer(mSummer_region5, mountCount_region5_red) # use the reduced mountain count created above (5.1)
names(meanSummerMountCount_region5) = c("1957":"2020", "mountCount") 


# Export the raster
RasDir_r1 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/", paste(paste("meanSummerMountCount_region1", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasDir_r2 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/", paste(paste("meanSummerMountCount_region2", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasDir_r3 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/", paste(paste("meanSummerMountCount_region3", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasDir_r4 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/", paste(paste("meanSummerMountCount_region4", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasDir_r5 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/", paste(paste("meanSummerMountCount_region5", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

# Write Raster
writeRaster(meanSummerMountCount_region1, RasDir_r1, format = "raster", overwrite = TRUE)
writeRaster(meanSummerMountCount_region2, RasDir_r2, format = "raster", overwrite = TRUE)
writeRaster(meanSummerMountCount_region3, RasDir_r3, format = "raster", overwrite = TRUE)
writeRaster(meanSummerMountCount_region4, RasDir_r4, format = "raster", overwrite = TRUE)
writeRaster(meanSummerMountCount_region5, RasDir_r5, format = "raster", overwrite = TRUE)

5.4 Beregne antall fjellpixler i hver “meanWinter” rasterpixel for hver region og slå de to lagene sammen til en raster

# Mean Winter temperature region 1
meanWinter_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/", full.names = TRUE, pattern = "region1.tif$")
meanWinter_region1 = stack(meanWinter_region1) 
mWinter_region1 = list()
for(i in 1:nlayers(meanWinter_region1)){
  mWinter_region1[[i]] = overlay(meanWinter_region1[[i]], mountCount_region1_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
mWinter_region1 = stack(mWinter_region1)
meanWinterMountCount_region1 = addLayer(mWinter_region1, mountCount_region1_red) # use the reduced mountain count created above (5.1)
winterSeason = c("1957":"2020")
wS = c()
for(i in 2:length(winterSeason)){
  wS = append(wS, paste(substr(x = winterSeason[i-1], start = 3, stop = 4), substr(x = winterSeason[i], start = 3, stop = 4), sep = "-"))
}
names(meanWinterMountCount_region1) = c(wS, "mountCount") # rename the layers

# Mean Winter temperature region 2
meanWinter_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/", full.names = TRUE, pattern = "region2.tif$")
meanWinter_region2 = stack(meanWinter_region2) 
mWinter_region2 = list()
for(i in 1:nlayers(meanWinter_region2)){
  mWinter_region2[[i]] = overlay(meanWinter_region2[[i]], MCount_region2_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
mWinter_region2 = stack(mWinter_region2)
meanWinterMountCount_region2 = addLayer(mWinter_region2, MCount_region2_red)
names(meanWinterMountCount_region2) = c(wS, "mountCount")

# Mean Winter temperature region 3
meanWinter_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/", full.names = TRUE, pattern = "region3.tif$")
meanWinter_region3 = stack(meanWinter_region3) 
mWinter_region3 = list()
for(i in 1:nlayers(meanWinter_region3)){
  mWinter_region3[[i]] = overlay(meanWinter_region3[[i]], MCount_region3_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
mWinter_region3 = stack(mWinter_region3)
meanWinterMountCount_region3 = addLayer(mWinter_region3, MCount_region3_red)
names(meanWinterMountCount_region3) = c(wS, "mountCount")

# Mean Winter temperature region 4
meanWinter_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/", full.names = TRUE, pattern = "region4.tif$")
meanWinter_region4 = stack(meanWinter_region4)
mWinter_region4 = list()
for(i in 1:nlayers(meanWinter_region4)){
  mWinter_region4[[i]] = overlay(meanWinter_region4[[i]], mountCount_region4_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
mWinter_region4 = stack(mWinter_region4)
meanWinterMountCount_region4 = addLayer(mWinter_region4, mountCount_region4_red)
names(meanWinterMountCount_region4) = c(wS, "mountCount")

# Mean Winter temperature region 5
meanWinter_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/", full.names = TRUE, pattern = "region5.tif$")
meanWinter_region5 = stack(meanWinter_region5) 
mWinter_region5 = list()
for(i in 1:nlayers(meanWinter_region5)){
  mWinter_region5[[i]] = overlay(meanWinter_region5[[i]], mountCount_region5_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
mWinter_region5 = stack(mWinter_region5)
meanWinterMountCount_region5 = addLayer(mWinter_region5, mountCount_region5_red)
names(meanWinterMountCount_region5) = c(wS, "mountCount")

# Export the raster
RasterLocation_r1 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/", paste(paste("meanWinterMountCount_region1", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r2 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/", paste(paste("meanWinterMountCount_region2", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r3 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/", paste(paste("meanWinterMountCount_region3", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r4 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/", paste(paste("meanWinterMountCount_region4", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r5 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/", paste(paste("meanWinterMountCount_region5", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

# Write Raster
writeRaster(meanWinterMountCount_region1, RasterLocation_r1, format = "raster", overwrite = TRUE)
writeRaster(meanWinterMountCount_region2, RasterLocation_r2, format = "raster", overwrite = TRUE)
writeRaster(meanWinterMountCount_region3, RasterLocation_r3, format = "raster", overwrite = TRUE)
writeRaster(meanWinterMountCount_region4, RasterLocation_r4, format = "raster", overwrite = TRUE)
writeRaster(meanWinterMountCount_region5, RasterLocation_r5, format = "raster", overwrite = TRUE)

5.5 Beregne antall fjellpixler i hver “daysSnowCover” rasterpixel for hver region og slå de to lagene sammen til en raster

# Days snow cover region 1
daysSnowCover_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", full.names = TRUE, pattern = "region1.tif$")
daysSnowCover_region1 = stack(daysSnowCover_region1)
daysSC_region1 = list()
for(i in 1:nlayers(daysSnowCover_region1)){
  daysSC_region1[[i]] = overlay(daysSnowCover_region1[[i]], mountCount_region1_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
daysSC_region1 = stack(daysSC_region1)
daysSnowCoverMountCount_region1 = addLayer(daysSC_region1, mountCount_region1_red) 
winterSeason = c("1957":"2020")
wS = c()
for(i in 2:length(winterSeason)){
  wS = append(wS, paste(substr(x = winterSeason[i-1], start = 3, stop = 4), substr(x = winterSeason[i], start = 3, stop = 4), sep = "-"))
}
names(daysSnowCoverMountCount_region1) = c(wS, "mountCount") # rename the layers

# Days snow cover region 2
daysSnowCover_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", full.names = TRUE, pattern = "region2.tif$")
daysSnowCover_region2 = stack(daysSnowCover_region2)
daysSC_region2 = list()
for(i in 1:nlayers(daysSnowCover_region2)){
  daysSC_region2[[i]] = overlay(daysSnowCover_region2[[i]], MCount_region2_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
daysSC_region2 = stack(daysSC_region2)
daysSnowCoverMountCount_region2 = addLayer(daysSC_region2, MCount_region2_red) 
names(daysSnowCoverMountCount_region2) = c(wS, "mountCount") 

# Days snow cover region 3
daysSnowCover_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", full.names = TRUE, pattern = "region3.tif$")
daysSnowCover_region3 = stack(daysSnowCover_region3) 
daysSC_region3 = list()
for(i in 1:nlayers(daysSnowCover_region3)){
  daysSC_region3[[i]] = overlay(daysSnowCover_region3[[i]], MCount_region3_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
daysSC_region3 = stack(daysSC_region3)
daysSnowCoverMountCount_region3 = addLayer(daysSC_region3, MCount_region3_red) 
names(daysSnowCoverMountCount_region3) = c(wS, "mountCount") 


# Days snow cover region 4
daysSnowCover_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", full.names = TRUE, pattern = "region4.tif$")
daysSnowCover_region4 = stack(daysSnowCover_region4) 
daysSC_region4 = list()
for(i in 1:nlayers(daysSnowCover_region4)){
  daysSC_region4[[i]] = overlay(daysSnowCover_region4[[i]], mountCount_region4_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
daysSC_region4 = stack(daysSC_region4)
daysSnowCoverMountCount_region4 = addLayer(daysSC_region4, mountCount_region4_red) 
names(daysSnowCoverMountCount_region4) = c(wS, "mountCount") 

# Days snow cover region 5
daysSnowCover_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", full.names = TRUE, pattern = "region5.tif$")
daysSnowCover_region5 = stack(daysSnowCover_region5) 
daysSC_region5 = list()
for(i in 1:nlayers(daysSnowCover_region5)){
  daysSC_region5[[i]] = overlay(daysSnowCover_region5[[i]], mountCount_region5_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
daysSC_region5 = stack(daysSC_region5)
daysSnowCoverMountCount_region5 = addLayer(daysSC_region5, mountCount_region5_red) 
names(daysSnowCoverMountCount_region5) = c(wS, "mountCount") 


# Export the raster
RasterLocation_r1 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", paste(paste("daysSnowCoverMountCount_region1", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r2 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", paste(paste("daysSnowCoverMountCount_region2", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r3 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", paste(paste("daysSnowCoverMountCount_region3", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r4 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", paste(paste("daysSnowCoverMountCount_region4", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r5 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", paste(paste("daysSnowCoverMountCount_region5", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

# Write Raster
writeRaster(daysSnowCoverMountCount_region1, RasterLocation_r1, format = "raster", overwrite = TRUE)
writeRaster(daysSnowCoverMountCount_region2, RasterLocation_r2, format = "raster", overwrite = TRUE)
writeRaster(daysSnowCoverMountCount_region3, RasterLocation_r3, format = "raster", overwrite = TRUE)
writeRaster(daysSnowCoverMountCount_region4, RasterLocation_r4, format = "raster", overwrite = TRUE)
writeRaster(daysSnowCoverMountCount_region5, RasterLocation_r5, format = "raster", overwrite = TRUE)

5.6 Beregne antall fjellpixler i hver “growthSeason” rasterpixel for hver region og slå de to lagene sammen til en raster

# Growth season region 1
growthSeason_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", full.names = TRUE, pattern = "region1.tif$")
growthSeason_region1 = stack(growthSeason_region1) 
growthS_region1 = list()
for(i in 1:nlayers(growthSeason_region1)){
  growthS_region1[[i]] = overlay(growthSeason_region1[[i]], mountCount_region1_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
growthS_region1 = stack(growthS_region1)
growthSeasonMountCount_region1 = addLayer(growthS_region1, mountCount_region1_red) 
names(growthSeasonMountCount_region1) = c("1957":"2020", "forCount") 

# Growth season region 2
growthSeason_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", full.names = TRUE, pattern = "region2.tif$")
growthSeason_region2 = stack(growthSeason_region2) 
growthS_region2 = list()
for(i in 1:nlayers(growthSeason_region2)){
  growthS_region2[[i]] = overlay(growthSeason_region2[[i]], MCount_region2_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
growthS_region2 = stack(growthS_region2)
growthSeasonMountCount_region2 = addLayer(growthS_region2, MCount_region2_red) 
names(growthSeasonMountCount_region2) = c("1957":"2020", "forCount") 

# Growth season region 3
growthSeason_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", full.names = TRUE, pattern = "region3.tif$")
growthSeason_region3 = stack(growthSeason_region3) 
growthS_region3 = list()
for(i in 1:nlayers(growthSeason_region3)){
  growthS_region3[[i]] = overlay(growthSeason_region3[[i]], MCount_region3_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
growthS_region3 = stack(growthS_region3)
growthSeasonMountCount_region3 = addLayer(growthS_region3, MCount_region3_red) 
names(growthSeasonMountCount_region3) = c("1957":"2020", "forCount") 

# Growth season region 4
growthSeason_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", full.names = TRUE, pattern = "region4.tif$")
growthSeason_region4 = stack(growthSeason_region4) 
growthS_region4 = list()
for(i in 1:nlayers(growthSeason_region4)){
  growthS_region4[[i]] = overlay(growthSeason_region4[[i]], mountCount_region4_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
growthS_region4 = stack(growthS_region4)
growthSeasonMountCount_region4 = addLayer(growthS_region4, mountCount_region4_red) 
names(growthSeasonMountCount_region4) = c("1957":"2020", "forCount") 

# Growth season region 5
growthSeason_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", full.names = TRUE, pattern = "region5.tif$")
growthSeason_region5 = stack(growthSeason_region5) 
growthS_region5 = list()
for(i in 1:nlayers(growthSeason_region5)){
  growthS_region5[[i]] = overlay(growthSeason_region5[[i]], mountCount_region5_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
growthS_region5 = stack(growthS_region5)
growthSeasonMountCount_region5 = addLayer(growthS_region5, mountCount_region5_red) 
names(growthSeasonMountCount_region5) = c("1957":"2020", "forCount") 


# Export the raster
RasterLocation_r1 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", paste(paste("growthSeasonMountCount_region1", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r2 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", paste(paste("growthSeasonMountCount_region2", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r3 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", paste(paste("growthSeasonMountCount_region3", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r4 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", paste(paste("growthSeasonMountCount_region4", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r5 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", paste(paste("growthSeasonMountCount_region5", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")


# NB - The first year 1957 only contain the months Sept - Dec and should not be used in analyses

# Write Raster
writeRaster(growthSeasonMountCount_region1, RasterLocation_r1, format = "raster", overwrite = TRUE)
writeRaster(growthSeasonMountCount_region2, RasterLocation_r2, format = "raster", overwrite = TRUE)
writeRaster(growthSeasonMountCount_region3, RasterLocation_r3, format = "raster", overwrite = TRUE)
writeRaster(growthSeasonMountCount_region4, RasterLocation_r4, format = "raster", overwrite = TRUE)
writeRaster(growthSeasonMountCount_region5, RasterLocation_r5, format = "raster", overwrite = TRUE)

5.7 Beregne antall fjellpixler i hver “winterRain” rasterpixel for hver region og slå de to lagene sammen til en raster

# Winter rain region 1
winterRain_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", full.names = TRUE, pattern = "region1.tif$")
winterRain_region1 = stack(winterRain_region1) 
winterR_region1 = list()
for(i in 1:nlayers(winterRain_region1)){
  winterR_region1[[i]] = overlay(winterRain_region1[[i]], mountCount_region1_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
winterR_region1 = stack(winterR_region1)
winterRainMountCount_region1 = addLayer(winterR_region1, mountCount_region1_red) 
names(winterRainMountCount_region1) = c("1957":"2020", "mountCount") 

# Winter rain region 2
winterRain_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", full.names = TRUE, pattern = "region2.tif$")
winterRain_region2 = stack(winterRain_region2) 
winterR_region2 = list()
for(i in 1:nlayers(winterRain_region2)){
  winterR_region2[[i]] = overlay(winterRain_region2[[i]], MCount_region2_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
winterR_region2 = stack(winterR_region2)
winterRainMountCount_region2 = addLayer(winterR_region2, MCount_region2_red) 
names(winterRainMountCount_region2) = c("1957":"2020", "mountCount") 

# Winter rain region 3
winterRain_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", full.names = TRUE, pattern = "region3.tif$")
winterRain_region3 = stack(winterRain_region3) 
winterR_region3 = list()
for(i in 1:nlayers(winterRain_region3)){
  winterR_region3[[i]] = overlay(winterRain_region3[[i]], MCount_region3_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
winterR_region3 = stack(winterR_region3)
winterRainMountCount_region3 = addLayer(winterR_region3, MCount_region3_red) 
names(winterRainMountCount_region3) = c("1957":"2020", "mountCount")

# Winter rain region 4
winterRain_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", full.names = TRUE, pattern = "region4.tif$")
winterRain_region4 = stack(winterRain_region4) 
winterR_region4 = list()
for(i in 1:nlayers(winterRain_region4)){
  winterR_region4[[i]] = overlay(winterRain_region4[[i]], mountCount_region4_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
winterR_region4 = stack(winterR_region4)
winterRainMountCount_region4 = addLayer(winterR_region4, mountCount_region4_red) 
names(winterRainMountCount_region4) = c("1957":"2020", "mountCount")

# Winter rain region 5
winterRain_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", full.names = TRUE, pattern = "region5.tif$")
winterRain_region5 = stack(winterRain_region5) 
winterR_region5 = list()
for(i in 1:nlayers(winterRain_region5)){
  winterR_region5[[i]] = overlay(winterRain_region5[[i]], mountCount_region5_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
winterR_region5 = stack(winterR_region5)
winterRainMountCount_region5 = addLayer(winterR_region5, mountCount_region5_red) 
names(winterRainMountCount_region5) = c("1957":"2020", "mountCount")


# Export the raster
RasterLocation_r1 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", paste(paste("winterRainMountCount_region1", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r2 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", paste(paste("winterRainMountCount_region2", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r3 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", paste(paste("winterRainMountCount_region3", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r4 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", paste(paste("winterRainMountCount_region4", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r5 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", paste(paste("winterRainMountCount_region5", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")


# Write Raster
writeRaster(winterRainMountCount_region1, RasterLocation_r1, format = "raster", overwrite = TRUE)
writeRaster(winterRainMountCount_region2, RasterLocation_r2, format = "raster", overwrite = TRUE)
writeRaster(winterRainMountCount_region3, RasterLocation_r3, format = "raster", overwrite = TRUE)
writeRaster(winterRainMountCount_region4, RasterLocation_r4, format = "raster", overwrite = TRUE)
writeRaster(winterRainMountCount_region5, RasterLocation_r5, format = "raster", overwrite = TRUE)

5.8 Beregne antall fjellpixler i hver “snowDepth” rasterpixel for hver region og slå de to lagene sammen til en raster

# Snow depth region 1
snowDepth_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/", full.names = TRUE, pattern = "region1.tif$")
snowDepth_region1 = stack(snowDepth_region1) 
snowD_region1 = list()
for(i in 1:nlayers(snowDepth_region1)){
  snowD_region1[[i]] = overlay(snowDepth_region1[[i]], mountCount_region1_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
snowD_region1 = stack(snowD_region1)
snowDepthMountCount_region1 = addLayer(snowD_region1, mountCount_region1_red)
winterSeason = c(1957:2020)
wS = c()
for(i in 2:length(winterSeason)){
  wS = append(wS, paste(substr(x = winterSeason[i-1], start = 3, stop = 4), substr(x = winterSeason[i], start = 3, stop = 4), sep = "-"))
}
names(snowDepthMountCount_region1) = c(wS, "mountCount") 

# Winter rain region 2
snowDepth_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/", full.names = TRUE, pattern = "region2.tif$")
snowDepth_region2 = stack(snowDepth_region2) 
snowD_region2 = list()
for(i in 1:nlayers(snowDepth_region2)){
  snowD_region2[[i]] = overlay(snowDepth_region2[[i]], MCount_region2_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
snowD_region2 = stack(snowD_region2)
snowDepthMountCount_region2 = addLayer(snowD_region2, MCount_region2_red) 
names(snowDepthMountCount_region2) = c(wS, "mountCount") 

# Winter rain region 3
snowDepth_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/", full.names = TRUE, pattern = "region3.tif$")
snowDepth_region3 = stack(snowDepth_region3) 
snowD_region3 = list()
for(i in 1:nlayers(snowDepth_region3)){
  snowD_region3[[i]] = overlay(snowDepth_region3[[i]], MCount_region3_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
snowD_region3 = stack(snowD_region3)
snowDepthMountCount_region3 = addLayer(snowD_region3, MCount_region3_red) 
names(snowDepthMountCount_region3) = c(wS, "mountCount")

# Winter rain region 4
snowDepth_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/", full.names = TRUE, pattern = "region4.tif$")
snowDepth_region4 = stack(snowDepth_region4) 
snowD_region4 = list()
for(i in 1:nlayers(snowDepth_region4)){
  snowD_region4[[i]] = overlay(snowDepth_region4[[i]], mountCount_region4_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
snowD_region4 = stack(snowD_region4)
snowDepthMountCount_region4 = addLayer(snowD_region4, mountCount_region4_red) 
names(snowDepthMountCount_region4) = c(wS, "mountCount")

# Winter rain region 5
snowDepth_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/", full.names = TRUE, pattern = "region5.tif$")
snowDepth_region5 = stack(snowDepth_region5) 
snowD_region5 = list()
for(i in 1:nlayers(snowDepth_region5)){
  snowD_region5[[i]] = overlay(snowDepth_region5[[i]], mountCount_region5_red, fun = function(x, y){
    x[y <= 21000] = NA; x # Use only cells that consist of more than 50% mountain
  })
}
snowD_region5 = stack(snowD_region5)
snowDepthMountCount_region5 = addLayer(snowD_region5, mountCount_region5_red) 
names(snowDepthMountCount_region5) = c(wS, "mountCount")


# Export the raster
RasterLocation_r1 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/", paste(paste("snowDepthMountCount_region1", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r2 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/", paste(paste("snowDepthMountCount_region2", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r3 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/", paste(paste("snowDepthMountCount_region3", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r4 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/", paste(paste("snowDepthMountCount_region4", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")

RasterLocation_r5 = paste("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/", paste(paste("snowDepthMountCount_region5", 1957, 2020, sep = "_"), ".grd", sep = ""), sep = "")


# Write Raster
writeRaster(snowDepthMountCount_region1, RasterLocation_r1, format = "raster", overwrite = TRUE)
writeRaster(snowDepthMountCount_region2, RasterLocation_r2, format = "raster", overwrite = TRUE)
writeRaster(snowDepthMountCount_region3, RasterLocation_r3, format = "raster", overwrite = TRUE)
writeRaster(snowDepthMountCount_region4, RasterLocation_r4, format = "raster", overwrite = TRUE)
writeRaster(snowDepthMountCount_region5, RasterLocation_r5, format = "raster", overwrite = TRUE)

6 Medianverdi

6.1 Beregne medianen for hvert år for sum nedbør

# Region 1
sumPrecipMountCount_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/", full.names = TRUE, pattern = "region1_1957_2020.grd")[1]
sumPrecipMountCount_region1 = stack(sumPrecipMountCount_region1)
years = c(1:64)
med = c()
for(i in years){
  temp_med = median(values(sumPrecipMountCount_region1[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
sumPrecip_val = data.frame(matrix(c(rep("Nord-Norge", 64), rep("nedbør", 64), rep("mm", 64), 1957:2020, med), nrow = 64, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(sumPrecip_val) = columnNames

# Region 2
sumPrecipMountCount_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/", full.names = TRUE, pattern = "region2_1957_2020.grd")[1]
sumPrecipMountCount_region2 = stack(sumPrecipMountCount_region2)
years = c(1:64)
med = c()
for(i in years){
  temp_med = median(values(sumPrecipMountCount_region2[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
sumPrecip_val_region2 = data.frame(matrix(c(rep("Midt-Norge", 64), rep("nedbør", 64), rep("mm", 64), 1957:2020, med), nrow = 64, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(sumPrecip_val_region2) = columnNames
sumPrecip_val = rbind(sumPrecip_val, sumPrecip_val_region2)

# Region 3
sumPrecipMountCount_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/", full.names = TRUE, pattern = "region3_1957_2020.grd")[1]
sumPrecipMountCount_region3 = stack(sumPrecipMountCount_region3)
years = c(1:64)
med = c()
for(i in years){
  temp_med = median(values(sumPrecipMountCount_region3[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
sumPrecip_val_region3 = data.frame(matrix(c(rep("Østlandet", 64), rep("nedbør", 64), rep("mm", 64), 1957:2020, med), nrow = 64, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(sumPrecip_val_region3) = columnNames
sumPrecip_val = rbind(sumPrecip_val, sumPrecip_val_region3)

# Region 4
sumPrecipMountCount_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/", full.names = TRUE, pattern = "region4_1957_2020.grd")[1]
sumPrecipMountCount_region4 = stack(sumPrecipMountCount_region4)
years = c(1:64)
med = c()
for(i in years){
  temp_med = median(values(sumPrecipMountCount_region4[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
sumPrecip_val_region4 = data.frame(matrix(c(rep("Vestlandet", 64), rep("nedbør", 64), rep("mm", 64), 1957:2020, med), nrow = 64, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(sumPrecip_val_region4) = columnNames
sumPrecip_val = rbind(sumPrecip_val, sumPrecip_val_region4)

# Region 5
sumPrecipMountCount_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/", full.names = TRUE, pattern = "region5_1957_2020.grd")[1]
sumPrecipMountCount_region5 = stack(sumPrecipMountCount_region5)
years = c(1:64)
med = c()
for(i in years){
  temp_med = median(values(sumPrecipMountCount_region5[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
sumPrecip_val_region5 = data.frame(matrix(c(rep("Sørlandet", 64), rep("nedbør", 64), rep("mm", 64), 1957:2020, med), nrow = 64, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(sumPrecip_val_region5) = columnNames
sumPrecip_val = rbind(sumPrecip_val, sumPrecip_val_region5)

write_xlsx(sumPrecip_val, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/medianNedbør.xlsx")

# Region 1
# normalSP_r1 = sumPrecipMountCount_region1[[5:34]]
# normalSP_r1_med = median(values(normalSP_r1), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(sumPrecipMountCount_region1[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalSP_r1_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nSP_r1_SD_plus = normalSP_r1_med + (2*sd(med[5:34]))
# nSP_r1_SD_min = normalSP_r1_med - (2*sd(med[5:34]))
# sumPrecip_med = data.frame(matrix(c("norge-norge", "sumPrecip", normalSP_r1_med, nSP_r1_SD_min, nSP_r1_SD_plus), nrow = 1))
# colnames(sumPrecip_med) = c("area", "variable", "norm_med", "norm_-2SD", "norm_+2SD")
# 
# sumPrecip_diff = data.frame(c(1957:2020),matrix(rep(NA, 64*5), nrow = 64, ncol = 5))
# columnNames = c("year","nord-norge", "midt-norge", "østlandet", "vestlandet", "sørlandet")
# colnames(sumPrecip_diff) = columnNames
# sumPrecip_diff[, 2] = diff

# Region 2
# normalSP_r2 = sumPrecipMountCount_region2[[5:34]]
# normalSP_r2_med = median(values(normalSP_r2), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(sumPrecipMountCount_region2[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalSP_r2_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nSP_r2_SD_plus = normalSP_r2_med + (2*sd(med[5:34]))
# nSP_r2_SD_min = normalSP_r2_med - (2*sd(med[5:34]))
# 
# sumPrecip_diff[, 3] = diff
# sumPrecip_med = rbind(sumPrecip_med, c("midt-norge", "sumPrecip", normalSP_r2_med, nSP_r2_SD_min, nSP_r2_SD_plus))
# 
# # Region 3
# normalSP_r3 = sumPrecipMountCount_region3[[5:34]]
# normalSP_r3_med = median(values(normalSP_r3), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(sumPrecipMountCount_region3[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalSP_r3_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nSP_r3_SD_plus = normalSP_r3_med + (2*sd(med[5:34]))
# nSP_r3_SD_min = normalSP_r3_med - (2*sd(med[5:34]))
# 
# sumPrecip_diff[, 4] = diff
# sumPrecip_med = rbind(sumPrecip_med, c("østlandet", "sumPrecip", normalSP_r3_med, nSP_r3_SD_min, nSP_r3_SD_plus))

# Region 4
# normalSP_r4 = sumPrecipMountCount_region4[[5:34]]
# normalSP_r4_med = median(values(normalSP_r4), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(sumPrecipMountCount_region4[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalSP_r4_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nSP_r4_SD_plus = normalSP_r4_med + (2*sd(med[5:34]))
# nSP_r4_SD_min = normalSP_r4_med - (2*sd(med[5:34]))
# 
# sumPrecip_diff[, 5] = diff
# sumPrecip_med = rbind(sumPrecip_med, c("vestlandet", "sumPrecip", normalSP_r4_med, nSP_r4_SD_min, nSP_r4_SD_plus))
# 
# # Region 5
# normalSP_r5 = sumPrecipMountCount_region5[[5:34]]
# normalSP_r5_med = median(values(normalSP_r5), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(sumPrecipMountCount_region5[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalSP_r5_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nSP_r5_SD_plus = normalSP_r5_med + (2*sd(med[5:34]))
# nSP_r5_SD_min = normalSP_r5_med - (2*sd(med[5:34]))
# 
# sumPrecip_diff[, 6] = diff
# sumPrecip_med = rbind(sumPrecip_med, c("sørlandet", "sumPrecip", normalSP_r5_med, nSP_r5_SD_min, nSP_r5_SD_plus))
# 
# sumPrecip_diff$variable = "sumPrecip"
#write_xlsx(sumPrecip_diff, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/sumPrecip_diff.xlsx")
#write_xlsx(sumPrecip_med, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Sum nedbør/sumPrecip_med.xlsx")

6.2 Beregne medianen for hvert år for dager med nedbør

# Region 1
daysPrecipMountCount_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/", full.names = TRUE, pattern = "region1_1957_2020.grd")[1]
daysPrecipMountCount_region1 = stack(daysPrecipMountCount_region1)
years = c(1:64)
med = c()
for(i in years){
  temp_med = median(values(daysPrecipMountCount_region1[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
daysPrecip_val = data.frame(matrix(c(rep("Nord-Norge", 64), rep("nedbør", 64), rep("dager", 64), 1957:2020, med), nrow = 64, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(daysPrecip_val) = columnNames

# Region 2
daysPrecipMountCount_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/", full.names = TRUE, pattern = "region2_1957_2020.grd")[1]
daysPrecipMountCount_region2 = stack(daysPrecipMountCount_region2)
years = c(1:64)
med = c()
for(i in years){
  temp_med = median(values(daysPrecipMountCount_region2[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
daysPrecip_val_region2 = data.frame(matrix(c(rep("Midt-Norge", 64), rep("nedbør", 64), rep("dager", 64), 1957:2020, med), nrow = 64, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(daysPrecip_val_region2) = columnNames
daysPrecip_val = rbind(daysPrecip_val, daysPrecip_val_region2)

# Region 3
daysPrecipMountCount_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/", full.names = TRUE, pattern = "region3_1957_2020.grd")[1]
daysPrecipMountCount_region3 = stack(daysPrecipMountCount_region3)
years = c(1:64)
med = c()
for(i in years){
  temp_med = median(values(daysPrecipMountCount_region3[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
daysPrecip_val_region3 = data.frame(matrix(c(rep("Østlandet", 64), rep("nedbør", 64), rep("dager", 64), 1957:2020, med), nrow = 64, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(daysPrecip_val_region3) = columnNames
daysPrecip_val = rbind(daysPrecip_val, daysPrecip_val_region3)

# Region 4
daysPrecipMountCount_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/", full.names = TRUE, pattern = "region4_1957_2020.grd")[1]
daysPrecipMountCount_region4 = stack(daysPrecipMountCount_region4)
years = c(1:64)
med = c()
for(i in years){
  temp_med = median(values(daysPrecipMountCount_region4[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
daysPrecip_val_region4 = data.frame(matrix(c(rep("Vestlandet", 64), rep("nedbør", 64), rep("dager", 64), 1957:2020, med), nrow = 64, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(daysPrecip_val_region4) = columnNames
daysPrecip_val = rbind(daysPrecip_val, daysPrecip_val_region4)

# Region 5
daysPrecipMountCount_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/", full.names = TRUE, pattern = "region5_1957_2020.grd")[1]
daysPrecipMountCount_region5 = stack(daysPrecipMountCount_region5)
years = c(1:64)
med = c()
for(i in years){
  temp_med = median(values(daysPrecipMountCount_region5[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
daysPrecip_val_region5 = data.frame(matrix(c(rep("Sørlandet", 64), rep("nedbør", 64), rep("dager", 64), 1957:2020, med), nrow = 64, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(daysPrecip_val_region5) = columnNames
daysPrecip_val = rbind(daysPrecip_val, daysPrecip_val_region5)

write_xlsx(daysPrecip_val, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/medianNedbør_dager.xlsx")
# 
# # Region 1
# normalDP_r1 = daysPrecipMountCount_region1[[5:34]]
# normalDP_r1_med = median(values(normalDP_r1), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(daysPrecipMountCount_region1[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalDP_r1_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nDP_r1_SD_plus = normalDP_r1_med + (2*sd(med[5:34]))
# nDP_r1_SD_min = normalDP_r1_med - (2*sd(med[5:34]))
# daysPrecip_med = data.frame(matrix(c("nord-norge", "daysPrecip", normalDP_r1_med, nDP_r1_SD_min, nDP_r1_SD_plus), nrow = 1))
# colnames(daysPrecip_med) = c("area", "variable", "norm_med", "norm_-2SD", "norm_+2SD")
# 
# daysPrecip_diff = data.frame(c(1957:2020),matrix(rep(NA, 64*5), nrow = 64, ncol = 5))
# columnNames = c("year","nord-norge", "midt-norge", "østlandet", "vestlandet", "sørlandet")
# colnames(daysPrecip_diff) = columnNames
# daysPrecip_diff[, 2] = diff
# 
# # Region 2
# normalDP_r2 = daysPrecipMountCount_region2[[5:34]]
# normalDP_r2_med = median(values(normalDP_r2), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(daysPrecipMountCount_region2[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalDP_r2_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nDP_r2_SD_plus = normalDP_r2_med + (2*sd(med[5:34]))
# nDP_r2_SD_min = normalDP_r2_med - (2*sd(med[5:34]))
# 
# daysPrecip_diff[, 3] = diff
# daysPrecip_med = rbind(daysPrecip_med, c("midt-norge", "daysPrecip", normalDP_r2_med, nDP_r2_SD_min, nDP_r2_SD_plus))
# 
# # Region 3
# normalDP_r3 = daysPrecipMountCount_region3[[5:34]]
# normalDP_r3_med = median(values(normalDP_r3), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(daysPrecipMountCount_region3[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalDP_r3_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nDP_r3_SD_plus = normalDP_r3_med + (2*sd(med[5:34]))
# nDP_r3_SD_min = normalDP_r3_med - (2*sd(med[5:34]))
# 
# daysPrecip_diff[, 4] = diff
# daysPrecip_med = rbind(daysPrecip_med, c("østlandet", "daysPrecip", normalDP_r3_med, nDP_r3_SD_min, nDP_r3_SD_plus))
# 
# # Region 4
# normalDP_r4 = daysPrecipMountCount_region4[[5:34]]
# normalDP_r4_med = median(values(normalDP_r4), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(daysPrecipMountCount_region4[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalDP_r4_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nDP_r4_SD_plus = normalDP_r4_med + (2*sd(med[5:34]))
# nDP_r4_SD_min = normalDP_r4_med - (2*sd(med[5:34]))
# 
# daysPrecip_diff[, 5] = diff
# daysPrecip_med = rbind(daysPrecip_med, c("vestlandet", "daysPrecip", normalDP_r4_med, nDP_r4_SD_min, nDP_r4_SD_plus))
# 
# # Region 5
# normalDP_r5 = daysPrecipMountCount_region5[[5:34]]
# normalDP_r5_med = median(values(normalDP_r5), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(daysPrecipMountCount_region5[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalDP_r5_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nDP_r5_SD_plus = normalDP_r5_med + (2*sd(med[5:34]))
# nDP_r5_SD_min = normalDP_r5_med - (2*sd(med[5:34]))
# 
# daysPrecip_diff[, 6] = diff
# daysPrecip_med = rbind(daysPrecip_med, c("sørlandet", "daysPrecip", normalDP_r5_med, nDP_r5_SD_min, nDP_r5_SD_plus))
# 
# daysPrecip_diff$variable = "daysPrecip"
#write_xlsx(daysPrecip_diff, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/daysPrecip_diff.xlsx")
#write_xlsx(daysPrecip_med, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Dager med nedbør/daysPrecip_med.xlsx")

6.3 Beregne medianen for hvert år for gjennomsnittlig sommertemperatur

# Region 1
meanSummerMountCount_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/", full.names = TRUE, pattern = "region1_1957_2020.grd")[1]
meanSummerMountCount_region1 = stack(meanSummerMountCount_region1)
years = c(1:64)
med = c()
for(i in years){
  temp_med = median(values(meanSummerMountCount_region1[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
meanSummer_val = data.frame(matrix(c(rep("Nord-Norge", 64), rep("gjennomsnittlig sommertemperatur", 64), rep("grader celsius", 64), 1957:2020, med), nrow = 64, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(meanSummer_val) = columnNames

# Region 2
meanSummerMountCount_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/", full.names = TRUE, pattern = "region2_1957_2020.grd")[1]
meanSummerMountCount_region2 = stack(meanSummerMountCount_region2)
years = c(1:64)
med = c()
for(i in years){
  temp_med = median(values(meanSummerMountCount_region2[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
meanSummer_val_region2 = data.frame(matrix(c(rep("Midt-Norge", 64), rep("gjennomsnittlig sommertemperatur", 64), rep("grader celsius", 64), 1957:2020, med), nrow = 64, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(meanSummer_val_region2) = columnNames
meanSummer_val = rbind(meanSummer_val, meanSummer_val_region2)

# Region 3
meanSummerMountCount_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/", full.names = TRUE, pattern = "region3_1957_2020.grd")[1]
meanSummerMountCount_region3 = stack(meanSummerMountCount_region3)
years = c(1:64)
med = c()
for(i in years){
  temp_med = median(values(meanSummerMountCount_region3[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
meanSummer_val_region3 = data.frame(matrix(c(rep("Østlandet", 64), rep("gjennomsnittlig sommertemperatur", 64), rep("grader celsius", 64), 1957:2020, med), nrow = 64, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(meanSummer_val_region3) = columnNames
meanSummer_val = rbind(meanSummer_val, meanSummer_val_region3)

# Region 4
meanSummerMountCount_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/", full.names = TRUE, pattern = "region4_1957_2020.grd")[1]
meanSummerMountCount_region4 = stack(meanSummerMountCount_region4)
years = c(1:64)
med = c()
for(i in years){
  temp_med = median(values(meanSummerMountCount_region4[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
meanSummer_val_region4 = data.frame(matrix(c(rep("Vestlandet", 64), rep("gjennomsnittlig sommertemperatur", 64), rep("grader celsius", 64), 1957:2020, med), nrow = 64, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(meanSummer_val_region4) = columnNames
meanSummer_val = rbind(meanSummer_val, meanSummer_val_region4)

# Region 5
meanSummerMountCount_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/", full.names = TRUE, pattern = "region5_1957_2020.grd")[1]
meanSummerMountCount_region5 = stack(meanSummerMountCount_region5)
years = c(1:64)
med = c()
for(i in years){
  temp_med = median(values(meanSummerMountCount_region5[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
meanSummer_val_region5 = data.frame(matrix(c(rep("Sørlandet", 64), rep("gjennomsnittlig sommertemperatur", 64), rep("grader celsius", 64), 1957:2020, med), nrow = 64, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(meanSummer_val_region5) = columnNames
meanSummer_val = rbind(meanSummer_val, meanSummer_val_region5)

write_xlsx(meanSummer_val, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/medianMeanSummer.xlsx")
# 
# # Region 1
# normalMS_r1 = meanSummerMountCount_region1[[5:34]]
# normalMS_r1_med = median(values(normalMS_r1), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(meanSummerMountCount_region1[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalMS_r1_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nMS_r1_SD_plus = normalMS_r1_med + (2*sd(med[5:34]))
# nMS_r1_SD_min = normalMS_r1_med - (2*sd(med[5:34]))
# meanSummer_med = data.frame(matrix(c("nord-norge", "meanSummer", normalMS_r1_med, nMS_r1_SD_min, nMS_r1_SD_plus), nrow = 1))
# colnames(meanSummer_med) = c("area", "variable", "norm_med", "norm_-2SD", "norm_+2SD")
# 
# meanSummer_diff = data.frame(c(1957:2020),matrix(rep(NA, 64*5), nrow = 64, ncol = 5))
# columnNames = c("year","nord-norge", "midt-norge", "østlandet", "vestlandet", "sørlandet")
# colnames(meanSummer_diff) = columnNames
# meanSummer_diff[, 2] = diff
# 
# # Region 2
# normalMS_r2 = meanSummerMountCount_region2[[5:34]]
# normalMS_r2_med = median(values(normalMS_r2), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(meanSummerMountCount_region2[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalMS_r2_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nMS_r2_SD_plus = normalMS_r2_med + (2*sd(med[5:34]))
# nMS_r2_SD_min = normalMS_r2_med - (2*sd(med[5:34]))
# 
# meanSummer_diff[, 3] = diff
# meanSummer_med = rbind(meanSummer_med, c("midt-norge", "meanSummer", normalMS_r2_med, nMS_r2_SD_min, nMS_r2_SD_plus))
# 
# # Region 3
# normalMS_r3 = meanSummerMountCount_region3[[5:34]]
# normalMS_r3_med = median(values(normalMS_r3), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(meanSummerMountCount_region3[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalMS_r3_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nMS_r3_SD_plus = normalMS_r3_med + (2*sd(med[5:34]))
# nMS_r3_SD_min = normalMS_r3_med - (2*sd(med[5:34]))
# 
# meanSummer_diff[, 4] = diff
# meanSummer_med = rbind(meanSummer_med, c("østlandet", "meanSummer", normalMS_r3_med, nMS_r3_SD_min, nMS_r3_SD_plus))
# 
# # Region 4
# normalMS_r4 = meanSummerMountCount_region4[[5:34]]
# normalMS_r4_med = median(values(normalMS_r4), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(meanSummerMountCount_region4[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalMS_r4_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nMS_r4_SD_plus = normalMS_r4_med + (2*sd(med[5:34]))
# nMS_r4_SD_min = normalMS_r4_med - (2*sd(med[5:34]))
# 
# meanSummer_diff[, 5] = diff
# meanSummer_med = rbind(meanSummer_med, c("vestlandet", "meanSummer", normalMS_r4_med, nMS_r4_SD_min, nMS_r4_SD_plus))
# 
# # Region 5
# normalMS_r5 = meanSummerMountCount_region5[[5:34]]
# normalMS_r5_med = median(values(normalMS_r5), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(meanSummerMountCount_region5[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalMS_r5_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nMS_r5_SD_plus = normalMS_r5_med + (2*sd(med[5:34]))
# nMS_r5_SD_min = normalMS_r5_med - (2*sd(med[5:34]))
# 
# meanSummer_diff[, 6] = diff
# meanSummer_med = rbind(meanSummer_med, c("sørlandet", "meanSummer", normalMS_r5_med, nMS_r5_SD_min, nMS_r5_SD_plus))
# 
# meanSummer_diff$variable = "meanSummer"
#write_xlsx(meanSummer_diff, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/meanSummer_diff.xlsx")
#write_xlsx(meanSummer_med, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt sommer/meanSummer_med.xlsx")

6.4 Beregne medianen for hvert år for gjennomsnittlig vintertemperatur

# Region 1
meanWinterMountCount_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/", full.names = TRUE, pattern = "region1_1957_2020.grd")[1]
meanWinterMountCount_region1 = stack(meanWinterMountCount_region1)
years = c(1:63)
med = c()
winterSeason = c(1957:2020)
ws = c()
for(i in years){
  temp_med = median(values(meanWinterMountCount_region1[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
  ws = append(ws, paste(winterSeason[i], winterSeason[i+1], sep = "-"))
}
meanWinter_val = data.frame(matrix(c(rep("Nord-Norge", 63), rep("gjennomsnittlig vintertemperatur", 63), rep("grader celsius", 63), ws, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(meanWinter_val) = columnNames

# Region 2
meanWinterMountCount_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/", full.names = TRUE, pattern = "region2_1957_2020.grd")[1]
meanWinterMountCount_region2 = stack(meanWinterMountCount_region2)
years = c(1:63)
med = c()
winterSeason = c(1957:2020)
ws = c()
for(i in years){
  temp_med = median(values(meanWinterMountCount_region2[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
  ws = append(ws, paste(winterSeason[i], winterSeason[i+1], sep = "-"))
}
meanWinter_val_region2 = data.frame(matrix(c(rep("Midt-Norge", 63), rep("gjennomsnittlig vintertemperatur", 63), rep("grader celsius", 63), ws, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(meanWinter_val_region2) = columnNames
meanWinter_val = rbind(meanWinter_val, meanWinter_val_region2)

# Region 3
meanWinterMountCount_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/", full.names = TRUE, pattern = "region3_1957_2020.grd")[1]
meanWinterMountCount_region3 = stack(meanWinterMountCount_region3)
years = c(1:63)
med = c()
winterSeason = c(1957:2020)
ws = c()
for(i in years){
  temp_med = median(values(meanWinterMountCount_region3[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
  ws = append(ws, paste(winterSeason[i], winterSeason[i+1], sep = "-"))
}
meanWinter_val_region3 = data.frame(matrix(c(rep("Østlandet", 63), rep("gjennomsnittlig vintertemperatur", 63), rep("grader celsius", 63), ws, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(meanWinter_val_region3) = columnNames
meanWinter_val = rbind(meanWinter_val, meanWinter_val_region3)

# Region 4
meanWinterMountCount_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/", full.names = TRUE, pattern = "region4_1957_2020.grd")[1]
meanWinterMountCount_region4 = stack(meanWinterMountCount_region4)
years = c(1:63)
med = c()
winterSeason = c(1957:2020)
ws = c()
for(i in years){
  temp_med = median(values(meanWinterMountCount_region4[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
  ws = append(ws, paste(winterSeason[i], winterSeason[i+1], sep = "-"))
}
meanWinter_val_region4 = data.frame(matrix(c(rep("Vestlandet", 63), rep("gjennomsnittlig vintertemperatur", 63), rep("grader celsius", 63), ws, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(meanWinter_val_region4) = columnNames
meanWinter_val = rbind(meanWinter_val, meanWinter_val_region4)

# Region 5
meanWinterMountCount_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/", full.names = TRUE, pattern = "region5_1957_2020.grd")[1]
meanWinterMountCount_region5 = stack(meanWinterMountCount_region5)
years = c(1:63)
med = c()
winterSeason = c(1957:2020)
ws = c()
for(i in years){
  temp_med = median(values(meanWinterMountCount_region5[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
  ws = append(ws, paste(winterSeason[i], winterSeason[i+1], sep = "-"))
}
meanWinter_val_region5 = data.frame(matrix(c(rep("Sørlandet", 63), rep("gjennomsnittlig vintertemperatur", 63), rep("grader celsius", 63), ws, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(meanWinter_val_region5) = columnNames
meanWinter_val = rbind(meanWinter_val, meanWinter_val_region5)

write_xlsx(meanWinter_val, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/medianmeanWinter.xlsx")

# # Region 1
# normalMW_r1 = meanWinterMountCount_region1[[5:34]]
# normalMW_r1_med = median(values(normalMW_r1), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(meanWinterMountCount_region1[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalMW_r1_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nMW_r1_SD_plus = normalMW_r1_med + (2*sd(med[5:34]))
# nMW_r1_SD_min = normalMW_r1_med - (2*sd(med[5:34]))
# meanWinter_med = data.frame(matrix(c("nord-norge", "meanWinter", normalMW_r1_med, nMW_r1_SD_min, nMW_r1_SD_plus), nrow = 1))
# colnames(meanWinter_med) = c("area", "variable", "norm_med", "norm_-2SD", "norm_+2SD")
# 
# meanWinter_diff = data.frame(c(1957:2020),matrix(rep(NA, 64*5), nrow = 64, ncol = 5))
# columnNames = c("year","nord-norge", "midt-norge", "østlandet", "vestlandet", "sørlandet")
# colnames(meanWinter_diff) = columnNames
# meanWinter_diff[, 2] = diff
# 
# # Region 2
# normalMW_r2 = meanWinterMountCount_region2[[5:34]]
# normalMW_r2_med = median(values(normalMW_r2), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(meanWinterMountCount_region2[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalMW_r2_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nMW_r2_SD_plus = normalMW_r2_med + (2*sd(med[5:34]))
# nMW_r2_SD_min = normalMW_r2_med - (2*sd(med[5:34]))
# 
# meanWinter_diff[, 3] = diff
# meanWinter_med = rbind(meanWinter_med, c("midt-norge", "meanWinter", normalMW_r2_med, nMW_r2_SD_min, nMW_r2_SD_plus))
# 
# # Region 3
# normalMW_r3 = meanWinterMountCount_region3[[5:34]]
# normalMW_r3_med = median(values(normalMW_r3), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(meanWinterMountCount_region3[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalMW_r3_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nMW_r3_SD_plus = normalMW_r3_med + (2*sd(med[5:34]))
# nMW_r3_SD_min = normalMW_r3_med - (2*sd(med[5:34]))
# 
# meanWinter_diff[, 4] = diff
# meanWinter_med = rbind(meanWinter_med, c("østlandet", "meanWinter", normalMW_r3_med, nMW_r3_SD_min, nMW_r3_SD_plus))
# 
# # Region 4
# normalMW_r4 = meanWinterMountCount_region4[[5:34]]
# normalMW_r4_med = median(values(normalMW_r4), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(meanWinterMountCount_region4[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalMW_r4_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nMW_r4_SD_plus = normalMW_r4_med + (2*sd(med[5:34]))
# nMW_r4_SD_min = normalMW_r4_med - (2*sd(med[5:34]))
# 
# meanWinter_diff[, 5] = diff
# meanWinter_med = rbind(meanWinter_med, c("vestlandet", "meanWinter", normalMW_r4_med, nMW_r4_SD_min, nMW_r4_SD_plus))
# 
# # Region 5
# normalMW_r5 = meanWinterMountCount_region5[[5:34]]
# normalMW_r5_med = median(values(normalMW_r5), na.rm = TRUE)
# 
# years = c(1:64)
# med = c()
# diff = c()
# for(i in years){
#   
#   temp_med = median(values(meanWinterMountCount_region5[[i]]), na.rm = TRUE)
#   med = append(med, temp_med)
#   
#   temp_diff = temp_med - normalMW_r5_med
#   diff = append(diff, temp_diff)
#   
# }
# 
# nMW_r5_SD_plus = normalMW_r5_med + (2*sd(med[5:34]))
# nMW_r5_SD_min = normalMW_r5_med - (2*sd(med[5:34]))
# 
# meanWinter_diff[, 6] = diff
# meanWinter_med = rbind(meanWinter_med, c("sørlandet", "meanWinter", normalMW_r5_med, nMW_r5_SD_min, nMW_r5_SD_plus))
# 
# meanWinter_diff$variable = "meanWinter"
#write_xlsx(meanWinter_diff, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/meanWinter_diff.xlsx")
#write_xlsx(meanWinter_med, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Gjennomsnitt vinter/meanWinter_med.xlsx")

6.5 Beregne medianen for hvert år for dager med snødekke

# Region 1
snowCoverMountCount_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", full.names = TRUE, pattern = "region1_1957_2020.grd")[1]
snowCoverMountCount_region1 = stack(snowCoverMountCount_region1)
years = c(1:63)
med = c()
winterSeason = c(1957:2020)
ws = c()
for(i in years){
  temp_med = median(values(snowCoverMountCount_region1[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
  ws = append(ws, paste(winterSeason[i], winterSeason[i+1], sep = "-"))
}
snowCover_val = data.frame(matrix(c(rep("Nord-Norge", 63), rep("snødekke", 63), rep("dager", 63), ws, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(snowCover_val) = columnNames

# Region 2
snowCoverMountCount_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", full.names = TRUE, pattern = "region2_1957_2020.grd")[1]
snowCoverMountCount_region2 = stack(snowCoverMountCount_region2)
years = c(1:63)
med = c()
winterSeason = c(1957:2020)
ws = c()
for(i in years){
  temp_med = median(values(snowCoverMountCount_region2[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
  ws = append(ws, paste(winterSeason[i], winterSeason[i+1], sep = "-"))
}
snowCover_val_region2 = data.frame(matrix(c(rep("Midt-Norge", 63), rep("snødekke", 63), rep("dager", 63), ws, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(snowCover_val_region2) = columnNames
snowCover_val = rbind(snowCover_val, snowCover_val_region2)

# Region 3
snowCoverMountCount_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", full.names = TRUE, pattern = "region3_1957_2020.grd")[1]
snowCoverMountCount_region3 = stack(snowCoverMountCount_region3)
years = c(1:63)
med = c()
winterSeason = c(1957:2020)
ws = c()
for(i in years){
  temp_med = median(values(snowCoverMountCount_region3[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
  ws = append(ws, paste(winterSeason[i], winterSeason[i+1], sep = "-"))
}
snowCover_val_region3 = data.frame(matrix(c(rep("Østlandet", 63), rep("snødekke", 63), rep("dager", 63), ws, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(snowCover_val_region3) = columnNames
snowCover_val = rbind(snowCover_val, snowCover_val_region3)

# Region 4
snowCoverMountCount_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", full.names = TRUE, pattern = "region4_1957_2020.grd")[1]
snowCoverMountCount_region4 = stack(snowCoverMountCount_region4)
years = c(1:63)
med = c()
winterSeason = c(1957:2020)
ws = c()
for(i in years){
  temp_med = median(values(snowCoverMountCount_region4[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
  ws = append(ws, paste(winterSeason[i], winterSeason[i+1], sep = "-"))
}
snowCover_val_region4 = data.frame(matrix(c(rep("Vestlandet", 63), rep("snødekke", 63), rep("dager", 63), ws, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(snowCover_val_region4) = columnNames
snowCover_val = rbind(snowCover_val, snowCover_val_region4)

# Region 5
snowCoverMountCount_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", full.names = TRUE, pattern = "region5_1957_2020.grd")[1]
snowCoverMountCount_region5 = stack(snowCoverMountCount_region5)
years = c(1:63)
med = c()
winterSeason = c(1957:2020)
ws = c()
for(i in years){
  temp_med = median(values(snowCoverMountCount_region5[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
  ws = append(ws, paste(winterSeason[i], winterSeason[i+1], sep = "-"))
}
snowCover_val_region5 = data.frame(matrix(c(rep("Sørlandet", 63), rep("snødekke", 63), rep("dager", 63), ws, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(snowCover_val_region5) = columnNames
snowCover_val = rbind(snowCover_val, snowCover_val_region5)

write_xlsx(snowCover_val, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/mediansnowCover.xlsx")

# # Load data
# daysSnowCoverMountCount_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", full.names = TRUE, pattern = "region1_1957_2020.grd")[1]
# daysSnowCoverMountCount_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", full.names = TRUE, pattern = "region2_1957_2020.grd")[1]
# daysSnowCoverMountCount_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", full.names = TRUE, pattern = "region3_1957_2020.grd")[1]
# daysSnowCoverMountCount_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", full.names = TRUE, pattern = "region4_1957_2020.grd")[1]
# daysSnowCoverMountCount_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/", full.names = TRUE, pattern = "region5_1957_2020.grd")[1]
# 
# daysSnowCoverMountCount_region1 = stack(daysSnowCoverMountCount_region1)
# daysSnowCoverMountCount_region2 = stack(daysSnowCoverMountCount_region2)
# daysSnowCoverMountCount_region3 = stack(daysSnowCoverMountCount_region3)
# daysSnowCoverMountCount_region4 = stack(daysSnowCoverMountCount_region4)
# daysSnowCoverMountCount_region5 = stack(daysSnowCoverMountCount_region5)
# snowCover_vec = c(daysSnowCoverMountCount_region1, daysSnowCoverMountCount_region2, daysSnowCoverMountCount_region3, daysSnowCoverMountCount_region4, daysSnowCoverMountCount_region5)
# 
# diff_vec = list()
# for(i in 1:5){
#   norm_snowCover = snowCover_vec[[i]][[5:34]] # normal period
#   curr_snowCover = snowCover_vec[[i]][[59:63]] # current period
#   temp_vec = c()
#   
#   for(j in 1:10000){
#     norm_samp = sample(c(1:30), 1)
#     curr_samp = sample(c(1:5), 1)
#     
#     temp_norm = median(values(norm_snowCover[[norm_samp]]), na.rm = TRUE) # median for random year in the normal period 
#     temp_curr = median(values(curr_snowCover[[curr_samp]]), na.rm = TRUE) # median for random year in the current period
#     
#     temp_diff = temp_curr - temp_norm
#     temp_vec = append(temp_vec, temp_diff)
#   }
#   diff_vec[[i]] = temp_vec
# }
# snowCover_r1 = diff_vec[[1]]
# snowCover_r2 = diff_vec[[2]]
# snowCover_r3 = diff_vec[[3]]
# snowCover_r4 = diff_vec[[4]]
# snowCover_r5 = diff_vec[[5]]
# snowCover_df = cbind(snowCover_r1, snowCover_r2, snowCover_r3, snowCover_r4, snowCover_r5)
# snowCover_df = data.frame(snowCover_df)
# snowCover_df = as_tibble(snowCover_df) %>% rename("nord-norge" = snowCover_r1, "midt-norge" = snowCover_r2, "østlandet" = snowCover_r3, "vestlandet" = snowCover_r4, "sørlandet" = snowCover_r5) %>% mutate(variable = "snowCover")
# #write_xlsx(snowCover_df, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødekning/snowCover_diff.xlsx")

6.6 Beregne medianen for hvert år for vekstsesong

# Region 1
growthSeasonMountCount_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", full.names = TRUE, pattern = "region1_1957_2020.grd")[1]
growthSeasonMountCount_region1 = stack(growthSeasonMountCount_region1)
years = c(1:64)
med = c()
for(i in 2:length(years)){
  temp_med = median(values(growthSeasonMountCount_region1[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
growthSeason_val = data.frame(matrix(c(rep("Nord-Norge", 63), rep("vekstsesong", 63), rep("dager", 63), 1958:2020, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(growthSeason_val) = columnNames

# Region 2
growthSeasonMountCount_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", full.names = TRUE, pattern = "region2_1957_2020.grd")[1]
growthSeasonMountCount_region2 = stack(growthSeasonMountCount_region2)
years = c(1:64)
med = c()
for(i in 2:length(years)){
  temp_med = median(values(growthSeasonMountCount_region2[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
growthSeason_val_region2 = data.frame(matrix(c(rep("Midt-Norge", 63), rep("vekstsesong", 63), rep("dager", 63), 1958:2020, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(growthSeason_val_region2) = columnNames
growthSeason_val = rbind(growthSeason_val, growthSeason_val_region2)

# Region 3
growthSeasonMountCount_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", full.names = TRUE, pattern = "region3_1957_2020.grd")[1]
growthSeasonMountCount_region3 = stack(growthSeasonMountCount_region3)
years = c(1:64)
med = c()
for(i in 2:length(years)){
  temp_med = median(values(growthSeasonMountCount_region3[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
growthSeason_val_region3 = data.frame(matrix(c(rep("Østlandet", 63), rep("vekstsesong", 63), rep("dager", 63), 1958:2020, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(growthSeason_val_region3) = columnNames
growthSeason_val = rbind(growthSeason_val, growthSeason_val_region3)

# Region 4
growthSeasonMountCount_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", full.names = TRUE, pattern = "region4_1957_2020.grd")[1]
growthSeasonMountCount_region4 = stack(growthSeasonMountCount_region4)
years = c(1:64)
med = c()
for(i in 2:length(years)){
  temp_med = median(values(growthSeasonMountCount_region4[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
growthSeason_val_region4 = data.frame(matrix(c(rep("Vestlandet", 63), rep("vekstsesong", 63), rep("dager", 63), 1958:2020, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(growthSeason_val_region4) = columnNames
growthSeason_val = rbind(growthSeason_val, growthSeason_val_region4)

# Region 5
growthSeasonMountCount_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", full.names = TRUE, pattern = "region5_1957_2020.grd")[1]
growthSeasonMountCount_region5 = stack(growthSeasonMountCount_region5)
years = c(1:64)
med = c()
for(i in 2:length(years)){
  temp_med = median(values(growthSeasonMountCount_region5[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
growthSeason_val_region5 = data.frame(matrix(c(rep("Sørlandet", 63), rep("vekstsesong", 63), rep("dager", 63), 1958:2020, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(growthSeason_val_region5) = columnNames
growthSeason_val = rbind(growthSeason_val, growthSeason_val_region5)

write_xlsx(growthSeason_val, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/growthSeason_med.xlsx")

# # Load data
# growthSeasonMountCount_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", full.names = TRUE, pattern = "region1_1957_2020.grd")[1]
# growthSeasonMountCount_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", full.names = TRUE, pattern = "region2_1957_2020.grd")[1]
# growthSeasonMountCount_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", full.names = TRUE, pattern = "region3_1957_2020.grd")[1]
# growthSeasonMountCount_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", full.names = TRUE, pattern = "region4_1957_2020.grd")[1]
# growthSeasonMountCount_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/", full.names = TRUE, pattern = "region5_1957_2020.grd")[1]
# 
# # Stack the rasters
# growthSeasonMountCount_region1 = stack(growthSeasonMountCount_region1)
# growthSeasonMountCount_region2 = stack(growthSeasonMountCount_region2)
# growthSeasonMountCount_region3 = stack(growthSeasonMountCount_region3)
# growthSeasonMountCount_region4 = stack(growthSeasonMountCount_region4)
# growthSeasonMountCount_region5 = stack(growthSeasonMountCount_region5)
# growthSeason_vec = c(growthSeasonMountCount_region1, growthSeasonMountCount_region2, growthSeasonMountCount_region3, growthSeasonMountCount_region4, growthSeasonMountCount_region5)
# 
# # Create normal period and current period for each region
# diff_vec = list()
# for(i in 1:5){
#   norm_growthSeason = growthSeason_vec[[i]][[5:34]] # normal period
#   curr_growthSeason = growthSeason_vec[[i]][[60:64]] # current period
#   temp_vec = c()
#   
#   for(j in 1:10000){
#     norm_samp = sample(c(1:30), 1)
#     curr_samp = sample(c(1:5), 1)
#     
#     temp_norm = median(values(norm_growthSeason[[norm_samp]]), na.rm = TRUE) # median for random year in the normal period 
#     temp_curr = median(values(curr_growthSeason[[curr_samp]]), na.rm = TRUE) # median for random year in the current period
#     
#     temp_diff = temp_curr - temp_norm
#     temp_vec = append(temp_vec, temp_diff)
#   }
#   diff_vec[[i]] = temp_vec
# }
# growthSeason_r1 = diff_vec[[1]]
# growthSeason_r2 = diff_vec[[2]]
# growthSeason_r3 = diff_vec[[3]]
# growthSeason_r4 = diff_vec[[4]]
# growthSeason_r5 = diff_vec[[5]]
# growthSeason_df = cbind(growthSeason_r1, growthSeason_r2, growthSeason_r3, growthSeason_r4, growthSeason_r5)
# growthSeason_df = data.frame(growthSeason_df)
# growthSeason_df = as_tibble(growthSeason_df) %>% rename("nord-norge" = growthSeason_r1, "midt-norge" = growthSeason_r2, "østlandet" = growthSeason_r3, "vestlandet" = growthSeason_r4, "sørlandet" = growthSeason_r5) %>% mutate(variable = "growthSeason")
# #write_xlsx(growthSeason_df, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vekstsesong/growthSeason_diff.xlsx")

6.7 Beregne medianen for hvert år for vinterregn

# Region 1
winterRainMountCount_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", full.names = TRUE, pattern = "region1_1957_2020.grd")[1]
winterRainMountCount_region1 = stack(winterRainMountCount_region1)
years = c(1:64)
med = c()
for(i in 2:length(years)){
  temp_med = median(values(winterRainMountCount_region1[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
winterRain_val = data.frame(matrix(c(rep("Nord-Norge", 63), rep("vinterregn", 63), rep("mm", 63), 1958:2020, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(winterRain_val) = columnNames

# Region 2
winterRainMountCount_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", full.names = TRUE, pattern = "region2_1957_2020.grd")[1]
winterRainMountCount_region2 = stack(winterRainMountCount_region2)
years = c(1:64)
med = c()
for(i in 2:length(years)){
  temp_med = median(values(winterRainMountCount_region2[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
winterRain_val_region2 = data.frame(matrix(c(rep("Midt-Norge", 63), rep("Vinterregn", 63), rep("mm", 63), 1958:2020, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(winterRain_val_region2) = columnNames
winterRain_val = rbind(winterRain_val, winterRain_val_region2)

# Region 3
winterRainMountCount_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", full.names = TRUE, pattern = "region3_1957_2020.grd")[1]
winterRainMountCount_region3 = stack(winterRainMountCount_region3)
years = c(1:64)
med = c()
for(i in 2:length(years)){
  temp_med = median(values(winterRainMountCount_region3[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
winterRain_val_region3 = data.frame(matrix(c(rep("Østlandet", 63), rep("Vinterregn", 63), rep("mm", 63), 1958:2020, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(winterRain_val_region3) = columnNames
winterRain_val = rbind(winterRain_val, winterRain_val_region3)

# Region 4
winterRainMountCount_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", full.names = TRUE, pattern = "region4_1957_2020.grd")[1]
winterRainMountCount_region4 = stack(winterRainMountCount_region4)
years = c(1:64)
med = c()
for(i in 2:length(years)){
  temp_med = median(values(winterRainMountCount_region4[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
winterRain_val_region4 = data.frame(matrix(c(rep("Vestlandet", 63), rep("Vinterregn", 63), rep("mm", 63), 1958:2020, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(winterRain_val_region4) = columnNames
winterRain_val = rbind(winterRain_val, winterRain_val_region4)

# Region 5
winterRainMountCount_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", full.names = TRUE, pattern = "region5_1957_2020.grd")[1]
winterRainMountCount_region5 = stack(winterRainMountCount_region5)
years = c(1:64)
med = c()
for(i in 2:length(years)){
  temp_med = median(values(winterRainMountCount_region5[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
winterRain_val_region5 = data.frame(matrix(c(rep("Sørlandet", 63), rep("Vinterregn", 63), rep("mm", 63), 1958:2020, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(winterRain_val_region5) = columnNames
winterRain_val = rbind(winterRain_val, winterRain_val_region5)

winterRain_val$value = as.numeric(winterRain_val$value) / 10 
write_xlsx(winterRain_val, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/winterRain_med.xlsx")

# # Load data
# winterRainMountCount_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", full.names = TRUE, pattern = "region1_1957_2020.grd")[1]
# winterRainMountCount_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", full.names = TRUE, pattern = "region2_1957_2020.grd")[1]
# winterRainMountCount_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", full.names = TRUE, pattern = "region3_1957_2020.grd")[1]
# winterRainMountCount_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", full.names = TRUE, pattern = "region4_1957_2020.grd")[1]
# winterRainMountCount_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/", full.names = TRUE, pattern = "region5_1957_2020.grd")[1]
# 
# # Stack the rasters
# winterRainMountCount_region1 = stack(winterRainMountCount_region1)
# winterRainMountCount_region2 = stack(winterRainMountCount_region2)
# winterRainMountCount_region3 = stack(winterRainMountCount_region3)
# winterRainMountCount_region4 = stack(winterRainMountCount_region4)
# winterRainMountCount_region5 = stack(winterRainMountCount_region5)
# 
# # Set condition of > 0mm precipitation and > 2 degrees centigrade
# for(i in 1:64){
#   vals = values(winterRainMountCount_region1[[i]])
#   vals_rows = which(vals == 0)
#   vals[vals_rows] = NA
#   values(winterRainMountCount_region1[[i]]) = vals
# }
# 
# for(i in 1:64){
#   vals = values(winterRainMountCount_region2[[i]])
#   vals_rows = which(vals == 0)
#   vals[vals_rows] = NA
#   values(winterRainMountCount_region2[[i]]) = vals
# }
# 
# for(i in 1:64){
#   vals = values(winterRainMountCount_region3[[i]])
#   vals_rows = which(vals == 0)
#   vals[vals_rows] = NA
#   values(winterRainMountCount_region3[[i]]) = vals
# }
# 
# for(i in 1:64){
#   vals = values(winterRainMountCount_region4[[i]])
#   vals_rows = which(vals == 0)
#   vals[vals_rows] = NA
#   values(winterRainMountCount_region4[[i]]) = vals
# }
# 
# for(i in 1:64){
#   vals = values(winterRainMountCount_region5[[i]])
#   vals_rows = which(vals == 0)
#   vals[vals_rows] = NA
#   values(winterRainMountCount_region5[[i]]) = vals
# }
# 
# 
# winterRain_vec = c(winterRainMountCount_region1, winterRainMountCount_region2, winterRainMountCount_region3, winterRainMountCount_region4, winterRainMountCount_region5)
# 
# # Create normal period and current period for each region
# diff_vec = list()
# for(i in 1:5){
#   norm_winterRain = winterRain_vec[[i]][[5:34]] # normal period
#   curr_winterRain = winterRain_vec[[i]][[60:64]] # current period
#   temp_vec = c()
#   
#   for(j in 1:10000){
#     norm_samp = sample(c(1:30), 1)
#     curr_samp = sample(c(1:5), 1)
#     
#     temp_norm = median(values(norm_winterRain[[norm_samp]]), na.rm = TRUE) # median for random year in the normal period 
#     temp_curr = median(values(curr_winterRain[[curr_samp]]), na.rm = TRUE) # median for random year in the current period
#     
#     temp_diff = temp_curr - temp_norm
#     temp_vec = append(temp_vec, temp_diff)
#   }
#   diff_vec[[i]] = temp_vec
# }
# winterRain_r1 = diff_vec[[1]]
# winterRain_r2 = diff_vec[[2]]
# winterRain_r3 = diff_vec[[3]]
# winterRain_r4 = diff_vec[[4]]
# winterRain_r5 = diff_vec[[5]]
# winterRain_df = cbind(winterRain_r1, winterRain_r2, winterRain_r3, winterRain_r4, winterRain_r5)
# winterRain_df = data.frame(winterRain_df)
# winterRain_df = as_tibble(winterRain_df) %>% rename("nord-norge" = winterRain_r1, "midt-norge" = winterRain_r2, "østlandet" = winterRain_r3, "vestlandet" = winterRain_r4, "sørlandet" = winterRain_r5) %>% mutate(variable = "winterRain")
#write_xlsx(winterRain_df, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Vinterregn/winterRain_diff.xlsx")

6.8 Beregne medianen for hvert år for snødybde

# Region 1
snowDepthMountCount_region1 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/", full.names = TRUE, pattern = "region1_1957_2020.grd")[1]
snowDepthMountCount_region1 = stack(snowDepthMountCount_region1)
years = c(1:63)
winterSeason = c(1957:2020)
ws = c()
med = c()
for(i in years){
  temp_med = median(values(snowDepthMountCount_region1[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
  ws = append(ws, paste(winterSeason[i], winterSeason[i+1], sep = "-"))
}
snowDepth_val = data.frame(matrix(c(rep("Nord-Norge", 63), rep("Snødybde", 63), rep("mm", 63), ws, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(snowDepth_val) = columnNames

# Region 2
snowDepthMountCount_region2 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/", full.names = TRUE, pattern = "region2_1957_2020.grd")[1]
snowDepthMountCount_region2 = stack(snowDepthMountCount_region2)
med = c()
for(i in years){
  temp_med = median(values(snowDepthMountCount_region2[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
snowDepth_val_region2 = data.frame(matrix(c(rep("Midt-Norge", 63), rep("Snødybde", 63), rep("mm", 63), ws, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(snowDepth_val_region2) = columnNames
snowDepth_val = rbind(snowDepth_val, snowDepth_val_region2)

# Region 3
snowDepthMountCount_region3 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/", full.names = TRUE, pattern = "region3_1957_2020.grd")[1]
snowDepthMountCount_region3 = stack(snowDepthMountCount_region3)
med = c()
for(i in years){
  temp_med = median(values(snowDepthMountCount_region3[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
snowDepth_val_region3 = data.frame(matrix(c(rep("Østlandet", 63), rep("Snødybde", 63), rep("mm", 63), ws, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(snowDepth_val_region3) = columnNames
snowDepth_val = rbind(snowDepth_val, snowDepth_val_region3)

# Region 4
snowDepthMountCount_region4 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/", full.names = TRUE, pattern = "region4_1957_2020.grd")[1]
snowDepthMountCount_region4 = stack(snowDepthMountCount_region4)
med = c()
for(i in years){
  temp_med = median(values(snowDepthMountCount_region4[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
snowDepth_val_region4 = data.frame(matrix(c(rep("Vestlandet", 63), rep("Snødybde", 63), rep("mm", 63), ws, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(snowDepth_val_region4) = columnNames
snowDepth_val = rbind(snowDepth_val, snowDepth_val_region4)

# Region 5
snowDepthMountCount_region5 = list.files("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/", full.names = TRUE, pattern = "region5_1957_2020.grd")[1]
snowDepthMountCount_region5 = stack(snowDepthMountCount_region5)
med = c()
for(i in years){
  temp_med = median(values(snowDepthMountCount_region5[[i]]), na.rm = TRUE)
  med = append(med, temp_med)
}
snowDepth_val_region5 = data.frame(matrix(c(rep("Sørlandet", 63), rep("Snødybde", 63), rep("mm", 63), ws, med), nrow = 63, ncol = 5))
columnNames = c("reg", "variable", "unit", "year", "value")
colnames(snowDepth_val_region5) = columnNames
snowDepth_val = rbind(snowDepth_val, snowDepth_val_region5)

write_xlsx(snowDepth_val, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/Klima/Snødybde/snowDepth_med.xlsx")