Raster with R

## .font200[<svg viewBox="0 0 581 512" style="height:1em;fill:currentColor;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M581 226.6C581 119.1 450.9 32 290.5 32S0 119.1 0 226.6C0 322.4 103.3 402 239.4 418.1V480h99.1v-61.5c24.3-2.7 47.6-7.4 69.4-13.9L448 480h112l-67.4-113.7c54.5-35.4 88.4-84.9 88.4-139.7zm-466.8 14.5c0-73.5 98.9-133 220.8-133s211.9 40.7 211.9 133c0 50.1-26.5 85-70.3 106.4-2.4-1.6-4.7-2.9-6.4-3.7-10.2-5.2-27.8-10.5-27.8-10.5s86.6-6.4 86.6-92.7-90.6-87.9-90.6-87.9h-199V361c-74.1-21.5-125.2-67.1-125.2-119.9zm225.1 38.3v-55.6c57.8 0 87.8-6.8 87.8 27.3 0 36.5-38.2 28.3-87.8 28.3zm-.9 72.5H365c10.8 0 18.9 11.7 24 19.2-16.1 1.9-33 2.8-50.6 2.9v-22.1z"></path></svg>asters]

<br>

.instructors[
  .font140[Manipulate and visualize spatial data with <svg viewBox="0 0 581 512" style="height:1em;fill:currentColor;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M581 226.6C581 119.1 450.9 32 290.5 32S0 119.1 0 226.6C0 322.4 103.3 402 239.4 418.1V480h99.1v-61.5c24.3-2.7 47.6-7.4 69.4-13.9L448 480h112l-67.4-113.7c54.5-35.4 88.4-84.9 88.4-139.7zm-466.8 14.5c0-73.5 98.9-133 220.8-133s211.9 40.7 211.9 133c0 50.1-26.5 85-70.3 106.4-2.4-1.6-4.7-2.9-6.4-3.7-10.2-5.2-27.8-10.5-27.8-10.5s86.6-6.4 86.6-92.7-90.6-87.9-90.6-87.9h-199V361c-74.1-21.5-125.2-67.1-125.2-119.9zm225.1 38.3v-55.6c57.8 0 87.8-6.8 87.8 27.3 0 36.5-38.2 28.3-87.8 28.3zm-.9 72.5H365c10.8 0 18.9 11.7 24 19.2-16.1 1.9-33 2.8-50.6 2.9v-22.1z"></path></svg>]
  <br><br><br>
  .authors140[David Beauchesne & Kevin Cazelles]
  <br><br>
  March 24, 2021
  <br>
  .font200[
  [<svg viewBox="0 0 496 512" style="height:1em;fill:currentColor;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M165.9 397.4c0 2-2.3 3.6-5.2 3.6-3.3.3-5.6-1.3-5.6-3.6 0-2 2.3-3.6 5.2-3.6 3-.3 5.6 1.3 5.6 3.6zm-31.1-4.5c-.7 2 1.3 4.3 4.3 4.9 2.6 1 5.6 0 6.2-2s-1.3-4.3-4.3-5.2c-2.6-.7-5.5.3-6.2 2.3zm44.2-1.7c-2.9.7-4.9 2.6-4.6 4.9.3 2 2.9 3.3 5.9 2.6 2.9-.7 4.9-2.6 4.6-4.6-.3-1.9-3-3.2-5.9-2.9zM244.8 8C106.1 8 0 113.3 0 252c0 110.9 69.8 205.8 169.5 239.2 12.8 2.3 17.3-5.6 17.3-12.1 0-6.2-.3-40.4-.3-61.4 0 0-70 15-84.7-29.8 0 0-11.4-29.1-27.8-36.6 0 0-22.9-15.7 1.6-15.4 0 0 24.9 2 38.6 25.8 21.9 38.6 58.6 27.5 72.9 20.9 2.3-16 8.8-27.1 16-33.7-55.9-6.2-112.3-14.3-112.3-110.5 0-27.5 7.6-41.3 23.6-58.9-2.6-6.5-11.1-33.3 2.6-67.9 20.9-6.5 69 27 69 27 20-5.6 41.5-8.5 62.8-8.5s42.8 2.9 62.8 8.5c0 0 48.1-33.6 69-27 13.7 34.7 5.2 61.4 2.6 67.9 16 17.7 25.8 31.5 25.8 58.9 0 96.5-58.9 104.2-114.8 110.5 9.2 7.9 17 22.9 17 46.4 0 33.7-.3 75.4-.3 83.6 0 6.5 4.6 14.4 17.3 12.1C428.2 457.8 496 362.9 496 252 496 113.3 383.5 8 244.8 8zM97.2 352.9c-1.3 1-1 3.3.7 5.2 1.6 1.6 3.9 2.3 5.2 1 1.3-1 1-3.3-.7-5.2-1.6-1.6-3.9-2.3-5.2-1zm-10.8-8.1c-.7 1.3.3 2.9 2.3 3.9 1.6 1 3.6.7 4.3-.7.7-1.3-.3-2.9-2.3-3.9-2-.6-3.6-.3-4.3.7zm32.4 35.6c-1.6 1.3-1 4.3 1.3 6.2 2.3 2.3 5.2 2.6 6.5 1 1.3-1.3.7-4.3-1.3-6.2-2.2-2.3-5.2-2.6-6.5-1zm-11.4-14.7c-1.6 1-1.6 3.6 0 5.9 1.6 2.3 4.3 3.3 5.6 2.3 1.6-1.3 1.6-3.9 0-6.2-1.4-2.3-4-3.3-5.6-2z"></path></svg>](https://github.com/inSilecoInc/wksh_raster_CWS)
  ]
]

[![deploy workshop](https://github.com/inSilecoInc/wksh_raster_CWS/actions/workflows/deploy.yml/badge.svg)](https://github.com/inSilecoInc/wksh_raster_CWS/actions/workflows/deploy.yml)

.instructors[Content under [<svg viewBox="0 0 496 512" style="height:1em;fill:currentColor;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M245.83 214.87l-33.22 17.28c-9.43-19.58-25.24-19.93-27.46-19.93-22.13 0-33.22 14.61-33.22 43.84 0 23.57 9.21 43.84 33.22 43.84 14.47 0 24.65-7.09 30.57-21.26l30.55 15.5c-6.17 11.51-25.69 38.98-65.1 38.98-22.6 0-73.96-10.32-73.96-77.05 0-58.69 43-77.06 72.63-77.06 30.72-.01 52.7 11.95 65.99 35.86zm143.05 0l-32.78 17.28c-9.5-19.77-25.72-19.93-27.9-19.93-22.14 0-33.22 14.61-33.22 43.84 0 23.55 9.23 43.84 33.22 43.84 14.45 0 24.65-7.09 30.54-21.26l31 15.5c-2.1 3.75-21.39 38.98-65.09 38.98-22.69 0-73.96-9.87-73.96-77.05 0-58.67 42.97-77.06 72.63-77.06 30.71-.01 52.58 11.95 65.56 35.86zM247.56 8.05C104.74 8.05 0 123.11 0 256.05c0 138.49 113.6 248 247.56 248 129.93 0 248.44-100.87 248.44-248 0-137.87-106.62-248-248.44-248zm.87 450.81c-112.54 0-203.7-93.04-203.7-202.81 0-105.42 85.43-203.27 203.72-203.27 112.53 0 202.82 89.46 202.82 203.26-.01 121.69-99.68 202.82-202.84 202.82z"></path></svg> <svg viewBox="0 0 496 512" style="height:1em;fill:currentColor;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M314.9 194.4v101.4h-28.3v120.5h-77.1V295.9h-28.3V194.4c0-4.4 1.6-8.2 4.6-11.3 3.1-3.1 6.9-4.7 11.3-4.7H299c4.1 0 7.8 1.6 11.1 4.7 3.1 3.2 4.8 6.9 4.8 11.3zm-101.5-63.7c0-23.3 11.5-35 34.5-35s34.5 11.7 34.5 35c0 23-11.5 34.5-34.5 34.5s-34.5-11.5-34.5-34.5zM247.6 8C389.4 8 496 118.1 496 256c0 147.1-118.5 248-248.4 248C113.6 504 0 394.5 0 256 0 123.1 104.7 8 247.6 8zm.8 44.7C130.2 52.7 44.7 150.6 44.7 256c0 109.8 91.2 202.8 203.7 202.8 103.2 0 202.8-81.1 202.8-202.8.1-113.8-90.2-203.3-202.8-203.3z"></path></svg>](https://creativecommons.org/licenses/by/4.0/) unless otherwise specified.]

---

---

# <svg viewBox="0 0 512 512" style="height:1em;fill:currentColor;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M256 8C119.043 8 8 119.083 8 256c0 136.997 111.043 248 248 248s248-111.003 248-248C504 119.083 392.957 8 256 8zm0 110c23.196 0 42 18.804 42 42s-18.804 42-42 42-42-18.804-42-42 18.804-42 42-42zm56 254c0 6.627-5.373 12-12 12h-88c-6.627 0-12-5.373-12-12v-24c0-6.627 5.373-12 12-12h12v-64h-12c-6.627 0-12-5.373-12-12v-24c0-6.627 5.373-12 12-12h64c6.627 0 12 5.373 12 12v100h12c6.627 0 12 5.373 12 12v24z"></path></svg> Introduction to `stars`

![:custom_hr]()

## <svg viewBox="0 0 512 512" style="height:1em;fill:currentColor;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M256 8C119 8 8 119 8 256s111 248 248 248 248-111 248-248S393 8 256 8zm0 448c-110.5 0-200-89.5-200-200S145.5 56 256 56s200 89.5 200 200-89.5 200-200 200zm61.8-104.4l-84.9-61.7c-3.1-2.3-4.9-5.9-4.9-9.7V116c0-6.6 5.4-12 12-12h32c6.6 0 12 5.4 12 12v141.7l66.8 48.6c5.4 3.9 6.5 11.4 2.6 16.8L334.6 349c-3.9 5.3-11.4 6.5-16.8 2.6z"></path></svg> ~20min

---

# Package [`stars`](https://CRAN.R-project.org/package=stars): SpatioTemporal ARrayS

### Extent `sf` to Raster and Vector Datacubes

.center[[![:scale 62%](https://raw.githubusercontent.com/r-spatial/stars/master/images/cube3.png)](https://r-spatial.github.io/stars/index.html)] 
.right[.font80[see https://r-spatial.github.io/stars/index.html]]

---
# Example file

![Accessed 16/12/2020](img/gebco.png)

- https://download.gebco.net/ 
- <svg viewBox="0 0 512 512" style="height:1em;fill:currentColor;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M326.612 185.391c59.747 59.809 58.927 155.698.36 214.59-.11.12-.24.25-.36.37l-67.2 67.2c-59.27 59.27-155.699 59.262-214.96 0-59.27-59.26-59.27-155.7 0-214.96l37.106-37.106c9.84-9.84 26.786-3.3 27.294 10.606.648 17.722 3.826 35.527 9.69 52.721 1.986 5.822.567 12.262-3.783 16.612l-13.087 13.087c-28.026 28.026-28.905 73.66-1.155 101.96 28.024 28.579 74.086 28.749 102.325.51l67.2-67.19c28.191-28.191 28.073-73.757 0-101.83-3.701-3.694-7.429-6.564-10.341-8.569a16.037 16.037 0 0 1-6.947-12.606c-.396-10.567 3.348-21.456 11.698-29.806l21.054-21.055c5.521-5.521 14.182-6.199 20.584-1.731a152.482 152.482 0 0 1 20.522 17.197zM467.547 44.449c-59.261-59.262-155.69-59.27-214.96 0l-67.2 67.2c-.12.12-.25.25-.36.37-58.566 58.892-59.387 154.781.36 214.59a152.454 152.454 0 0 0 20.521 17.196c6.402 4.468 15.064 3.789 20.584-1.731l21.054-21.055c8.35-8.35 12.094-19.239 11.698-29.806a16.037 16.037 0 0 0-6.947-12.606c-2.912-2.005-6.64-4.875-10.341-8.569-28.073-28.073-28.191-73.639 0-101.83l67.2-67.19c28.239-28.239 74.3-28.069 102.325.51 27.75 28.3 26.872 73.934-1.155 101.96l-13.087 13.087c-4.35 4.35-5.769 10.79-3.783 16.612 5.864 17.194 9.042 34.999 9.69 52.721.509 13.906 17.454 20.446 27.294 10.606l37.106-37.106c59.271-59.259 59.271-155.699.001-214.959z"></path></svg> [data_and_script.zip](https://github.com/inSilecoInc/workshop_raster/raw/main/data_and_script.zip)

---
# Read with [`stars`](https://CRAN.R-project.org/package=stars)

```r
library(stars)
ras <- read_stars("data/bathy.tif")   # argument `driver` to specify the driver
class(ras)
#> [1] "stars"
ras
#> stars object with 2 dimensions and 1 attribute
#> attribute(s), summary of first 1e+05 cells:
#>  bathy.tif [m]   
#>  Min.   :-435.0  
#>  1st Qu.:  23.0  
#>  Median : 479.0  
#>  Mean   : 299.5  
#>  3rd Qu.: 540.0  
#>  Max.   : 812.0  
#> dimension(s):
#>   from   to offset       delta refsys point values x/y
#> x    1 3840    -71  0.00416667 WGS 84 FALSE   NULL [x]
#> y    1 2880     55 -0.00416667 WGS 84 FALSE   NULL [y]
```

---
# What do we have?

```r
dim(ras)
#>    x    y 
#> 3840 2880
st_bbox(ras)
#> xmin ymin xmax ymax 
#>  -71   43  -55   55
```

---
# What do we have?

```r
ras_crs <- st_crs(ras)
class(ras_crs)
#> [1] "crs"
ras_crs
#> Coordinate Reference System:
#>   User input: WGS 84 
#>   wkt:
#> GEOGCRS["WGS 84",
#>     DATUM["World Geodetic System 1984",
#>         ELLIPSOID["WGS 84",6378137,298.257223563,
#>             LENGTHUNIT["metre",1]]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433]],
#>     CS[ellipsoidal,2],
#>         AXIS["geodetic latitude (Lat)",north,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         AXIS["geodetic longitude (Lon)",east,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>     ID["EPSG",4326]]
```

---
# What do we have?

```r
ras_crs$input
#> [1] "WGS 84"
ras_crs$proj4string
#> [1] "+proj=longlat +datum=WGS84 +no_defs"
ras_crs$epsg
#> [1] 4326
```
--

- https://epsg.io/
- https://epsg.io/4326

---
# Conversion between [`stars`](https://CRAN.R-project.org/package=stars) and [`raster`](https://CRAN.R-project.org/package=raster)

```r
library(raster)
rar <- raster("data/bathy.tif")
ras_c <- st_as_stars(rar) 
class(ras_c)
#> [1] "stars"
```

```r
rar_c <- as(ras, "Raster") 
class(rar_c)
#> [1] "RasterLayer"
#> attr(,"package")
#> [1] "raster"
```

---
# Basic cell-based computations with [`stars`](https://CRAN.R-project.org/package=stars)

```r
class(ras[[1]])
#> [1] "units"
class(ras[["bathy.tif"]])
#> [1] "units"
class(ras$bathy.tif)
#> [1] "units"
```

---
# Use units with [`units`](https://CRAN.R-project.org/package=units)

```r
library(units)
as_units(1, "mm") + as_units(1, "cm")
#> 11 [mm]
```

```r
units::as_units(0, "celsius") + units::as_units(0, "fahrenheit")
#> -17.77778 [°C]
```

---
# Basic cell-based computations with [`stars`](https://CRAN.R-project.org/package=stars)

```r
ras_nounit <- drop_units(ras)
class(ras_nounit[[1]])
#> [1] "matrix" "array"
```

---
# Basic cell-based computations with [`stars`](https://CRAN.R-project.org/package=stars)

```r
# extract value
ras[[1]][1, 1]
#> 457 [m]
ras[[1]][1:10, 11:20]
#> Units: [m]
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#>  [1,]  455  475  479  476  474  473  488  529  525   521
#>  [2,]  453  472  473  471  476  474  494  533  533   518
#>  [3,]  454  471  471  471  477  475  493  539  536   515
#>  [4,]  454  469  470  472  472  473  487  530  528   511
#>  [5,]  456  469  470  475  472  476  493  519  520   508
#>  [6,]  459  474  471  481  478  488  508  523  525   507
#>  [7,]  459  484  473  494  487  502  519  537  537   513
#>  [8,]  467  497  479  514  495  502  525  554  542   522
#>  [9,]  483  506  483  526  497  495  531  554  535   525
#> [10,]  508  507  485  520  490  485  530  543  527   513
```

---
# Basic cell-based computations with [`stars`](https://CRAN.R-project.org/package=stars)

```r
# mean along latitude
m_lat <- apply(ras[[1]], 2, mean)
m_lat[1:10]
#>  [1] 294.2492 295.1260 295.9008 297.1172 296.7591 295.6539 294.2279 295.1870
#>  [9] 296.0227 296.5232
# mean along longitude
m_lon <- apply(ras[[1]], 1, mean)
m_lon[1:10]
#>  [1] 478.7497 479.6403 480.5201 480.8931 480.7000 480.7611 480.9503 481.4774
#>  [9] 482.5569 483.8781
```

---
# Basic cell-based computations with [`stars`](https://CRAN.R-project.org/package=stars)

```r
lon <- st_get_dimension_values(ras, 1)
lat <- st_get_dimension_values(ras, 2)
```

---
# Basic cell-based computations with [`stars`](https://CRAN.R-project.org/package=stars)

```r
plot(lat, m_lat)
```

---
# Basic cell-based computations with [`stars`](https://CRAN.R-project.org/package=stars)

```r
range(ras[[1]])  # min/max
#> Units: [m]
#> [1] -4490  1502
diff(range(ras[[1]]))
#> 5992 [m]
mean(ras[[1]]) 
#> 59.65972 [m]
quantile(ras[[1]])
#> Units: [m]
#>    0%   25%   50%   75%  100% 
#> -4490   -94   114   439  1502
```

---
# Basic plot with [`stars`](https://CRAN.R-project.org/package=stars)

```r
plot(ras)
```

---
# Basic plot with [`stars`](https://CRAN.R-project.org/package=stars)

```r
par(mfrow = c(1, 2))
hist(ras[[1]])
plot(density(ras[[1]]))
```

---
# Available drivers

### - GeoTIFF <svg viewBox="0 0 512 512" style="height:1em;fill:currentColor;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M173.898 439.404l-166.4-166.4c-9.997-9.997-9.997-26.206 0-36.204l36.203-36.204c9.997-9.998 26.207-9.998 36.204 0L192 312.69 432.095 72.596c9.997-9.997 26.207-9.997 36.204 0l36.203 36.204c9.997 9.997 9.997 26.206 0 36.204l-294.4 294.401c-9.998 9.997-26.207 9.997-36.204-.001z"></path></svg>

### - Other formats, how do we know?

```r
st_dr <- sf::st_drivers(what = "raster")
head(st_dr)
#>            name                                  long_name write  copy
#> VRT         VRT                             Virtual Raster  TRUE  TRUE
#> DERIVED DERIVED Derived datasets using VRT pixel functions FALSE FALSE
#> GTiff     GTiff                                    GeoTIFF  TRUE  TRUE
#> NITF       NITF       National Imagery Transmission Format  TRUE  TRUE
#> RPFTOC   RPFTOC           Raster Product Format TOC format FALSE FALSE
#> ECRGTOC ECRGTOC                            ECRG TOC format FALSE FALSE
#>         is_raster is_vector   vsi
#> VRT          TRUE     FALSE  TRUE
#> DERIVED      TRUE     FALSE FALSE
#> GTiff        TRUE     FALSE  TRUE
#> NITF         TRUE     FALSE  TRUE
#> RPFTOC       TRUE     FALSE  TRUE
#> ECRGTOC      TRUE     FALSE  TRUE
```

---
# Check read/write capabilities with

```r
st_dr[which(st_dr$name == "GTiff"), ]
#>        name long_name write copy is_raster is_vector  vsi
#> GTiff GTiff   GeoTIFF  TRUE TRUE      TRUE     FALSE TRUE
```

--
.right[

.font80[<svg viewBox="0 0 512 512" style="height:1em;fill:currentColor;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M256 8C119.043 8 8 119.083 8 256c0 136.997 111.043 248 248 248s248-111.003 248-248C504 119.083 392.957 8 256 8zm0 110c23.196 0 42 18.804 42 42s-18.804 42-42 42-42-18.804-42-42 18.804-42 42-42zm56 254c0 6.627-5.373 12-12 12h-88c-6.627 0-12-5.373-12-12v-24c0-6.627 5.373-12 12-12h12v-64h-12c-6.627 0-12-5.373-12-12v-24c0-6.627 5.373-12 12-12h64c6.627 0 12 5.373 12 12v100h12c6.627 0 12 5.373 12 12v24z"></path></svg> https://gdal.org/drivers/raster/index.html]

]

---
# Write raster files

##  .center[.font200[<svg viewBox="0 0 581 512" style="height:1em;fill:currentColor;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M581 226.6C581 119.1 450.9 32 290.5 32S0 119.1 0 226.6C0 322.4 103.3 402 239.4 418.1V480h99.1v-61.5c24.3-2.7 47.6-7.4 69.4-13.9L448 480h112l-67.4-113.7c54.5-35.4 88.4-84.9 88.4-139.7zm-466.8 14.5c0-73.5 98.9-133 220.8-133s211.9 40.7 211.9 133c0 50.1-26.5 85-70.3 106.4-2.4-1.6-4.7-2.9-6.4-3.7-10.2-5.2-27.8-10.5-27.8-10.5s86.6-6.4 86.6-92.7-90.6-87.9-90.6-87.9h-199V361c-74.1-21.5-125.2-67.1-125.2-119.9zm225.1 38.3v-55.6c57.8 0 87.8-6.8 87.8 27.3 0 36.5-38.2 28.3-87.8 28.3zm-.9 72.5H365c10.8 0 18.9 11.7 24 19.2-16.1 1.9-33 2.8-50.6 2.9v-22.1z"></path></svg> <svg viewBox="0 0 448 512" style="height:1em;fill:currentColor;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M190.5 66.9l22.2-22.2c9.4-9.4 24.6-9.4 33.9 0L441 239c9.4 9.4 9.4 24.6 0 33.9L246.6 467.3c-9.4 9.4-24.6 9.4-33.9 0l-22.2-22.2c-9.5-9.5-9.3-25 .4-34.3L311.4 296H24c-13.3 0-24-10.7-24-24v-32c0-13.3 10.7-24 24-24h287.4L190.9 101.2c-9.8-9.3-10-24.8-.4-34.3z"></path></svg> <svg viewBox="0 0 384 512" style="height:1em;fill:currentColor;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M369.9 97.9L286 14C277 5 264.8-.1 252.1-.1H48C21.5 0 0 21.5 0 48v416c0 26.5 21.5 48 48 48h288c26.5 0 48-21.5 48-48V131.9c0-12.7-5.1-25-14.1-34zM332.1 128H256V51.9l76.1 76.1zM48 464V48h160v104c0 13.3 10.7 24 24 24h104v288H48z"></path></svg>]]

## - Share the output file

## - Use another software (e.g. QGIS)

<br>

### <svg viewBox="0 0 512 512" style="height:1em;fill:currentColor;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M504 256c0 136.997-111.043 248-248 248S8 392.997 8 256C8 119.083 119.043 8 256 8s248 111.083 248 248zm-248 50c-25.405 0-46 20.595-46 46s20.595 46 46 46 46-20.595 46-46-20.595-46-46-46zm-43.673-165.346l7.418 136c.347 6.364 5.609 11.346 11.982 11.346h48.546c6.373 0 11.635-4.982 11.982-11.346l7.418-136c.375-6.874-5.098-12.654-11.982-12.654h-63.383c-6.884 0-12.356 5.78-11.981 12.654z"></path></svg> Create an `output` directory

```r
dir.create("output", showWarnings = FALSE)
```

---
# Write with [`stars`](https://CRAN.R-project.org/package=stars)

```r
sort(st_dr$name[st_dr$write])
#>  [1] "ADRG"      "BMP"       "BT"        "BYN"       "CTable2"   "EHdr"     
#>  [7] "ELAS"      "ENVI"      "ERS"       "FITS"      "GPKG"      "GS7BG"    
#> [13] "GSBG"      "GTiff"     "GTX"       "HDF4Image" "HFA"       "IDA"      
#> [19] "ILWIS"     "INGR"      "ISCE"      "ISIS2"     "ISIS3"     "KRO"      
#> [25] "LAN"       "Leveller"  "MBTiles"   "MEM"       "MFF"       "MFF2"     
#> [31] "MRF"       "netCDF"    "NGW"       "NITF"      "NTv2"      "NWT_GRD"  
#> [37] "PAux"      "PCIDSK"    "PCRaster"  "PDF"       "PDS4"      "PNM"      
#> [43] "RMF"       "ROI_PAC"   "RRASTER"   "RST"       "SAGA"      "SGI"      
#> [49] "Terragen"  "VRT"
```

---
# Write with [`stars`](https://CRAN.R-project.org/package=stars)

```r
write_stars(ras, dsn = "output/ras.gpkg", driver = "GPKG")
```

---
# <svg viewBox="0 0 512 512" style="height:1em;fill:currentColor;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M243.2 189.9V258c26.1 5.9 49.3 15.6 73.6 22.3v-68.2c-26-5.8-49.4-15.5-73.6-22.2zm223.3-123c-34.3 15.9-76.5 31.9-117 31.9C296 98.8 251.7 64 184.3 64c-25 0-47.3 4.4-68 12 2.8-7.3 4.1-15.2 3.6-23.6C118.1 24 94.8 1.2 66.3 0 34.3-1.3 8 24.3 8 56c0 19 9.5 35.8 24 45.9V488c0 13.3 10.7 24 24 24h16c13.3 0 24-10.7 24-24v-94.4c28.3-12.1 63.6-22.1 114.4-22.1 53.6 0 97.8 34.8 165.2 34.8 48.2 0 86.7-16.3 122.5-40.9 8.7-6 13.8-15.8 13.8-26.4V95.9c.1-23.3-24.2-38.8-45.4-29zM169.6 325.5c-25.8 2.7-50 8.2-73.6 16.6v-70.5c26.2-9.3 47.5-15 73.6-17.4zM464 191c-23.6 9.8-46.3 19.5-73.6 23.9V286c24.8-3.4 51.4-11.8 73.6-26v70.5c-25.1 16.1-48.5 24.7-73.6 27.1V286c-27 3.7-47.9 1.5-73.6-5.6v67.4c-23.9-7.4-47.3-16.7-73.6-21.3V258c-19.7-4.4-40.8-6.8-73.6-3.8v-70c-22.4 3.1-44.6 10.2-73.6 20.9v-70.5c33.2-12.2 50.1-19.8 73.6-22v71.6c27-3.7 48.4-1.3 73.6 5.7v-67.4c23.7 7.4 47.2 16.7 73.6 21.3v68.4c23.7 5.3 47.6 6.9 73.6 2.7V143c27-4.8 52.3-13.6 73.6-22.5z"></path></svg> Recap

| Action            | *raster*                 | *stars*               |
|:------------------|:-------------------------|:----------------------|
| list drivers      | `rgdal::gdalDrivers()`   | `sf::st_drivers()`    |
| read files        | `raster()`               | `read_stars()`        |
| write files       | `writeRaster()`          | `write_stars()`       |

---
# Create a raster from scratch: Matrix

```r
M1 <- matrix(runif(36), 6, 6)
M1
#>              [,1]       [,2]      [,3]      [,4]       [,5]      [,6]
#> [1,] 0.4133132407 0.50870563 0.5326004 0.4920551 0.79999252 0.5981404
#> [2,] 0.0007539436 0.92675278 0.9368617 0.4082894 0.71711866 0.4173992
#> [3,] 0.0477994368 0.08819536 0.1527133 0.6701650 0.11999926 0.6179882
#> [4,] 0.3670485087 0.28122962 0.4955913 0.1019189 0.30700703 0.7223401
#> [5,] 0.0418033798 0.34368128 0.4580708 0.7249584 0.07742749 0.7319797
#> [6,] 0.0844074509 0.62711792 0.8952811 0.5294131 0.88966825 0.1493172
```

---
# Create a raster from scratch: Matrix

```r
ras_M1 <- st_as_stars(M1) 
ras_M1
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#>       A1            
#>  Min.   :0.0007539  
#>  1st Qu.:0.1518642  
#>  Median :0.4750630  
#>  Mean   :0.4521418  
#>  3rd Qu.:0.6819034  
#>  Max.   :0.9368617  
#> dimension(s):
#>    from to offset delta refsys point values x/y
#> X1    1  6      0     1     NA FALSE   NULL [x]
#> X2    1  6      0     1     NA FALSE   NULL [y]
st_crs(ras_M1)
#> Coordinate Reference System: NA
```

---
# Create a raster from scratch: Matrix

```r
plot(ras_M1)
```

---
# Create a raster from scratch: Matrix

```r
st_crs(ras_M1) <- st_crs(4326)
st_crs(ras_M1)
#> Coordinate Reference System:
#>   User input: EPSG:4326 
#>   wkt:
#> GEOGCRS["WGS 84",
#>     DATUM["World Geodetic System 1984",
#>         ELLIPSOID["WGS 84",6378137,298.257223563,
#>             LENGTHUNIT["metre",1]]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433]],
#>     CS[ellipsoidal,2],
#>         AXIS["geodetic latitude (Lat)",north,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         AXIS["geodetic longitude (Lon)",east,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>     USAGE[
#>         SCOPE["unknown"],
#>         AREA["World"],
#>         BBOX[-90,-180,90,180]],
#>     ID["EPSG",4326]]
```

---
# Create a raster from scratch: Data Frame

```r
D1 <- expand.grid(lon = 1:10, lat = 1:10) 
D1$val <- runif(100)
head(D1)
#>   lon lat       val
#> 1   1   1 0.5675234
#> 2   2   1 0.7456320
#> 3   3   1 0.7996474
#> 4   4   1 0.2249119
#> 5   5   1 0.5699844
#> 6   6   1 0.1766411
```

---
# Create a raster from scratch: Data Frame

```r
ras_D1 <- st_as_stars(D1, coords = c("lon", "lat"))
st_crs(ras_D1) <- st_crs(4326)
ras_D1
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#>       val          
#>  Min.   :0.005091  
#>  1st Qu.:0.266516  
#>  Median :0.440904  
#>  Mean   :0.473069  
#>  3rd Qu.:0.673801  
#>  Max.   :0.992526  
#> dimension(s):
#>     from to offset delta refsys point values x/y
#> lon    1 10    0.5     1 WGS 84    NA   NULL [x]
#> lat    1 10   10.5    -1 WGS 84    NA   NULL [y]
```

* See `?st_as_stars` for more details.

---

![:custom_hr]()

---
# Create a modified copy

```r
ras2 <- 2*ras 
ras2
#> stars object with 2 dimensions and 1 attribute
#> attribute(s), summary of first 1e+05 cells:
#>  bathy.tif [m] 
#>  Min.   :-870  
#>  1st Qu.:  46  
#>  Median : 958  
#>  Mean   : 599  
#>  3rd Qu.:1080  
#>  Max.   :1624  
#> dimension(s):
#>   from   to offset       delta refsys point values x/y
#> x    1 3840    -71  0.00416667 WGS 84 FALSE   NULL [x]
#> y    1 2880     55 -0.00416667 WGS 84 FALSE   NULL [y]
```

---
# Simple cell-based filter with [`stars`](https://CRAN.R-project.org/package=stars)

```r
ras_f1 <- ras # create a copy
ras_f1[[1]][ras[[1]] > units::as_units(0, "m")] <- NA # filter
plot(ras_f1)
```

---
# Simple cell-based filter with [`stars`](https://CRAN.R-project.org/package=stars)

```r
ras_f2 <- units::drop_units(ras) # create a copy without unit
ras_f2[[1]][ras_f2[[1]] < 0] <- NA # filter
plot(ras_f2)
```

---
# Change map projection

- Reminder

- [<svg viewBox="0 0 640 512" style="height:1em;fill:currentColor;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M640 51.2l-.3 12.2c-28.1.8-45 15.8-55.8 40.3-25 57.8-103.3 240-155.3 358.6H415l-81.9-193.1c-32.5 63.6-68.3 130-99.2 193.1-.3.3-15 0-15-.3C172 352.3 122.8 243.4 75.8 133.4 64.4 106.7 26.4 63.4.2 63.7c0-3.1-.3-10-.3-14.2h161.9v13.9c-19.2 1.1-52.8 13.3-43.3 34.2 21.9 49.7 103.6 240.3 125.6 288.6 15-29.7 57.8-109.2 75.3-142.8-13.9-28.3-58.6-133.9-72.8-160-9.7-17.8-36.1-19.4-55.8-19.7V49.8l142.5.3v13.1c-19.4.6-38.1 7.8-29.4 26.1 18.9 40 30.6 68.1 48.1 104.7 5.6-10.8 34.7-69.4 48.1-100.8 8.9-20.6-3.9-28.6-38.6-29.4.3-3.6 0-10.3.3-13.6 44.4-.3 111.1-.3 123.1-.6v13.6c-22.5.8-45.8 12.8-58.1 31.7l-59.2 122.8c6.4 16.1 63.3 142.8 69.2 156.7L559.2 91.8c-8.6-23.1-36.4-28.1-47.2-28.3V49.6l127.8 1.1.2.5z"></path></svg> List of map projections](https://en.wikipedia.org/wiki/List_of_map_projections)

- [Interactive examples](https://mathigon.org/course/circles/spheres-cones-cylinders#sphere-maps)

- Let's use the spherical Mercator 
  - https://epsg.io/
  - https://epsg.io/3857

---
# Change map projection with [`stars`](https://CRAN.R-project.org/package=stars)

### Just as it is done with `sf`

```r
ras_t <- st_transform(ras, crs = 3857)
st_crs(ras_t)
#> Coordinate Reference System:
#>   User input: EPSG:3857 
#>   wkt:
#> PROJCRS["WGS 84 / Pseudo-Mercator",
#>     BASEGEOGCRS["WGS 84",
#>         DATUM["World Geodetic System 1984",
#>             ELLIPSOID["WGS 84",6378137,298.257223563,
#>                 LENGTHUNIT["metre",1]]],
#>         PRIMEM["Greenwich",0,
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         ID["EPSG",4326]],
#>     CONVERSION["Popular Visualisation Pseudo-Mercator",
#>         METHOD["Popular Visualisation Pseudo Mercator",
#>             ID["EPSG",1024]],
#>         PARAMETER["Latitude of natural origin",0,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8801]],
#>         PARAMETER["Longitude of natural origin",0,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8802]],
#>         PARAMETER["False easting",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8806]],
#>         PARAMETER["False northing",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8807]]],
#>     CS[Cartesian,2],
#>         AXIS["easting (X)",east,
#>             ORDER[1],
#>             LENGTHUNIT["metre",1]],
#>         AXIS["northing (Y)",north,
#>             ORDER[2],
#>             LENGTHUNIT["metre",1]],
#>     USAGE[
#>         SCOPE["unknown"],
#>         AREA["World - 85°S to 85°N"],
#>         BBOX[-85.06,-180,85.06,180]],
#>     ID["EPSG",3857]]
```

---
# Crop with [`stars`](https://CRAN.R-project.org/package=raster)

---
# Crop with [`stars`](https://CRAN.R-project.org/package=raster)

```r
ext <- st_bbox(c(xmin = -65, xmax = -60, ymin = 45, ymax = 50), crs = 4326)
ras_c <- st_crop(ras, ext)
plot(ras_c)
```

---
# Mask

Let's use a shapefile of the Gulf of the St-Lawrence

```r
(stl <- sf::st_read("data/st_laurence.geojson"))
#> Reading layer `data_stl2' from data source `/__w/wksh_raster_CWS/wksh_raster_CWS/data/st_laurence.geojson' using driver `GeoJSON'
#> Simple feature collection with 1 feature and 10 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -74.86481 ymin: 44.9585 xmax: -54.70345 ymax: 52.22242
#> Geodetic CRS:  WGS 84
#> Simple feature collection with 1 feature and 10 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -74.86481 ymin: 44.9585 xmax: -54.70345 ymax: 52.22242
#> Geodetic CRS:  WGS 84
#>                   name id longitude latitude     min_x   min_y     max_x
#> 1 Gulf of St. Lawrence 24 -61.56657 48.38072 -74.86481 44.9585 -54.70345
#>      max_y   area mrgid                       geometry
#> 1 52.22242 290874  4290 POLYGON ((-64.67242 47.7168...
```

---
# Mask

Values outside of a given geometry set to `NA`

```r
plot(sf::st_geometry(stl))
```

---
# Mask with [`stars`](https://CRAN.R-project.org/package=stars)

```r
ras_m <- ras[stl]
plot(ras_m)
```

---
# Mask with [`stars`](https://CRAN.R-project.org/package=stars)

```r
stl_i <- sf::st_difference(sf::st_as_sfc(st_bbox(ras)), stl)
plot(stl_i, col = 2)
```

---
# Mask with [`stars`](https://CRAN.R-project.org/package=stars)

```r
ras_mi <- ras[stl_i]
plot(ras_mi)
```

---
# Warping/Resampling with [`stars`](https://CRAN.R-project.org/package=stars)

```r
ras_template1 <- st_as_stars(st_bbox(ras), nx = 21, ny = 21, 
  values = runif(21 * 21)) # create template
ras_w1 <- st_warp(ras, ras_template1) 
plot(ras_w1)
```

???

Work with rasters of different resolutions

---
# Warping/Resampling with [`stars`](https://CRAN.R-project.org/package=stars)

```r
ras_template1 <- st_as_stars(st_bbox(ras), nx = 21, ny = 21, 
  values = runif(21 * 21)) # create template
ras_w1 <- st_warp(ras, st_transform(ras_template1, crs = 3857)) 
plot(ras_w1)
```

---
# Warping/Resampling with [`stars`](https://CRAN.R-project.org/package=stars)

```r
ras_w2 <- st_warp(ras, cellsize = 0.25, crs = st_crs(ras)) 
plot(ras_w2)
```

---
# Rasterize with [`stars`](https://CRAN.R-project.org/package=stars)

```r
ras_template1[[1]] <- NA
stl_r <- st_rasterize(stl, ras_template1) # use a template
plot(stl_r)
```

---
# Rasterize with [`stars`](https://CRAN.R-project.org/package=stars)

```r
stl_s <- st_rasterize(stl, dy = .1, dx = .1) # use cell-size
plot(stl_s)
```

---
# Rasterize faster

- Checkout [fasterize](https://github.com/ecohealthalliance/fasterize)

- Check out [terra](https://cran.r-project.org/web/packages/terra/index.html)

---
# Mask faster with `st_rasterize()`

### Solution proposed by Dewey Dunnington

```r
# load files
gulf_region <- sf::read_sf("data/st_laurence.geojson")
strs <- stars::read_stars("data/bathy.tif")
# create raster mask
template <- strs
template[[1]][] <- NA
gulf_region$val_ras <- 1

rasterized <- stars::st_rasterize(gulf_region["val_ras"], template = template)
```

---
# Mask faster  with `st_rasterize()`

### Solution proposed by Dewey Dunnington

```r
plot(rasterized)
```

---
# Polygonize

```r
ras_p <- units::drop_units(ras)
ras_p[[1]][ras_p[[1]] > 0] <- NA # filter
ras_p[[1]][ras_p[[1]] < 0] <- 1
pol <- st_as_sf(ras_p, as_points = FALSE, merge = TRUE)
```

---
# Polygonize

```r
plot(pol)  
```

---
# Polygonize

```r
plot(pol[which.max(st_area(pol)),])  
```

---
# Dealing with large rasters

- When a raster is large (takes fair amount of your RAM), `stars` creates a `stars_proxy` object.

> Stars proxy objects take another approach: upon creation they contain no data at all, but only pointers to where the data can be read. Data is only read when it is needed, and only as much as is needed [...]. ([Stars proxy objects](https://r-spatial.github.io/stars/articles/stars2.html))

- You might want to try out [terra](https://cran.r-project.org/web/packages/terra/index.html).

---
# Stack

### - Collection of rasters

### - Same extent/resolution

### - Multi-band rasters (satellite images)

### - Handle several rasters at once

### - `rasterStack` and `rasterBrick` in package `raster`