generated from jtr13/quarto-edav-template
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathch3_a.qmd
279 lines (192 loc) · 10.9 KB
/
ch3_a.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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
---
title-block-banner: true
bibliography: references.bib
editor:
markdown:
wrap: sentence
---
## Raster Arithmetic Operations
Arithmetic, comparison and logical operations on `SpatRasters` use the same operators as those used for simple vector-like operations.
These operators are listed below:
**a) Arithmetic operators**
{fig-align="center" width="336"}
**b) Comparison operators**
{fig-align="center" width="331"}
**c) Logical operators**
{fig-align="center" width="331"}
------------------------------------------------------------------------
### Single Raster Operations
The application of these operators on rasters can be conceptualized as element-wise operations on a matrix.
For example, assume ***A*** is a raster with values shown in the matrix form as below:
{width="300"}
Then, ***A*** + 1 will be:
{width="400"}
Note how each element of the matrix gets an addition of 1.
Similarly, for ***A*** x 2, each element of ***A*** will be multiplied by 2 as:
{width="400"}
and, ***A*** \^ 2 will be given as:
{width="400"}
Logical operations can also be carried out in a similar fashion.
For example, if we were to test if the cell values of ***A*** are greater than 3 (i.e. ***A*** \> 3), it will be written as:
{width="400"}
Now let us try a few examples of arithmetic, comparison and logical operations on `SpatRasters`.
Before we begin this chapter, let us ensure that we have the sample dataset necessary for the analysis, otherwise, download the files from the GitHub repository using the following codes:
```{r, message = FALSE, warning = FALSE}
#~~~ Check if sample data exists on disk. If not, download from GitHub repository
if (dir.exists("SampleData-master")==FALSE){
download.file(url = "https://github.com/Vinit-Sehgal/SampleData/archive/master.zip",
destfile = "SampleData-master.zip")
# Unzip the downloaded .zip file
unzip(zipfile = "SampleData-master.zip")
list.files("./SampleData-master") # List folder contents
}
```
Taking example of SMAP soil moisture raster let's practice the application of arithmetic operations on `SpatRasters`.
We start by importing the `SpatRaster` in the current R environment.
```{r, message = FALSE, warning = FALSE}
library(terra) # Import library
# Import SMAP soil moisture raster from the downloaded folder
sm=terra::rast("./SampleData-master/raster_files/SMAP_SM.tif")
print(sm)
# Add 1 to raster values
sm2=sm+1
print(sm2) #Notice the max and min values have increased by 1.
# Multiply raster values by 2
sm2=sm*2
plot(sm2, main= "sm*2") # Try sm2=sm*10, or sm2=sm^2 and see the difference in sm2 values
```
Similarly, logical operations can also be carried out on rasters with ease in R.
Logical operators can be a convenient tool for raster data manipulation.
For example, any subset of values can be manipulated based on a user-defined logical criteria.
```{r, message = FALSE, warning = FALSE}
# Are cell values of sm raster greater than 0.3?
sm=terra::rast("./SampleData-master/raster_files/SMAP_SM.tif")
sm2=sm>0.3
plot(sm2, main="is cell value >0.3") # Notice that the raster has only True and False values
# Replace all values of sm>0.3 with 0.5
sm=terra::rast("./SampleData-master/raster_files/SMAP_SM.tif")
sm[sm>0.3]=0.5
plot(sm, main="replace sm>0.3 with 0.5")
# Or replace all values of sm>0.3 with NA
sm=terra::rast("./SampleData-master/raster_files/SMAP_SM.tif")
sm[sm>0.3]=NA
plot(sm, main="replace sm>0.3 with NA")
# Reaplace all NA values with 0
sm=terra::rast("./SampleData-master/raster_files/SMAP_SM.tif")
sm[is.na(sm)]=0
plot(sm)
```
Satellites record observations in regularly spaced swaths as they orbit around the globe.
In the case of SMAP satellite, it takes 2-3 days to cover the Earth.
So far, we used a pre-processed soil moisture raster that provided 3-day averaged soil moisture using rasters from individual days.
Now, let us evaluate daily raster available from SMAP (pixel recorded at 6 AM local time).
Notice that this raster is patchy and uses a fill value (values assigned to missing observations) of -9999.
Replacing these fill values with `NA` prior to statistical analysis are highly useful.
```{r, message = FALSE, warning = FALSE}
# Reaplace fill values with NA
sm=rast("./SampleData-master/raster_files/SMAP_L3_SM_P_20150401_R17000_002_Soil_Moisture_Retrieval_Data_AM_soil_moisture_867a5e36.tif")
smFilter=sm # Copy original raster before manipulation
smFilter[smFilter==-9999]=NA # Replace negative values with -9999
# Plot original and manipulated raster
layout(matrix(1:2, ncol = 2)) # Specify layout for plots
plot(sm, main= "raw raster with fill values")
plot(smFilter, main="Post filtering -9999")
```
------------------------------------------------------------------------
### Multi-Raster Operations
So far we saw examples of arithmetic/ logical operations on a single raster.
The tenets of arithmetic operations discussed earlier are also applicable for operations on two (or more) rasters.
These operations are carried out cell-wise between the corresponding pixels of the raster layers.
Suppose, another raster ***B*** is given as:
{width="300"}
then, ***A***+***B*** will be:
{width="500"}
Similarly, ***A***x***B*** will also be a cell-wise multiplication of respective data values:
{width="500"}
For this illustration, we will use rootzone (0-100 cm depth) soil moisture from NASA's Catchment Land Surface Model (provided as SMAP science Level 4 product) for Contiguous U.S. over a period of three consecutive days i.e.
December 13-15, 2021.
```{r, message = FALSE, warning = FALSE}
SMday1=rast("./SampleData-master/SMAPL4_rasters/SMAP_L4_SM_gph_20211213T103000_Vv6030_001_Geophysical_Data_sm_rootzone_2052db45.tif")
SMday2=rast("./SampleData-master/SMAPL4_rasters/SMAP_L4_SM_gph_20211214T103000_Vv6030_001_Geophysical_Data_sm_rootzone_22156243.tif")
SMday3=rast("./SampleData-master/SMAPL4_rasters/SMAP_L4_SM_gph_20211215T103000_Vv6030_001_Geophysical_Data_sm_rootzone_23bd7fe2.tif")
layout(matrix(1:3, ncol = 1)) # Specify layout
plot(SMday1, main="Rootzone SM day1")
plot(SMday2, main="Rootzone SM day2")
plot(SMday3, main="Rootzone SM day3")
```
Note that we can not make out the differences in the rasters as the colors are skewed by the high negative fill value for missing data in the rasters.
So, we first remove the fill values (i.e. cells valued at -9999) with `NA`.
```{r, message = FALSE, warning = FALSE}
SMday1[SMday1==-9999]=NA # Replace negative values with -9999
SMday2[SMday2==-9999]=NA # Replace negative values with -9999
SMday3[SMday3==-9999]=NA # Replace negative values with -9999
layout(matrix(1:3, ncol = 1)) # Specify layout
plot(SMday1, main="SM day1, filtered")
plot(SMday2, main="SM day2, filtered")
plot(SMday3, main="SM day3, filtered")
```
Now, let's calculate the change in soil moisture values between days 1-2 and 1-3.
Can we identify a region which showed anomalous wetting during the three days of observation?
```{r, message = FALSE, warning = FALSE}
delta12=SMday1-SMday2
delta13=SMday1-SMday3
layout(matrix(1:2, ncol = 1)) # Specify layout
plot(delta12, main="SM difference (day 1 & 2")
plot(delta13, main="SM difference (day 1 & 3")
```
Curious to know why Southern California showed such high levels of wetting?
Check out: [']('Storm%20of%20the%20season'%20dumps%20record-breaking%20rainfall%20on%20SoCal%20and%20snow%20in%20the%20mountains%20-%20Los%20Angeles%20Times%20(latimes.com))[Storm of the season' dumps record-breaking rainfall on SoCal and snow in the mountains - Los Angeles Times (latimes.com)](https://www.latimes.com/california/story/2021-12-14/storm-of-the-season-delivers-rain-snow-and-wind-to-socal).
------------------------------------------------------------------------
## Raster Statistics
We are all familiar with common statistical functions available in base R for summarizing arrays, such as: `min`, `max`, `range`, `sum`, `stdev`, `median`, `mean`, `modal`.
These functions are applied using the `global` function, which calculates the "*global*" statistics based on the values of all cells within a raster layer.
`Global` is a versatile function and can also be used to apply more complex user-defined operations on rasters.
Let us now look at some examples:
```{r, message = FALSE, warning = FALSE}
library(terra)
# Import SMAP soil moisture raster from the downloaded folder
sm=terra::rast("./SampleData-master/raster_files/SMAP_SM.tif")
# Summary statistics
global(sm, mean, na.rm = T) # Calculate mean, while ignoring NA values
# This is equivalent to:
sm_val=values(sm) # Create array of cell values of raster
mean(sm_val, na.rm = T) # Calculate mean of cell values
```
Use of `global` function is equivalent to applying a function on all cell values
```{r, message = FALSE, warning = FALSE}
global(sm, sd, na.rm = T) # Calculate standard deviation, while ignoring NA values
# OR
sd(values(sm), na.rm = T) # Calculate mean of cell values
# Similarly, for finding quantiles of a raster layer
global(sm, quantile, probs = c(0.25, 0.75), na.rm = T) # Calculate 25th and 75th percentiles of the ratser layer
# OR
quantile(values(sm), probs = c(0.25, 0.75), na.rm = T)
```
We can also generate common summary statistics plots using functions such as hist (histogram), barplot (box and whisker plot), etc.
```{r, message = FALSE, warning = FALSE}
sm=terra::rast("./SampleData-master/raster_files/SMAP_SM.tif")
# Plot summary using standard statistical functions
layout(matrix(1:2, ncol = 2)) # Specify layout
hist(sm)
boxplot(sm)
```
### User-defined Functions for Raster Statistics
User-defined functions can also be applied to raster layer using the `global` function for specific operations.
A user can define their own function with specific instructions which will be executed everytime the function is called.
A simple template of a typical R function is given below:
{width="400"}
Let us revisit the application of the `quantile` function by passing it through a user-defined function with `global`.
In the `quant_fun` shown below function, `myRas` is the raster whose values are being summarized.
Notice how instruction for `na.rm` is passed through the `ignoreNA` argument.
```{r, message = FALSE, warning = FALSE}
# User-defined statistics by defining own function.
quant_fun = function(myRas, ignoreNA=TRUE){
p=quantile(myRas, probs = c(0.25, 0.75), na.rm=ignoreNA)
return(p)
}
global(sm, quant_fun) # 25th, and 75th percentile of each layer
# This is equivalent to:
quant_fun(values(sm))
```
------------------------------------------------------------------------
{fig-align="center" width="600"}