/
raster-vector.html
556 lines (537 loc) · 63.7 KB
/
raster-vector.html
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
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
<!DOCTYPE html>
<html lang="en">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1, shrink-to-fit=no">
<title>Chapter 6 Raster-vector interactions | Geocomputation with R</title>
<meta name="author" content="Robin Lovelace, Jakub Nowosad, Jannes Muenchow">
<meta name="description" content="Prerequisites This chapter requires the following packages: library(sf) library(terra) library(dplyr) 6.1 Introduction This Chapter focuses on interactions between raster and vector geographic...">
<meta name="generator" content="bookdown 0.39 with bs4_book()">
<meta property="og:title" content="Chapter 6 Raster-vector interactions | Geocomputation with R">
<meta property="og:type" content="book">
<meta property="og:url" content="https://r.geocompx.org/raster-vector.html">
<meta property="og:image" content="https://r.geocompx.org/images/cover.png">
<meta property="og:description" content="Prerequisites This chapter requires the following packages: library(sf) library(terra) library(dplyr) 6.1 Introduction This Chapter focuses on interactions between raster and vector geographic...">
<meta name="twitter:card" content="summary">
<meta name="twitter:title" content="Chapter 6 Raster-vector interactions | Geocomputation with R">
<meta name="twitter:description" content="Prerequisites This chapter requires the following packages: library(sf) library(terra) library(dplyr) 6.1 Introduction This Chapter focuses on interactions between raster and vector geographic...">
<meta name="twitter:image" content="https://r.geocompx.org/images/cover.png">
<!-- JS --><script src="https://cdnjs.cloudflare.com/ajax/libs/clipboard.js/2.0.6/clipboard.min.js" integrity="sha256-inc5kl9MA1hkeYUt+EC3BhlIgyp/2jDIyBLS6k3UxPI=" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/fuse.js/6.4.6/fuse.js" integrity="sha512-zv6Ywkjyktsohkbp9bb45V6tEMoWhzFzXis+LrMehmJZZSys19Yxf1dopHx7WzIKxr5tK2dVcYmaCk2uqdjF4A==" crossorigin="anonymous"></script><script src="https://kit.fontawesome.com/6ecbd6c532.js" crossorigin="anonymous"></script><script src="libs/jquery-3.6.0/jquery-3.6.0.min.js"></script><meta name="viewport" content="width=device-width, initial-scale=1, shrink-to-fit=no">
<link href="libs/bootstrap-4.6.0/bootstrap.min.css" rel="stylesheet">
<script src="libs/bootstrap-4.6.0/bootstrap.bundle.min.js"></script><link href="libs/Lato-0.4.9/font.css" rel="stylesheet">
<link href="libs/Roboto_Mono-0.4.9/font.css" rel="stylesheet">
<link href="libs/Montserrat-0.4.9/font.css" rel="stylesheet">
<script src="libs/bs3compat-0.7.0/transition.js"></script><script src="libs/bs3compat-0.7.0/tabs.js"></script><script src="libs/bs3compat-0.7.0/bs3compat.js"></script><link href="libs/bs4_book-1.0.0/bs4_book.css" rel="stylesheet">
<script src="libs/bs4_book-1.0.0/bs4_book.js"></script><meta name="citation_title" content="Chapter 6 Raster-vector interactions | Geocomputation with R">
<meta name="citation_author" content="Robin Lovelace">
<meta name="citation_author" content="Jakub Nowosad">
<meta name="citation_author" content="Jannes Muenchow">
<meta name="citation_publication_date" content="2019">
<meta name="citation_isbn" content="9780203730058">
<link href="libs/htmltools-fill-0.5.8.1/fill.css" rel="stylesheet">
<script src="libs/htmlwidgets-1.6.4/htmlwidgets.js"></script><link href="libs/leaflet-1.3.1/leaflet.css" rel="stylesheet">
<script src="libs/leaflet-1.3.1/leaflet.js"></script><link href="libs/leafletfix-1.0.0/leafletfix.css" rel="stylesheet">
<script src="libs/proj4-2.6.2/proj4.min.js"></script><script src="libs/Proj4Leaflet-1.0.1/proj4leaflet.js"></script><link href="libs/rstudio_leaflet-1.3.1/rstudio_leaflet.css" rel="stylesheet">
<script src="libs/leaflet-binding-2.2.2/leaflet.js"></script><script src="libs/kePrint-0.0.1/kePrint.js"></script><link href="libs/lightable-0.0.1/lightable.css" rel="stylesheet">
<link rel="icon" type="image/png" sizes="32x32" href="images/favicon-32x32.png">
<link rel="icon" type="image/png" sizes="16x16" href="images/favicon-16x16.png">
<script>
(function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){
(i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),
m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)
})(window,document,'script','https://www.google-analytics.com/analytics.js','ga');
ga('create', 'UA-99618359-1', 'auto');
ga('send', 'pageview');
</script><!-- Google tag (gtag.js) --><script async src="https://www.googletagmanager.com/gtag/js?id=G-VDC2S0ZNH5"></script><script>
window.dataLayer = window.dataLayer || [];
function gtag(){dataLayer.push(arguments);}
gtag('js', new Date());
gtag('config', 'G-VDC2S0ZNH5');
</script><script src="https://cdnjs.cloudflare.com/ajax/libs/autocomplete.js/0.38.0/autocomplete.jquery.min.js" integrity="sha512-GU9ayf+66Xx2TmpxqJpliWbT5PiGYxpaG8rfnBEk1LL8l1KGkRShhngwdXK1UgqhAzWpZHSiYPc09/NwDQIGyg==" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/mark.js/8.11.1/mark.min.js" integrity="sha512-5CYOlHXGh6QpOFA/TeTylKLWfB3ftPsde7AnmhuitiTX4K5SqCLBeKro6sPS8ilsz1Q4NRx3v8Ko2IBiszzdww==" crossorigin="anonymous"></script><!-- CSS --><link rel="stylesheet" href="style/style.css">
</head>
<body data-spy="scroll" data-target="#toc">
<div class="container-fluid">
<div class="row">
<header class="col-sm-12 col-lg-3 sidebar sidebar-book"><a class="sr-only sr-only-focusable" href="#content">Skip to main content</a>
<div class="d-flex align-items-start justify-content-between">
<h2>
<a href="index.html" title="">Geocomputation with R</a>
</h2>
<button class="btn btn-outline-primary d-lg-none ml-2 mt-1" type="button" data-toggle="collapse" data-target="#main-nav" aria-expanded="true" aria-controls="main-nav"><i class="fas fa-bars"></i><span class="sr-only">Show table of contents</span></button>
</div>
<div id="main-nav" class="collapse-lg">
<form role="search">
<input id="search" class="form-control" type="search" placeholder="Search" aria-label="Search">
</form>
<nav aria-label="Table of contents"><h2>Table of contents</h2>
<ul class="book-toc list-unstyled">
<li><a class="" href="index.html">Welcome</a></li>
<li><a class="" href="foreword-1st-edition.html">Foreword (1st Edition)</a></li>
<li><a class="" href="foreword-2nd-edition.html">Foreword (2nd Edition)</a></li>
<li><a class="" href="preface.html">Preface</a></li>
<li><a class="" href="intro.html"><span class="header-section-number">1</span> Introduction</a></li>
<li class="book-part">Foundations</li>
<li><a class="" href="spatial-class.html"><span class="header-section-number">2</span> Geographic data in R</a></li>
<li><a class="" href="attr.html"><span class="header-section-number">3</span> Attribute data operations</a></li>
<li><a class="" href="spatial-operations.html"><span class="header-section-number">4</span> Spatial data operations</a></li>
<li><a class="" href="geometry-operations.html"><span class="header-section-number">5</span> Geometry operations</a></li>
<li><a class="active" href="raster-vector.html"><span class="header-section-number">6</span> Raster-vector interactions</a></li>
<li><a class="" href="reproj-geo-data.html"><span class="header-section-number">7</span> Reprojecting geographic data</a></li>
<li><a class="" href="read-write.html"><span class="header-section-number">8</span> Geographic data I/O</a></li>
<li class="book-part">Extensions</li>
<li><a class="" href="adv-map.html"><span class="header-section-number">9</span> Making maps with R</a></li>
<li><a class="" href="gis.html"><span class="header-section-number">10</span> Bridges to GIS software</a></li>
<li><a class="" href="algorithms.html"><span class="header-section-number">11</span> Scripts, algorithms and functions</a></li>
<li><a class="" href="spatial-cv.html"><span class="header-section-number">12</span> Statistical learning</a></li>
<li class="book-part">Applications</li>
<li><a class="" href="transport.html"><span class="header-section-number">13</span> Transportation</a></li>
<li><a class="" href="location.html"><span class="header-section-number">14</span> Geomarketing</a></li>
<li><a class="" href="eco.html"><span class="header-section-number">15</span> Ecology</a></li>
<li><a class="" href="conclusion.html"><span class="header-section-number">16</span> Conclusion</a></li>
<li><a class="" href="references.html">References</a></li>
</ul>
<div class="book-extra">
<p><a id="book-repo" href="https://github.com/geocompx/geocompr">View book source <i class="fab fa-github"></i></a></p>
</div>
</nav>
</div>
</header><main class="col-sm-12 col-md-9 col-lg-7" id="content"><div id="raster-vector" class="section level1" number="6">
<h1>
<span class="header-section-number">6</span> Raster-vector interactions<a class="anchor" aria-label="anchor" href="#raster-vector"><i class="fas fa-link"></i></a>
</h1>
<div id="prerequisites-4" class="section level2 unnumbered">
<h2>Prerequisites<a class="anchor" aria-label="anchor" href="#prerequisites-4"><i class="fas fa-link"></i></a>
</h2>
<ul>
<li>This chapter requires the following packages:</li>
</ul>
<div class="sourceCode" id="cb199"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="kw"><a href="https://rdrr.io/r/base/library.html">library</a></span><span class="op">(</span><span class="va"><a href="https://r-spatial.github.io/sf/">sf</a></span><span class="op">)</span></span>
<span><span class="kw"><a href="https://rdrr.io/r/base/library.html">library</a></span><span class="op">(</span><span class="va"><a href="https://rspatial.org/">terra</a></span><span class="op">)</span></span>
<span><span class="kw"><a href="https://rdrr.io/r/base/library.html">library</a></span><span class="op">(</span><span class="va"><a href="https://dplyr.tidyverse.org">dplyr</a></span><span class="op">)</span></span></code></pre></div>
</div>
<div id="introduction-3" class="section level2" number="6.1">
<h2>
<span class="header-section-number">6.1</span> Introduction<a class="anchor" aria-label="anchor" href="#introduction-3"><i class="fas fa-link"></i></a>
</h2>
<p>
This Chapter focuses on interactions between raster and vector geographic data models, introduced in Chapter <a href="spatial-class.html#spatial-class">2</a>.
It includes several main techniques:
raster cropping and masking using vector objects (Section <a href="raster-vector.html#raster-cropping">6.2</a>);
extracting raster values using different types of vector data (Section <a href="raster-vector.html#raster-extraction">6.3</a>);
and raster-vector conversion (Sections <a href="raster-vector.html#rasterization">6.4</a> and <a href="raster-vector.html#spatial-vectorization">6.5</a>).
The above concepts are demonstrated using data used in previous chapters to understand their potential real-world applications.</p>
</div>
<div id="raster-cropping" class="section level2" number="6.2">
<h2>
<span class="header-section-number">6.2</span> Raster cropping<a class="anchor" aria-label="anchor" href="#raster-cropping"><i class="fas fa-link"></i></a>
</h2>
<p>
Many geographic data projects involve integrating data from many different sources, such as remote sensing images (rasters) and administrative boundaries (vectors).
Often the extent of input raster datasets is larger than the area of interest.
In this case, raster <strong>cropping</strong> and <strong>masking</strong> are useful for unifying the spatial extent of input data.
Both operations reduce object memory use and associated computational resources for subsequent analysis steps, and may be a necessary preprocessing step before creating attractive maps involving raster data.</p>
<p>We will use two objects to illustrate raster cropping:</p>
<ul>
<li>A <code>SpatRaster</code> object <code>srtm</code> representing elevation (meters above sea level) in south-western Utah</li>
<li>A vector (<code>sf</code>) object <code>zion</code> representing Zion National Park</li>
</ul>
<p>Both target and cropping objects must have the same projection.
The following code chunk therefore not only reads the datasets from the <strong>spDataLarge</strong> package installed in Chapter <a href="spatial-class.html#spatial-class">2</a>, it also ‘reprojects’ <code>zion</code> (a topic covered in Chapter <a href="reproj-geo-data.html#reproj-geo-data">7</a>):</p>
<div class="sourceCode" id="cb200"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">srtm</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/rast.html">rast</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/system.file.html">system.file</a></span><span class="op">(</span><span class="st">"raster/srtm.tif"</span>, package <span class="op">=</span> <span class="st">"spDataLarge"</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="va">zion</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_read.html">read_sf</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/system.file.html">system.file</a></span><span class="op">(</span><span class="st">"vector/zion.gpkg"</span>, package <span class="op">=</span> <span class="st">"spDataLarge"</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="va">zion</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_transform.html">st_transform</a></span><span class="op">(</span><span class="va">zion</span>, <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_crs.html">st_crs</a></span><span class="op">(</span><span class="va">srtm</span><span class="op">)</span><span class="op">)</span></span></code></pre></div>
<p>We use <code><a href="https://rspatial.github.io/terra/reference/crop.html">crop()</a></code> from the <strong>terra</strong> package to crop the <code>srtm</code> raster.
The function reduces the rectangular extent of the object passed to its first argument based on the extent of the object passed to its second argument.
This functionality is demonstrated in the command below, which generates Figure <a href="raster-vector.html#fig:cropmask">6.1</a>(B).</p>
<div class="sourceCode" id="cb201"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">srtm_cropped</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/crop.html">crop</a></span><span class="op">(</span><span class="va">srtm</span>, <span class="va">zion</span><span class="op">)</span></span></code></pre></div>
<p>
Related to <code><a href="https://rspatial.github.io/terra/reference/crop.html">crop()</a></code> is the <strong>terra</strong> function <code><a href="https://rspatial.github.io/terra/reference/mask.html">mask()</a></code>, which sets values outside of the bounds of the object passed to its second argument to <code>NA</code>.
The following command therefore masks every cell outside of the Zion National Park boundaries (Figure <a href="raster-vector.html#fig:cropmask">6.1</a>(C)).</p>
<div class="sourceCode" id="cb202"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">srtm_masked</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/mask.html">mask</a></span><span class="op">(</span><span class="va">srtm</span>, <span class="va">zion</span><span class="op">)</span></span></code></pre></div>
<p>Importantly, we want to use both <code><a href="https://rspatial.github.io/terra/reference/crop.html">crop()</a></code> and <code><a href="https://rspatial.github.io/terra/reference/mask.html">mask()</a></code> together in most cases.
This combination of functions would (a) limit the raster’s extent to our area of interest and then (b) replace all of the values outside of the area to NA.<a class="footnote-ref" tabindex="0" data-toggle="popover" data-content="<p>These two operations can be combined into a single step with <code>terra::crop(srtm, zion, mask = TRUE)</code>, but we prefer to keep them separate for clarity.</p>"><sup>28</sup></a></p>
<div class="sourceCode" id="cb203"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">srtm_cropped</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/crop.html">crop</a></span><span class="op">(</span><span class="va">srtm</span>, <span class="va">zion</span><span class="op">)</span></span>
<span><span class="va">srtm_final</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/mask.html">mask</a></span><span class="op">(</span><span class="va">srtm_cropped</span>, <span class="va">zion</span><span class="op">)</span></span></code></pre></div>
<p>Changing the settings of <code><a href="https://rspatial.github.io/terra/reference/mask.html">mask()</a></code> yields different results.
Setting <code>inverse = TRUE</code> will mask everything <em>inside</em> the bounds of the park (see <code><a href="https://rspatial.github.io/terra/reference/mask.html">?mask</a></code> for details) (Figure <a href="raster-vector.html#fig:cropmask">6.1</a>(D)), while setting <code>updatevalue = 0</code> will set all pixels outside the national park to 0.</p>
<div class="sourceCode" id="cb204"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">srtm_inv_masked</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/mask.html">mask</a></span><span class="op">(</span><span class="va">srtm</span>, <span class="va">zion</span>, inverse <span class="op">=</span> <span class="cn">TRUE</span><span class="op">)</span></span></code></pre></div>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:cropmask"></span>
<img src="figures/cropmask-1.png" alt="Illustration of raster cropping and raster masking." width="100%"><p class="caption">
FIGURE 6.1: Illustration of raster cropping and raster masking.
</p>
</div>
</div>
<div id="raster-extraction" class="section level2" number="6.3">
<h2>
<span class="header-section-number">6.3</span> Raster extraction<a class="anchor" aria-label="anchor" href="#raster-extraction"><i class="fas fa-link"></i></a>
</h2>
<p>
Raster extraction is the process of identifying and returning the values associated with a ‘target’ raster at specific locations, based on a (typically vector) geographic ‘selector’ object.
The results depend on the type of selector used (points, lines or polygons) and arguments passed to the <code><a href="https://rspatial.github.io/terra/reference/extract.html">terra::extract()</a></code> function.
The reverse of raster extraction — assigning raster cell values based on vector objects — is rasterization, described in Section <a href="raster-vector.html#rasterization">6.4</a>.</p>
<p>
The basic example is of extracting the value of a raster cell at specific <strong>points</strong>.
For this purpose, we will use <code>zion_points</code>, which contain a sample of 30 locations within the Zion National Park (Figure <a href="raster-vector.html#fig:pointextr">6.2</a>).
The following command extracts elevation values from <code>srtm</code> and creates a data frame with points’ IDs (one value per vector’s row) and related <code>srtm</code> values for each point.
Now, we can add the resulting object to our <code>zion_points</code> dataset with the <code><a href="https://rdrr.io/r/base/cbind.html">cbind()</a></code> function:</p>
<div class="sourceCode" id="cb205"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/data.html">data</a></span><span class="op">(</span><span class="st">"zion_points"</span>, package <span class="op">=</span> <span class="st">"spDataLarge"</span><span class="op">)</span></span>
<span><span class="va">elevation</span> <span class="op">=</span> <span class="fu">terra</span><span class="fu">::</span><span class="fu"><a href="https://rspatial.github.io/terra/reference/extract.html">extract</a></span><span class="op">(</span><span class="va">srtm</span>, <span class="va">zion_points</span><span class="op">)</span></span>
<span><span class="va">zion_points</span> <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/cbind.html">cbind</a></span><span class="op">(</span><span class="va">zion_points</span>, <span class="va">elevation</span><span class="op">)</span></span></code></pre></div>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:pointextr"></span>
<img src="figures/pointextr-1.png" alt="Locations of points used for raster extraction." width="100%"><p class="caption">
FIGURE 6.2: Locations of points used for raster extraction.
</p>
</div>
<p>
Raster extraction also works with <strong>line</strong> selectors.
Then, it extracts one value for each raster cell touched by a line.
However, the line extraction approach is not recommended to obtain values along the transects as it is hard to get the correct distance between each pair of extracted raster values.</p>
<p>In this case, a better approach is to split the line into many points and then extract the values for these points.
To demonstrate this, the code below creates <code>zion_transect</code>, a straight line going from northwest to southeast of the Zion National Park, illustrated in Figure <a href="raster-vector.html#fig:lineextr">6.3</a>(A) (see Section <a href="spatial-class.html#vector-data">2.2</a> for a recap on the vector data model):</p>
<div class="sourceCode" id="cb206"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">zion_transect</span> <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/cbind.html">cbind</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span><span class="op">(</span><span class="op">-</span><span class="fl">113.2</span>, <span class="op">-</span><span class="fl">112.9</span><span class="op">)</span>, <span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span><span class="op">(</span><span class="fl">37.45</span>, <span class="fl">37.2</span><span class="op">)</span><span class="op">)</span> <span class="op">|></span></span>
<span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st.html">st_linestring</a></span><span class="op">(</span><span class="op">)</span> <span class="op">|></span> </span>
<span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/sfc.html">st_sfc</a></span><span class="op">(</span>crs <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/crs.html">crs</a></span><span class="op">(</span><span class="va">srtm</span><span class="op">)</span><span class="op">)</span> <span class="op">|></span></span>
<span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/sf.html">st_sf</a></span><span class="op">(</span>geometry <span class="op">=</span> <span class="va">_</span><span class="op">)</span></span></code></pre></div>
<p>The utility of extracting heights from a linear selector is illustrated by imagining that you are planning a hike.
The method demonstrated below provides an ‘elevation profile’ of the route (the line does not need to be straight), useful for estimating how long it will take due to long climbs.</p>
<p>The first step is to add a unique <code>id</code> for each transect.
Next, with the <code><a href="https://r-spatial.github.io/sf/reference/geos_unary.html">st_segmentize()</a></code> function we can add points along our line(s) with a provided density (<code>dfMaxLength</code>) and convert them into points with <code><a href="https://r-spatial.github.io/sf/reference/st_cast.html">st_cast()</a></code>.</p>
<div class="sourceCode" id="cb207"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">zion_transect</span><span class="op">$</span><span class="va">id</span> <span class="op">=</span> <span class="fl">1</span><span class="op">:</span><span class="fu"><a href="https://rspatial.github.io/terra/reference/dimensions.html">nrow</a></span><span class="op">(</span><span class="va">zion_transect</span><span class="op">)</span></span>
<span><span class="va">zion_transect</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/geos_unary.html">st_segmentize</a></span><span class="op">(</span><span class="va">zion_transect</span>, dfMaxLength <span class="op">=</span> <span class="fl">250</span><span class="op">)</span></span>
<span><span class="va">zion_transect</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_cast.html">st_cast</a></span><span class="op">(</span><span class="va">zion_transect</span>, <span class="st">"POINT"</span><span class="op">)</span></span></code></pre></div>
<p>Now, we have a large set of points, and we want to derive a distance between the first point in our transects and each of the subsequent points.
In this case, we only have one transect, but the code, in principle, should work on any number of transects:</p>
<div class="sourceCode" id="cb208"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">zion_transect</span> <span class="op">=</span> <span class="va">zion_transect</span> <span class="op">|></span> </span>
<span> <span class="fu"><a href="https://dplyr.tidyverse.org/reference/group_by.html">group_by</a></span><span class="op">(</span><span class="va">id</span><span class="op">)</span> <span class="op">|></span> </span>
<span> <span class="fu"><a href="https://dplyr.tidyverse.org/reference/mutate.html">mutate</a></span><span class="op">(</span>dist <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/geos_measures.html">st_distance</a></span><span class="op">(</span><span class="va">geometry</span><span class="op">)</span><span class="op">[</span>, <span class="fl">1</span><span class="op">]</span><span class="op">)</span> </span></code></pre></div>
<p>Finally, we can extract elevation values for each point in our transects and combine this information with our main object.</p>
<div class="sourceCode" id="cb209"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">zion_elev</span> <span class="op">=</span> <span class="fu">terra</span><span class="fu">::</span><span class="fu"><a href="https://rspatial.github.io/terra/reference/extract.html">extract</a></span><span class="op">(</span><span class="va">srtm</span>, <span class="va">zion_transect</span><span class="op">)</span></span>
<span><span class="va">zion_transect</span> <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/cbind.html">cbind</a></span><span class="op">(</span><span class="va">zion_transect</span>, <span class="va">zion_elev</span><span class="op">)</span></span></code></pre></div>
<p>The resulting <code>zion_transect</code> can be used to create elevation profiles, as illustrated in Figure <a href="raster-vector.html#fig:lineextr">6.3</a>(B).</p>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:lineextr"></span>
<img src="figures/lineextr-1.png" alt="Location of a line used for raster extraction (left) and the elevation along this line (right)." width="100%"><p class="caption">
FIGURE 6.3: Location of a line used for raster extraction (left) and the elevation along this line (right).
</p>
</div>
<p>
The final type of geographic vector object for raster extraction is <strong>polygons</strong>.
Like lines, polygons tend to return many raster values per polygon.
This is demonstrated in the command below, which results in a data frame with column names <code>ID</code> (the row number of the polygon) and <code>srtm</code> (associated elevation values):</p>
<div class="sourceCode" id="cb210"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">zion_srtm_values</span> <span class="op">=</span> <span class="fu">terra</span><span class="fu">::</span><span class="fu"><a href="https://rspatial.github.io/terra/reference/extract.html">extract</a></span><span class="op">(</span>x <span class="op">=</span> <span class="va">srtm</span>, y <span class="op">=</span> <span class="va">zion</span><span class="op">)</span></span></code></pre></div>
<p>Such results can be used to generate summary statistics for raster values per polygon, for example to characterize a single region or to compare many regions.
This is shown in the code below, which creates the object <code>zion_srtm_df</code> containing summary statistics for elevation values in Zion National Park (see Figure <a href="raster-vector.html#fig:polyextr">6.4</a>(A)):</p>
<div class="sourceCode" id="cb211"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://dplyr.tidyverse.org/reference/group_by.html">group_by</a></span><span class="op">(</span><span class="va">zion_srtm_values</span>, <span class="va">ID</span><span class="op">)</span> <span class="op">|></span> </span>
<span> <span class="fu"><a href="https://dplyr.tidyverse.org/reference/summarise.html">summarize</a></span><span class="op">(</span><span class="fu"><a href="https://dplyr.tidyverse.org/reference/across.html">across</a></span><span class="op">(</span><span class="va">srtm</span>, <span class="fu"><a href="https://rdrr.io/r/base/list.html">list</a></span><span class="op">(</span>min <span class="op">=</span> <span class="va">min</span>, mean <span class="op">=</span> <span class="va">mean</span>, max <span class="op">=</span> <span class="va">max</span><span class="op">)</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="co">#> # A tibble: 1 × 4</span></span>
<span><span class="co">#> ID srtm_min srtm_mean srtm_max</span></span>
<span><span class="co">#> <dbl> <int> <dbl> <int></span></span>
<span><span class="co">#> 1 1 1122 1818. 2661</span></span></code></pre></div>
<p>The preceding code chunk used <strong>dplyr</strong> to provide summary statistics for cell values per polygon ID, as described in Chapter <a href="attr.html#attr">3</a>.
The results provide useful summaries, for example that the maximum height in the park is around 2,661 meters above sea level (other summary statistics, such as standard deviation, can also be calculated in this way).
Because there is only one polygon in the example, a data frame with a single row is returned; however, the method works when multiple selector polygons are used.</p>
<p>A similar approach works for counting occurrences of categorical raster values within polygons.
This is illustrated with a land cover dataset (<code>nlcd</code>) from the <strong>spDataLarge</strong> package in Figure <a href="raster-vector.html#fig:polyextr">6.4</a>(B), and demonstrated in the code below:</p>
<div class="sourceCode" id="cb212"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">nlcd</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/rast.html">rast</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/system.file.html">system.file</a></span><span class="op">(</span><span class="st">"raster/nlcd.tif"</span>, package <span class="op">=</span> <span class="st">"spDataLarge"</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="va">zion2</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_transform.html">st_transform</a></span><span class="op">(</span><span class="va">zion</span>, <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_crs.html">st_crs</a></span><span class="op">(</span><span class="va">nlcd</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="va">zion_nlcd</span> <span class="op">=</span> <span class="fu">terra</span><span class="fu">::</span><span class="fu"><a href="https://rspatial.github.io/terra/reference/extract.html">extract</a></span><span class="op">(</span><span class="va">nlcd</span>, <span class="va">zion2</span><span class="op">)</span></span>
<span><span class="va">zion_nlcd</span> <span class="op">|></span> </span>
<span> <span class="fu"><a href="https://dplyr.tidyverse.org/reference/group_by.html">group_by</a></span><span class="op">(</span><span class="va">ID</span>, <span class="va">levels</span><span class="op">)</span> <span class="op">|></span></span>
<span> <span class="fu"><a href="https://dplyr.tidyverse.org/reference/count.html">count</a></span><span class="op">(</span><span class="op">)</span></span>
<span><span class="co">#> # A tibble: 7 × 3</span></span>
<span><span class="co">#> # Groups: ID, levels [7]</span></span>
<span><span class="co">#> ID levels n</span></span>
<span><span class="co">#> <dbl> <fct> <int></span></span>
<span><span class="co">#> 1 1 Developed 4205</span></span>
<span><span class="co">#> 2 1 Barren 98285</span></span>
<span><span class="co">#> 3 1 Forest 298299</span></span>
<span><span class="co">#> 4 1 Shrubland 203700</span></span>
<span><span class="co">#> # ℹ 3 more rows</span></span></code></pre></div>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:polyextr"></span>
<img src="figures/polyextr-1.png" alt="Area used for continuous (left) and categorical (right) raster extraction." width="100%"><p class="caption">
FIGURE 6.4: Area used for continuous (left) and categorical (right) raster extraction.
</p>
</div>
<p>
Although the <strong>terra</strong> package offers rapid extraction of raster values within polygons, <code><a href="https://rspatial.github.io/terra/reference/extract.html">extract()</a></code> can still be a bottleneck when processing large polygon datasets.
The <strong>exactextractr</strong> package offers a <a href="https://github.com/geocompx/geocompr/issues/813">significantly faster alternative</a> for extracting pixel values through the <code>exact_extract()</code> function.
The <code>exact_extract()</code> function also computes, by default, the fraction of each raster cell overlapped by the polygon, which is more precise (see note below for details).</p>
<div class="rmdnote">
Polygons usually have irregular shapes, and therefore, a polygon can overlap only some parts of a raster’s cells.
To get more detailed results, the <code><a href="https://rspatial.github.io/terra/reference/extract.html">terra::extract()</a></code> function has an argument called <code>exact</code>.
With <code>exact = TRUE</code>, we get one more column <code>fraction</code> in the output data frame, which represents a fraction of each cell that is covered by the polygon.
This could be useful to calculate, for example, a weighted mean for continuous rasters or more precise coverage for categorical rasters.
By default, it is <code>FALSE</code> as this operation requires more computations.
The <code><a href="https://isciences.gitlab.io/exactextractr/reference/exact_extract.html">exactextractr::exact_extract()</a></code> function always computes the coverage fraction of the polygon in each cell.
</div>
</div>
<div id="rasterization" class="section level2" number="6.4">
<h2>
<span class="header-section-number">6.4</span> Rasterization<a class="anchor" aria-label="anchor" href="#rasterization"><i class="fas fa-link"></i></a>
</h2>
<p>
Rasterization is the conversion of vector objects into their representation in raster objects.
Usually, the output raster is then used for quantitative analysis (e.g., analysis of terrain) or modeling.
As we saw in Chapter <a href="spatial-class.html#spatial-class">2</a>, the raster data model has some characteristics that make it conducive to certain methods.
Furthermore, the process of rasterization can help simplify datasets because the resulting values all have the same spatial resolution: rasterization can be seen as a special type of geographic data aggregation.</p>
<p>The <strong>terra</strong> package contains the function <code><a href="https://rspatial.github.io/terra/reference/rasterize.html">rasterize()</a></code> for doing this work.
Its first two arguments are, <code>x</code>, vector object to be rasterized and, <code>y</code>, a ‘template raster’ object defining the extent, resolution and CRS of the output.
The geographic resolution of the input raster has a major impact on the results: if it is too low (cell size is too large), the result may miss the full geographic variability of the vector data; if it is too high, computational times may be excessive.
There are no simple rules to follow when deciding an appropriate geographic resolution, which is heavily dependent on the intended use of the results.
Often the target resolution is imposed on the user, for example when the output of rasterization needs to be aligned to some other existing raster.</p>
<p>
To demonstrate rasterization in action, we will use a template raster that has the same extent and CRS as the input vector data <code>cycle_hire_osm_projected</code> (a dataset on cycle hire points in London is illustrated in Figure <a href="raster-vector.html#fig:vector-rasterization1">6.5</a>(A)) and spatial resolution of 1000 meters:</p>
<div class="sourceCode" id="cb213"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">cycle_hire_osm</span> <span class="op">=</span> <span class="fu">spData</span><span class="fu">::</span><span class="va"><a href="https://jakubnowosad.com/spData/reference/cycle_hire_osm.html">cycle_hire_osm</a></span></span>
<span><span class="va">cycle_hire_osm_projected</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_transform.html">st_transform</a></span><span class="op">(</span><span class="va">cycle_hire_osm</span>, <span class="st">"EPSG:27700"</span><span class="op">)</span></span>
<span><span class="va">raster_template</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/rast.html">rast</a></span><span class="op">(</span><span class="fu"><a href="https://rspatial.github.io/terra/reference/ext.html">ext</a></span><span class="op">(</span><span class="va">cycle_hire_osm_projected</span><span class="op">)</span>, resolution <span class="op">=</span> <span class="fl">1000</span>,</span>
<span> crs <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/crs.html">crs</a></span><span class="op">(</span><span class="va">cycle_hire_osm_projected</span><span class="op">)</span><span class="op">)</span></span></code></pre></div>
<p>Rasterization is a very flexible operation: the results depend not only on the nature of the template raster, but also on the type of input vector (e.g., points, polygons) and a variety of arguments taken by the <code><a href="https://rspatial.github.io/terra/reference/rasterize.html">rasterize()</a></code> function.</p>
<p>To illustrate this flexibility we will try three different approaches to rasterization.
First, we create a raster representing the presence or absence of cycle hire points (known as presence/absence rasters).
In this case <code><a href="https://rspatial.github.io/terra/reference/rasterize.html">rasterize()</a></code> requires no argument in addition to <code>x</code> and <code>y</code>, the aforementioned vector and raster objects (results illustrated Figure <a href="raster-vector.html#fig:vector-rasterization1">6.5</a>(B)).</p>
<div class="sourceCode" id="cb214"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">ch_raster1</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/rasterize.html">rasterize</a></span><span class="op">(</span><span class="va">cycle_hire_osm_projected</span>, <span class="va">raster_template</span><span class="op">)</span></span></code></pre></div>
<p>The <code>fun</code> argument specifies summary statistics used to convert multiple observations in close proximity into associate cells in the raster object.
By default <code>fun = "last"</code> is used but other options such as <code>fun = "length"</code> can be used, in this case to count the number of cycle hire points in each grid cell (the results of this operation are illustrated in Figure <a href="raster-vector.html#fig:vector-rasterization1">6.5</a>(C)).</p>
<div class="sourceCode" id="cb215"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">ch_raster2</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/rasterize.html">rasterize</a></span><span class="op">(</span><span class="va">cycle_hire_osm_projected</span>, <span class="va">raster_template</span>, </span>
<span> fun <span class="op">=</span> <span class="st">"length"</span><span class="op">)</span></span></code></pre></div>
<p>The new output, <code>ch_raster2</code>, shows the number of cycle hire points in each grid cell.
The cycle hire locations have different numbers of bicycles described by the <code>capacity</code> variable, raising the question, what’s the capacity in each grid cell?
To calculate that we must <code>sum</code> the field (<code>"capacity"</code>), resulting in output illustrated in Figure <a href="raster-vector.html#fig:vector-rasterization1">6.5</a>(D), calculated with the following command (other summary functions such as <code>mean</code> could be used).</p>
<div class="sourceCode" id="cb216"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">ch_raster3</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/rasterize.html">rasterize</a></span><span class="op">(</span><span class="va">cycle_hire_osm_projected</span>, <span class="va">raster_template</span>, </span>
<span> field <span class="op">=</span> <span class="st">"capacity"</span>, fun <span class="op">=</span> <span class="va">sum</span>, na.rm <span class="op">=</span> <span class="cn">TRUE</span><span class="op">)</span></span></code></pre></div>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:vector-rasterization1"></span>
<img src="figures/vector-rasterization1-1.png" alt="Examples of point rasterization." width="100%"><p class="caption">
FIGURE 6.5: Examples of point rasterization.
</p>
</div>
<p>
Another dataset based on California’s polygons and borders (created below) illustrates rasterization of lines.
After casting the polygon objects into a multilinestring, a template raster is created with a resolution of a 0.5 degree:</p>
<div class="sourceCode" id="cb217"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">california</span> <span class="op">=</span> <span class="fu">dplyr</span><span class="fu">::</span><span class="fu"><a href="https://dplyr.tidyverse.org/reference/filter.html">filter</a></span><span class="op">(</span><span class="va">us_states</span>, <span class="va">NAME</span> <span class="op">==</span> <span class="st">"California"</span><span class="op">)</span></span>
<span><span class="va">california_borders</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_cast.html">st_cast</a></span><span class="op">(</span><span class="va">california</span>, <span class="st">"MULTILINESTRING"</span><span class="op">)</span></span>
<span><span class="va">raster_template2</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/rast.html">rast</a></span><span class="op">(</span><span class="fu"><a href="https://rspatial.github.io/terra/reference/ext.html">ext</a></span><span class="op">(</span><span class="va">california</span><span class="op">)</span>, resolution <span class="op">=</span> <span class="fl">0.5</span>,</span>
<span> crs <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_crs.html">st_crs</a></span><span class="op">(</span><span class="va">california</span><span class="op">)</span><span class="op">$</span><span class="va">wkt</span><span class="op">)</span></span></code></pre></div>
<p>When considering line or polygon rasterization, one useful additional argument is <code>touches</code>.
By default it is <code>FALSE</code>, but when changed to <code>TRUE</code> – all cells that are touched by a line or polygon border get a value.
Line rasterization with <code>touches = TRUE</code> is demonstrated in the code below (Figure <a href="raster-vector.html#fig:vector-rasterization2">6.6</a>(A)).</p>
<div class="sourceCode" id="cb218"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">california_raster1</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/rasterize.html">rasterize</a></span><span class="op">(</span><span class="va">california_borders</span>, <span class="va">raster_template2</span>,</span>
<span> touches <span class="op">=</span> <span class="cn">TRUE</span><span class="op">)</span></span></code></pre></div>
<p>Compare it to a polygon rasterization, with <code>touches = FALSE</code> by default, which selects only raster cells whose centroids are inside the selector polygon, as illustrated in Figure <a href="raster-vector.html#fig:vector-rasterization2">6.6</a>(B).</p>
<div class="sourceCode" id="cb219"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">california_raster2</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/rasterize.html">rasterize</a></span><span class="op">(</span><span class="va">california</span>, <span class="va">raster_template2</span><span class="op">)</span> </span></code></pre></div>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:vector-rasterization2"></span>
<img src="figures/vector-rasterization2-1.png" alt="Examples of line and polygon rasterizations." width="100%"><p class="caption">
FIGURE 6.6: Examples of line and polygon rasterizations.
</p>
</div>
</div>
<div id="spatial-vectorization" class="section level2" number="6.5">
<h2>
<span class="header-section-number">6.5</span> Spatial vectorization<a class="anchor" aria-label="anchor" href="#spatial-vectorization"><i class="fas fa-link"></i></a>
</h2>
<p>
Spatial vectorization is the counterpart of rasterization (Section <a href="raster-vector.html#rasterization">6.4</a>), but in the opposite direction.
It involves converting spatially continuous raster data into spatially discrete vector data such as points, lines or polygons.</p>
<div class="rmdnote">
Be careful with the wording!
In R, vectorization usually refers to the possibility of replacing <code>for</code>-loops and alike by doing things like <code>1:10 / 2</code> (see also <span class="citation">Wickham (<a href="references.html#ref-wickham_advanced_2019">2019</a>)</span>).
</div>
<p>
The simplest form of vectorization is to convert the centroids of raster cells into points.
<code><a href="https://rspatial.github.io/terra/reference/as.points.html">as.points()</a></code> does exactly this for all non-<code>NA</code> raster grid cells (Figure <a href="raster-vector.html#fig:raster-vectorization1">6.7</a>).
Note, here we also used <code><a href="https://r-spatial.github.io/sf/reference/st_as_sf.html">st_as_sf()</a></code> to convert the resulting object to the <code>sf</code> class.</p>
<div class="sourceCode" id="cb220"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">elev</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/rast.html">rast</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/system.file.html">system.file</a></span><span class="op">(</span><span class="st">"raster/elev.tif"</span>, package <span class="op">=</span> <span class="st">"spData"</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="va">elev_point</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/as.points.html">as.points</a></span><span class="op">(</span><span class="va">elev</span><span class="op">)</span> <span class="op">|></span> </span>
<span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_as_sf.html">st_as_sf</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:raster-vectorization1"></span>
<img src="figures/raster-vectorization1-1.png" alt="Raster and point representation of the elev object." width="100%"><p class="caption">
FIGURE 6.7: Raster and point representation of the elev object.
</p>
</div>
<p>
Another common type of spatial vectorization is the creation of contour lines representing lines of continuous height or temperatures (isotherms) for example.
We will use a real-world digital elevation model (DEM) because the artificial raster <code>elev</code> produces parallel lines (task for the reader: verify this and explain why this happens).
Contour lines can be created with the <strong>terra</strong> function <code><a href="https://rspatial.github.io/terra/reference/contour.html">as.contour()</a></code>, which is itself a wrapper around the built-in R function <code><a href="https://rdrr.io/r/graphics/filled.contour.html">filled.contour()</a></code>, as demonstrated below (not shown):</p>
<div class="sourceCode" id="cb221"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">dem</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/rast.html">rast</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/system.file.html">system.file</a></span><span class="op">(</span><span class="st">"raster/dem.tif"</span>, package <span class="op">=</span> <span class="st">"spDataLarge"</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="va">cl</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/contour.html">as.contour</a></span><span class="op">(</span><span class="va">dem</span><span class="op">)</span> <span class="op">|></span> </span>
<span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_as_sf.html">st_as_sf</a></span><span class="op">(</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rspatial.github.io/terra/reference/plot.html">plot</a></span><span class="op">(</span><span class="va">dem</span>, axes <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rspatial.github.io/terra/reference/plot.html">plot</a></span><span class="op">(</span><span class="va">cl</span>, add <span class="op">=</span> <span class="cn">TRUE</span><span class="op">)</span></span></code></pre></div>
<p>Contours can also be added to existing plots with functions such as <code><a href="https://rspatial.github.io/terra/reference/contour.html">contour()</a></code>, <code>rasterVis::contourplot()</code>.
<!-- or `tmap::tm_iso()` (not yet implemented as of 2023-11-24) -->
As illustrated in Figure <a href="raster-vector.html#fig:contour-tmap">6.8</a>, isolines can be labelled.</p>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:contour-tmap"></span>
<img src="figures/06-contour-tmap.png" alt="DEM with hillshading, showing the southern flank of Mt. Mongón overlaid with contour lines." width="100%"><p class="caption">
FIGURE 6.8: DEM with hillshading, showing the southern flank of Mt. Mongón overlaid with contour lines.
</p>
</div>
<p>
The final type of vectorization involves conversion of rasters to polygons.
This can be done with <code><a href="https://rspatial.github.io/terra/reference/as.polygons.html">terra::as.polygons()</a></code>, which converts each raster cell into a polygon consisting of five coordinates, all of which are stored in memory (explaining why rasters are often fast compared with vectors!).</p>
<p>This is illustrated below by converting the <code>grain</code> object into polygons and subsequently dissolving borders between polygons with the same attribute values (also see the <code>dissolve</code> argument in <code><a href="https://rspatial.github.io/terra/reference/as.polygons.html">as.polygons()</a></code>).</p>
<div class="sourceCode" id="cb222"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">grain</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/rast.html">rast</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/system.file.html">system.file</a></span><span class="op">(</span><span class="st">"raster/grain.tif"</span>, package <span class="op">=</span> <span class="st">"spData"</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="va">grain_poly</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/as.polygons.html">as.polygons</a></span><span class="op">(</span><span class="va">grain</span><span class="op">)</span> <span class="op">|></span> </span>
<span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_as_sf.html">st_as_sf</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<pre><code>#> Some legend items or map compoments do not fit well (e.g. due to the specified font size). They are rescaled Set the tmap option component.autoscale to FALSE to disable this.FALSE</code></pre>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:06-raster-vector-40"></span>
<img src="figures/06-raster-vector-40-1.png" alt="Illustration of vectorization of raster (left) into polygons (dissolve = FALSE; center) and aggregated polygons (dissolve = TRUE; right)." width="100%"><p class="caption">
FIGURE 6.9: Illustration of vectorization of raster (left) into polygons (dissolve = FALSE; center) and aggregated polygons (dissolve = TRUE; right).
</p>
</div>
<p>The aggregated polygons of the <code>grain</code> dataset have rectilinear boundaries which arise from being defined by connecting rectangular pixels.
The <strong>smoothr</strong> package described in Chapter <a href="geometry-operations.html#geometry-operations">5</a> can be used to smooth the edges of the polygons.
As smoothing removes sharp edges in the polygon boundaries, the smoothed polygons will not have the same exact spatial coverage as the original pixels.
Caution should therefore be taken when using the smoothed polygons for further analysis.</p>
</div>
<div id="exercises-4" class="section level2" number="6.6">
<h2>
<span class="header-section-number">6.6</span> Exercises<a class="anchor" aria-label="anchor" href="#exercises-4"><i class="fas fa-link"></i></a>
</h2>
<p>Some of the following exercises use a vector (<code>zion_points</code>) and raster dataset (<code>srtm</code>) from the <strong>spDataLarge</strong> package.
They also use a polygonal ‘convex hull’ derived from the vector dataset (<code>ch</code>) to represent the area of interest:</p>
<div class="sourceCode" id="cb224"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="kw"><a href="https://rdrr.io/r/base/library.html">library</a></span><span class="op">(</span><span class="va"><a href="https://r-spatial.github.io/sf/">sf</a></span><span class="op">)</span></span>
<span><span class="kw"><a href="https://rdrr.io/r/base/library.html">library</a></span><span class="op">(</span><span class="va"><a href="https://rspatial.org/">terra</a></span><span class="op">)</span></span>
<span><span class="kw"><a href="https://rdrr.io/r/base/library.html">library</a></span><span class="op">(</span><span class="va"><a href="https://jakubnowosad.com/spData/">spData</a></span><span class="op">)</span></span>
<span><span class="va">zion_points_path</span> <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/system.file.html">system.file</a></span><span class="op">(</span><span class="st">"vector/zion_points.gpkg"</span>, package <span class="op">=</span> <span class="st">"spDataLarge"</span><span class="op">)</span></span>
<span><span class="va">zion_points</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_read.html">read_sf</a></span><span class="op">(</span><span class="va">zion_points_path</span><span class="op">)</span></span>
<span><span class="va">srtm</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/rast.html">rast</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/system.file.html">system.file</a></span><span class="op">(</span><span class="st">"raster/srtm.tif"</span>, package <span class="op">=</span> <span class="st">"spDataLarge"</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="va">ch</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/geos_combine.html">st_combine</a></span><span class="op">(</span><span class="va">zion_points</span><span class="op">)</span> <span class="op">|></span></span>
<span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/geos_unary.html">st_convex_hull</a></span><span class="op">(</span><span class="op">)</span> <span class="op">|></span> </span>
<span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_as_sf.html">st_as_sf</a></span><span class="op">(</span><span class="op">)</span></span></code></pre></div>
<p>E1. Crop the <code>srtm</code> raster using (1) the <code>zion_points</code> dataset and (2) the <code>ch</code> dataset.
Are there any differences in the output maps?
Next, mask <code>srtm</code> using these two datasets.
Can you see any difference now?
How can you explain that?</p>
<p>E2. Firstly, extract values from <code>srtm</code> at the points represented in <code>zion_points</code>.
Next, extract average values of <code>srtm</code> using a 90 buffer around each point from <code>zion_points</code> and compare these two sets of values.
When would extracting values by buffers be more suitable than by points alone?</p>
<ul>
<li>Bonus: Implement extraction using the <strong>exactextractr</strong> package and compare the results.</li>
</ul>
<p>E3. Subset points higher than 3100 meters in New Zealand (the <code>nz_height</code> object) and create a template raster with a resolution of 3 km for the extent of the new point dataset.
Using these two new objects:</p>
<ul>
<li>Count numbers of the highest points in each grid cell.</li>
<li>Find the maximum elevation in each grid cell.</li>
</ul>
<p>E4. Aggregate the raster counting high points in New Zealand (created in the previous exercise), reduce its geographic resolution by half (so cells are 6 by 6 km) and plot the result.</p>
<ul>
<li>Resample the lower resolution raster back to the original resolution of 3 km. How have the results changed?</li>
<li>Name two advantages and disadvantages of reducing raster resolution.</li>
</ul>
<p>E5. Polygonize the <code>grain</code> dataset and filter all squares representing clay.</p>
<ul>
<li>Name two advantages and disadvantages of vector data over raster data.</li>
<li>When would it be useful to convert rasters to vectors in your work?</li>
</ul>
</div>
</div>
<div class="chapter-nav">
<div class="prev"><a href="geometry-operations.html"><span class="header-section-number">5</span> Geometry operations</a></div>
<div class="next"><a href="reproj-geo-data.html"><span class="header-section-number">7</span> Reprojecting geographic data</a></div>
</div></main><div class="col-md-3 col-lg-2 d-none d-md-block sidebar sidebar-chapter">
<h2>Note: Second Edition is under final polishing 🏗</h2>
<!--<p>Now is a great time to provide feedback</p>-->
<ul class="list-unstyled">
<!--<li><a href="https://forms.gle/nq9RmbxJyZXQgc948">Provide feedback (5 min)</a></li>--><li><a href="https://geocompx.org/">geocompx 🌐</a></li>
<li><a href="https://r.geocompx.org/#reproducibility">Install updated packages</a></li>
<li><a href="https://github.com/geocompx/geocompr/issues">Open an issue <i class="fas fa-question"></i></a></li>
<li><a href="https://discord.gg/PMztXYgNxp">Chat on Discord <i class="fab fa-discord"></i></a></li>
<li><a href="https://r.geocompx.org/solutions/">Check exercise solutions <i class="fa fa-check"></i></a></li>
<li><a href="https://supportukrainenow.org/">Support Ukraine 🇺🇦
</a></li>
</ul>
<hr>
<nav id="toc" data-toggle="toc" aria-label="On this page"><h2>On this page</h2>
<ul class="nav navbar-nav">
<li><a class="nav-link" href="#raster-vector"><span class="header-section-number">6</span> Raster-vector interactions</a></li>
<li><a class="nav-link" href="#prerequisites-4">Prerequisites</a></li>
<li><a class="nav-link" href="#introduction-3"><span class="header-section-number">6.1</span> Introduction</a></li>
<li><a class="nav-link" href="#raster-cropping"><span class="header-section-number">6.2</span> Raster cropping</a></li>
<li><a class="nav-link" href="#raster-extraction"><span class="header-section-number">6.3</span> Raster extraction</a></li>
<li><a class="nav-link" href="#rasterization"><span class="header-section-number">6.4</span> Rasterization</a></li>
<li><a class="nav-link" href="#spatial-vectorization"><span class="header-section-number">6.5</span> Spatial vectorization</a></li>
<li><a class="nav-link" href="#exercises-4"><span class="header-section-number">6.6</span> Exercises</a></li>
</ul>
<div class="book-extra">
<ul class="list-unstyled">
<li><a id="book-source" href="https://github.com/geocompx/geocompr/blob/main/06-raster-vector.Rmd">View source <i class="fab fa-github"></i></a></li>
<li><a id="book-edit" href="https://github.com/geocompx/geocompr/edit/main/06-raster-vector.Rmd">Edit this page <i class="fab fa-github"></i></a></li>
</ul>
</div>
</nav>
</div>
</div>
</div> <!-- .container -->
<footer class="bg-primary text-light mt-5"><div class="container"><div class="row">
<div class="col-12 col-md-6 mt-3">
<p>"<strong>Geocomputation with R</strong>" was written by Robin Lovelace, Jakub Nowosad, Jannes Muenchow. It was last built on 2024-04-26.</p>
</div>
<div class="col-12 col-md-6 mt-3">
<p>This book was built by the <a class="text-light" href="https://bookdown.org">bookdown</a> R package.</p>
</div>
</div></div>
</footer><!-- dynamically load mathjax for compatibility with self-contained --><script>
(function () {
var script = document.createElement("script");
script.type = "text/javascript";
var src = "true";
if (src === "" || src === "true") src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-MML-AM_CHTML";
if (location.protocol !== "file:")
if (/^https?:/.test(src))
src = src.replace(/^https?:/, '');
script.src = src;
document.getElementsByTagName("head")[0].appendChild(script);
})();
</script><script type="text/x-mathjax-config">const popovers = document.querySelectorAll('a.footnote-ref[data-toggle="popover"]');
for (let popover of popovers) {
const div = document.createElement('div');
div.setAttribute('style', 'position: absolute; top: 0, left:0; width:0, height:0, overflow: hidden; visibility: hidden;');
div.innerHTML = popover.getAttribute('data-content');
var has_math = div.querySelector("span.math");
if (has_math) {
document.body.appendChild(div);
MathJax.Hub.Queue(["Typeset", MathJax.Hub, div]);
MathJax.Hub.Queue(function() {
popover.setAttribute('data-content', div.innerHTML);
document.body.removeChild(div);
})
}
}
</script>
</body>
</html>