-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathspatial_raster.qmd
125 lines (84 loc) · 4.65 KB
/
spatial_raster.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
---
title: "Raster Data"
engine: knitr
---
## Library Loading
Before we can dive into "actual" spatial work we'll need to load some libraries.
:::panel-tabset
## [{{< fa brands r-project >}} R]{.r}
Load `terra` to work with raster data.
```{r r-libs, warning = F, message = F}
# Load needed libraries
library(terra)
```
## [{{< fa brands python >}} Python]{.py}
Load the `rioxarray` and `rasterio` libraries to work with raster data. We'll also load the `os` library to deal with file path issues.
```{python py-libs}
# Load needed libraries
import os
import rasterio as rio
import rioxarray as rxr
```
:::
To demonstrate raster data operations we'll use NASA's [Shuttle Radar Topography Mission Global 3 arc second](https://lpdaac.usgs.gov/products/srtmgl3v003/) elevation data. Note that some minor preparatory work was necessary to get the data ready for our purposes here and is preserved in [this folder](https://github.com/njlyon0/collab_bilingualism/tree/main/dev) of the website's {{< fa brands github >}} GitHub repository.
## Loading Rasters
We can begin by reading in the data.
:::panel-tabset
## [{{< fa brands r-project >}} R]{.r}
There are several [{{< fa brands r-project >}} R]{.r} packages for working with raster data but we'll focus on `terra`.
```{r r-rast-read}
# Load a raster of elevation data
rast_r <- terra::rast(x = file.path("data", "elevation.tif"))
# Check the class of this object
class(rast_r)
```
## [{{< fa brands python >}} Python]{.py}
The `rioxarray` library will be our focal library for raster operations with [{{< fa brands python >}} Python]{.py}.
```{python py-rast-read}
# Load a raster of elevation data
rast_py = rxr.open_rasterio(os.path.join("data", "elevation.tif"), masked = True).squeeze()
# Check the type of this variable
type(rast_py.rio.crs)
```
:::
## Checking Raster CRS
As noted earlier, the very first thing we should do after reading in spatial data is _check the coordinate reference system!_
:::panel-tabset
## [{{< fa brands r-project >}} R]{.r}
We can check the CRS with the `crs` function (from `terra`).
```{r r-crs-check}
# Check the CRS of the elevation data
terra::crs(x = rast_r)
```
## [{{< fa brands python >}} Python]{.py}
[{{< fa brands python >}} Python]{.py} stores raster CRS as an attribute.
```{python py-crs-check}
# Check the CRS of the elevation data
rast_py.rio.crs
```
:::
You may note in the above operations that [{{< fa brands python >}} Python]{.py} and [{{< fa brands r-project >}} R]{.r} seem to be returning different information for the same data. In fact they are saying the same thing just in slightly different ways. [{{< fa brands r-project >}} R]{.r} is telling us the name of the coordinate reference system (World Geodetic System 1984) while [{{< fa brands python >}} Python]{.py} is giving us the four-digit "EPSG" code. For context, the European Petroleum Survey Group (EPSG) is a major authority on accepted coordinate reference systems. Each CRS is given a unique, four-digit code. As you may now be able to guess, "4326" is the EPSG code for the World Geodetic System 1984 CRS!
This CRS (WGS84) is an especially common one so you may eventually commit its EPSG code to memory but often you'll encounter either CRS names or EPSG codes with which you are not familiar. It may not sound like sage advice but it can be simplest to just use whichever piece of information your coding language returns to Google the "missing" information.
## Transforming Raster CRS
Once we know the current CRS, we can transform into a different coordinate reference system as needed. Let's transform from WGS84 into another commonly-used CRS: Albers Equal Area (EPSG code 3083).
Quick warning, transforming rasters can take a **long** time and is very computationally intensive. It'll work for a relatively small raster like the one we're working with here but in general you should be careful with attempting to re-project a raster into a different CRS.
:::panel-tabset
## [{{< fa brands r-project >}} R]{.r}
For CRS transformations in [{{< fa brands r-project >}} R]{.r} we can use the `project` function (pronounced like the verb not the noun).
```{r r-crs-change}
# Transform the raster CRS
rast_alber_r <- terra::project(x = rast_r, y = "epsg:3083")
# Re-check the CRS to confirm it worked
terra::crs(rast_alber_r)
```
## [{{< fa brands python >}} Python]{.py}
[{{< fa brands python >}} Python]{.py} CRS re-projections have slightly more steps to work through.
```{python py-crs-change}
# Create a CRS variable for the desired ending CRS
wgs84_crs = rio.crs.CRS.from_string("EPSG:3083")
# Reproject the raster we had
rast_alber_py = rast_py.rio.reproject(wgs84_crs)
# Re-check the CRS to confirm it worked
rast_alber_py.rio.crs
```
:::