Skip to content

Commit

Permalink
update Day 3
Browse files Browse the repository at this point in the history
  • Loading branch information
@mkkallio
mkkallio committed Nov 3, 2021
1 parent b5c8d66 commit 1faa861
Show file tree
Hide file tree
Showing 4 changed files with 234 additions and 43 deletions.
29 changes: 25 additions & 4 deletions Day 3/CSC_D3S1 - Raster data in R.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,16 +33,18 @@ library(paletteer)
side <- 15
empty_mat <- matrix(0, nrow = side, ncol = side)
empty_mat
mat1 <- empty_mat
for(i in 1:side) {
mat1[,i] <- i
}
# look at the matrix, convert it to a raster and plot the raster
mat1
r1 <- terra::rast(mat1)
r1 <- terra::rast(t(mat1))
r1
terra::ext(r1) <- c(0,side,0,side) # set extent
r1
plot(r1)
```

Expand All @@ -53,7 +55,7 @@ Let's do another kind of raster, this time filled with random numbers.
mat2 <- matrix(runif(length(empty_mat),0,10), nrow=side, ncol=side)
r2 <- terra::rast(mat2)
terra::ext(r2) <- c(0,side,0,side) # set extent
plot(r2, col = paletteer::paletteer_c("grDevices::Viridis", 10))
plot(r2, col = paletteer::paletteer_c("grDevices::Viridis", 100))
```


Expand Down Expand Up @@ -81,6 +83,7 @@ r1^2
1 / r1
stack <- c(r1, r1+5, r1^2, 1/r1)
stack
names(stack) <- c("add 5", "square root", "square", "division")
plot(stack)
Expand All @@ -98,11 +101,12 @@ Let's read in the elevation model of Finland included in the data and inspect it
elevation <- terra::rast("../Data/raster/Elevation_10km.tif")
summary(elevation)
ext(elevation) # this is the boundinx box of the raster
ext(elevation) # this is the bounding box of the raster
plot(elevation, col = paletteer::paletteer_c("viridis::viridis", 100))
elevation[elevation == 0] <- NA
plot(elevation, col = paletteer::paletteer_c("viridis::viridis", 100))
summary(elevation)
```

We can compute statistics from the elevation layer using *global()* function.
Expand All @@ -114,6 +118,7 @@ global(elevation, median, na.rm=TRUE)
global(elevation, range, na.rm=TRUE)
global(elevation, sd, na.rm=TRUE)
plot(stack)
global(stack, mean)
```
Expand All @@ -131,6 +136,15 @@ remember: for CSC notebooks replace ".." with "/home/rstudio/r-spatial-course"
* plot a histogram using *hist()*

```{r}
t01 <- rast("../Data/raster/Tmean_Monthly_2005-2014/Tmon_20050101.tif")
t07 <- rast("../Data/raster/Tmean_Monthly_2005-2014/Tmon_20050701.tif")
t10 <- rast("../Data/raster/Tmean_Monthly_2005-2014/Tmon_20051001.tif")
t <- c(t01, t07, t10)
global(t, mean, na.rm=TRUE)
```

Expand All @@ -142,6 +156,13 @@ remember: for CSC notebooks replace ".." with "/home/rstudio/r-spatial-course"


```{r}
plot(t)
brk <- seq(-11,20, length.out = 10)
brk
plot(t, col = paletteer_c("scico::turku", 100), breaks = brk)
```

110 changes: 99 additions & 11 deletions Day 3/CSC_D3S2 - Raster_data_manipulation.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,13 @@ library(tmap)
# for CSC notebooks replace ".." with "/home/rstudio/r-spatial-course"
admin <- sf::read_sf("../Data/raster/SuomenKuntajako_2021_100k.shp")
plot(admin)
helsinki <- dplyr::filter(admin, NAMEFIN == "Helsinki")
plot(admin)
plot(helsinki)
```


Expand All @@ -31,24 +35,58 @@ Read in a list of files.
# for CSC notebooks replace ".." with "/home/rstudio/r-spatial-course"
files <- list.files("../Data/raster/Tmean_Monthly_2005-2014/",
full.names = TRUE)
files
temperature <- lapply(files, terra::rast) #read in all files, output is a list
length(temperature)
temperature[[1]]
temperature[1:5]
temperature <- do.call("c", temperature) # stack all rasters together
temperature
plot(temperature)
temperature <- terra::rast(files)
plot(temperature)
```


Before we do anything, let's make sure both datasets are in the same coordinate reference system. We'll use the CRS of temperature raster stack, because tranforming a raster changes its values.

```{r}
sf::st_crs(helsinki)
( sfcrs <- sf::st_crs(helsinki) )
terra::crs(temperature)
class(sfcrs)
str(sfcrs)
sfcrs$wkt
helsinki <- sf::st_transform(helsinki, terra::crs(temperature))
st_crs(helsinki)
# we could project the raster with:
# temperature <- terra::project(temperature, sf::st_crs(helsinki))
temperature <- terra::project(temperature, sf::st_crs(helsinki)$wkt)
temp <- temperature[[1]]
temp
plot(temp)
temp_projected <- project(temp, "epsg:4326" )
temp_reprojected <- project(temp_projected, temp )
plot( c(temp, temp_reprojected) )
c(temp, temp_reprojected)
all.equal(temp, temp_reprojected)
plot(temp_projected)
plot(temp_reprojected)
```


Expand All @@ -59,6 +97,8 @@ Now we can find out all the values in the raster stack for Helsinki! to do this
```{r}
helsinki_terra <- terra::vect(helsinki)
helsinki_temps <- terra::extract(temperature, helsinki_terra)
helsinki_temps
helsinki_temps <- dplyr::summarise_all(helsinki_temps, mean, na.rm=TRUE) %>%
tidyr::gather(date, temperature, -ID) %>%
Expand All @@ -67,6 +107,7 @@ helsinki_temps <- dplyr::summarise_all(helsinki_temps, mean, na.rm=TRUE) %>%
helsinki_temps
plot(helsinki_temps[,3])
lines(helsinki_temps[,3])
```


Expand All @@ -75,8 +116,10 @@ lines(helsinki_temps[,3])
Aggregation is done by providing the function a factor - how many cells in each direction are used to aggregate. The temperature has a resolution of 10km. Aggregating by factor 4 yields a raster with resolution of 40km (by 40km).

```{r}
temp_40km <- terra::aggregate(temperature, fact = 4, fun = mean)
temp_40km <- terra::aggregate(temperature, fact = 4, fun = mean, na.rm=TRUE)
plot(temp_40km)
temp_40km
```

We can also disaggregate - but the original values cannot be returned!
Expand All @@ -97,6 +140,8 @@ Resampling can achieve both. But are they the same - no!
temp_40km_resampled <- terra::resample(temperature, temp_40km, method = "bilinear")
plot(temp_40km_resampled)
temp_40km_resampled
temp_40km
plot(temp_40km - temp_40km_resampled)
```
Expand Down Expand Up @@ -129,9 +174,8 @@ However, it only runs for a single layer
# ntice the message saying it only used the first layer
temperature_smoothed <- terra::focal(temperature,
w = 3,
fun = mean,
na.rm=TRUE)
w = 21,
fun = max)
plot(temperature_smoothed)
plot(temperature[[1]])
Expand All @@ -141,9 +185,8 @@ plot(temperature[[1]])
We can also define the neighbourhood in other ways:
```{r}
temperature_smoothed <- terra::focal(temperature[[1]],
w = c(3,ncol(temperature)-1),
fun = mean,
na.rm=TRUE)
w = c(nrow(temperature)-1, 3),
fun = mean)
plot(temperature_smoothed)
# remove values from the sea
Expand All @@ -158,6 +201,7 @@ or with matrices:
```{r}
weight_matrix <- matrix(1, nrow = nrow(temperature)-1, ncol = 1)
weight_matrix
temperature_smoothed <- terra::focal(temperature[[1]],
w = weight_matrix,
Expand All @@ -180,22 +224,32 @@ The final part of this lesson shows reclassification. Here we divide Finland int
But before that, lets compute the mean temeprature at each cell over the 10-year period 2005-2014, and find "optimal" class intervals.

```{r}
plot(temperature)
temperature_mean <- mean(temperature)
temperature_mean
plot(temperature_mean)
library(classInt)
intervals <- classIntervals(terra::values(temperature_mean),
intervals <- classInt::classIntervals(terra::values(temperature_mean),
n = 4,
style = "quantile")
# quantile(values(temperature_mean), probs = c(0.2, 0.4, 0.6, 0.8), na.rm=TRUE)
intervals
reclass_matrix <- matrix(c(-Inf, intervals$brks[1], 1,
intervals$brks[1], intervals$brks[2], 2,
intervals$brks[2], intervals$brks[3], 3,
intervals$brks[3], intervals$brks[4], 4,
intervals$brks[4], Inf, 5),
ncol = 3,
byrow = TRUE)
reclass_matrix
temperature_classified <- terra::classify(temperature_mean,
reclass_matrix)
Expand All @@ -209,7 +263,8 @@ And finally, let's save that raster for later use.
```{r}
# for CSC notebooks replace ".." with "/home/rstudio/r-spatial-course"
terra::writeRaster(temperature_classified,
filename = "../Data/raster/temperature_classified.tif")
filename = "../Data/raster/temperature_classified.tif",
overwrite = TRUE)
```
Expand All @@ -225,6 +280,18 @@ terra::writeRaster(temperature_classified,
* find the mean, min, and max temperatures for that municipality over the period 2005-2014 (the full raster stack)

```{r}
kolari <- admin %>%
filter(NAMEFIN == "Kolari")
temp_kolari <- crop(temperature, kolari)
values <- global(temp_kolari, mean)
mean_temp <- mean(values$mean)
values <- terra::extract(temp_kolari, vect(kolari))
mean(as.matrix(values[,-1]),na.rm=TRUE)
```

Expand All @@ -237,6 +304,15 @@ terra::writeRaster(temperature_classified,


```{r}
focal_mean <- focal(temperature_mean,
w = matrix(c(0,1,0,1,-4,1,0,1,0), nrow=3))
plot(focal_mean)
focal_mean <- crop(focal_mean, kolari)
plot(c(focal_mean, mean(temp_kolari)))
plot(focal_mean - mean(temp_kolari))
```

Expand All @@ -248,6 +324,18 @@ terra::writeRaster(temperature_classified,

```{r}
reclass_matrix <- matrix( c(-Inf, mean(values(temp_kolari), na.rm=TRUE), factor("cold", levels = c("hot", "cold")),
mean(values(temp_kolari), na.rm=TRUE), Inf, factor("hot", levels = c("hot", "cold"))),
ncol = 3,
byrow = TRUE)
reclass_matrix
temp_kolari_classified <- classify(mean(temp_kolari),
reclass_matrix)
plot(mean(temp_kolari))
plot(temp_kolari_classified)
```


Loading

0 comments on commit 1faa861

Please sign in to comment.