/
spatial-operations.html
1054 lines (1033 loc) · 179 KB
/
spatial-operations.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
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!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 4 Spatial data operations | Geocomputation with R</title>
<meta name="author" content="Robin Lovelace, Jakub Nowosad, Jannes Muenchow">
<meta name="description" content="Prerequisites This chapter requires the same packages used in Chapter 3: library(sf) library(terra) library(dplyr) library(spData) 4.1 Introduction Spatial operations, including spatial joins...">
<meta name="generator" content="bookdown 0.39 with bs4_book()">
<meta property="og:title" content="Chapter 4 Spatial data operations | Geocomputation with R">
<meta property="og:type" content="book">
<meta property="og:url" content="https://r.geocompx.org/spatial-operations.html">
<meta property="og:image" content="https://r.geocompx.org/images/cover.png">
<meta property="og:description" content="Prerequisites This chapter requires the same packages used in Chapter 3: library(sf) library(terra) library(dplyr) library(spData) 4.1 Introduction Spatial operations, including spatial joins...">
<meta name="twitter:card" content="summary">
<meta name="twitter:title" content="Chapter 4 Spatial data operations | Geocomputation with R">
<meta name="twitter:description" content="Prerequisites This chapter requires the same packages used in Chapter 3: library(sf) library(terra) library(dplyr) library(spData) 4.1 Introduction Spatial operations, including spatial joins...">
<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 4 Spatial data operations | 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="active" 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="" 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="spatial-operations" class="section level1" number="4">
<h1>
<span class="header-section-number">4</span> Spatial data operations<a class="anchor" aria-label="anchor" href="#spatial-operations"><i class="fas fa-link"></i></a>
</h1>
<div id="prerequisites-2" class="section level2 unnumbered">
<h2>Prerequisites<a class="anchor" aria-label="anchor" href="#prerequisites-2"><i class="fas fa-link"></i></a>
</h2>
<ul>
<li>This chapter requires the same packages used in Chapter <a href="attr.html#attr">3</a>:</li>
</ul>
<div class="sourceCode" id="cb111"><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>
<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></code></pre></div>
</div>
<div id="introduction-1" class="section level2" number="4.1">
<h2>
<span class="header-section-number">4.1</span> Introduction<a class="anchor" aria-label="anchor" href="#introduction-1"><i class="fas fa-link"></i></a>
</h2>
<p>Spatial operations, including spatial joins between vector datasets and local and focal operations on raster datasets, are a vital part of geocomputation.
This chapter shows how spatial objects can be modified in a multitude of ways based on their location and shape.
Many spatial operations have a non-spatial (attribute) equivalent, so concepts such as subsetting and joining datasets demonstrated in the previous chapter are applicable here.
This is especially true for <em>vector</em> operations: Section <a href="attr.html#vector-attribute-manipulation">3.2</a> on vector attribute manipulation provides the basis for understanding its spatial counterpart, namely spatial subsetting (covered in Section <a href="spatial-operations.html#spatial-subsetting">4.2.1</a>).
Spatial joining (Sections <a href="spatial-operations.html#spatial-joining">4.2.5</a>, <a href="spatial-operations.html#non-overlapping-joins">4.2.6</a> and <a href="spatial-operations.html#incongruent">4.2.8</a>) and aggregation (Section <a href="spatial-operations.html#spatial-aggr">4.2.7</a>) also have non-spatial counterparts, covered in the previous chapter.</p>
<p>Spatial operations differ from non-spatial operations in a number of ways, however:
spatial joins, for example, can be done in a number of ways — including matching entities that intersect with or are within a certain distance of the target dataset — while the attribution joins discussed in Section <a href="attr.html#vector-attribute-joining">3.2.4</a> in the previous chapter can only be done in one way (except when using fuzzy joins, as described in the documentation of the <a href="https://cran.r-project.org/package=fuzzyjoin"><strong>fuzzyjoin</strong></a> package).
Different <em>types</em> of spatial relationship between objects, including intersects and disjoint, are described in Sections <a href="spatial-operations.html#topological-relations">4.2.2</a> and <a href="spatial-operations.html#DE-9IM-strings">4.2.4</a>.
Another unique aspect of spatial objects is distance: all spatial objects are related through space, and distance calculations can be used to explore the strength of this relationship, as described in the context of vector data in Section <a href="spatial-operations.html#distance-relations">4.2.3</a>.</p>
<p>Spatial operations on raster objects include subsetting — covered in Section <a href="spatial-operations.html#spatial-raster-subsetting">4.3.1</a>.
<em>Map algebra</em> covers a range of operations that modify raster cell values, with or without reference to surrounding cell values.
The concept of map algebra, vital for many applications, is introduced in Section <a href="spatial-operations.html#map-algebra">4.3.2</a>; local, focal and zonal map algebra operations are covered in sections <a href="spatial-operations.html#local-operations">4.3.3</a>, <a href="spatial-operations.html#focal-operations">4.3.4</a>, and <a href="spatial-operations.html#zonal-operations">4.3.5</a>, respectively.
Global map algebra operations, which generate summary statistics representing an entire raster dataset, and distance calculations on rasters, are discussed in Section <a href="spatial-operations.html#global-operations-and-distances">4.3.6</a>.
Next, relation between map algebra and vector operations are discussed in Section <a href="spatial-operations.html#map-algebra-counterparts-in-vector-processing">4.3.7</a>.
In the final section before the exercises (<a href="spatial-operations.html#merging-rasters">4.3.8</a>) the process of merging two raster datasets is discussed and demonstrated with reference to a reproducible example.</p>
<div class="rmdnote">
It is important to note that spatial operations that use two spatial objects rely on both objects having the same coordinate reference system, a topic that was introduced in Section <a href="spatial-class.html#crs-intro">2.4</a> and which will be covered in more depth in Chapter <a href="reproj-geo-data.html#reproj-geo-data">7</a>.
</div>
</div>
<div id="spatial-vec" class="section level2" number="4.2">
<h2>
<span class="header-section-number">4.2</span> Spatial operations on vector data<a class="anchor" aria-label="anchor" href="#spatial-vec"><i class="fas fa-link"></i></a>
</h2>
<p>This section provides an overview of spatial operations on vector geographic data represented as simple features in the <strong>sf</strong> package.
Section <a href="spatial-operations.html#spatial-ras">4.3</a> presents spatial operations on raster datasets using classes and functions from the <strong>terra</strong> package.</p>
<div id="spatial-subsetting" class="section level3" number="4.2.1">
<h3>
<span class="header-section-number">4.2.1</span> Spatial subsetting<a class="anchor" aria-label="anchor" href="#spatial-subsetting"><i class="fas fa-link"></i></a>
</h3>
<p>Spatial subsetting is the process of taking a spatial object and returning a new object containing only features that <em>relate</em> in space to another object.
Analogous to <em>attribute subsetting</em> (covered in Section <a href="attr.html#vector-attribute-subsetting">3.2.1</a>), subsets of <code>sf</code> data frames can be created with square bracket (<code>[</code>) operator using the syntax <code>x[y, , op = st_intersects]</code>, where <code>x</code> is an <code>sf</code> object from which a subset of rows will be returned, <code>y</code> is the ‘subsetting object’ and <code>, op = st_intersects</code> is an optional argument that specifies the topological relation (also known as the binary predicate) used to do the subsetting.
The default topological relation used when an <code>op</code> argument is not provided is <code><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_intersects()</a></code>: the command <code>x[y, ]</code> is identical to <code>x[y, , op = st_intersects]</code> shown above but not <code>x[y, , op = st_disjoint]</code> (the meaning of these and other topological relations is described in the next section).
The <code><a href="https://dplyr.tidyverse.org/reference/filter.html">filter()</a></code> function from the <strong>tidyverse</strong> can also be used but this approach is more verbose, as we will see in the examples below.
</p>
<p>To demonstrate spatial subsetting, we will use the <code>nz</code> and <code>nz_height</code> datasets in the <strong>spData</strong> package, which contain geographic data on the 16 main regions and 101 highest points in New Zealand, respectively (Figure <a href="spatial-operations.html#fig:nz-subset">4.1</a>), in a projected coordinate reference system.
The following code chunk creates an object representing Canterbury, then uses spatial subsetting to return all high points in the region.</p>
<div class="sourceCode" id="cb112"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">canterbury</span> <span class="op">=</span> <span class="va">nz</span> <span class="op">|></span> <span class="fu"><a href="https://dplyr.tidyverse.org/reference/filter.html">filter</a></span><span class="op">(</span><span class="va">Name</span> <span class="op">==</span> <span class="st">"Canterbury"</span><span class="op">)</span></span>
<span><span class="va">canterbury_height</span> <span class="op">=</span> <span class="va">nz_height</span><span class="op">[</span><span class="va">canterbury</span>, <span class="op">]</span></span></code></pre></div>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:nz-subset"></span>
<img src="figures/nz-subset-1.png" alt="Illustration of spatial subsetting with red triangles representing 101 high points in New Zealand, clustered near the central Canterbuy region (left). The points in Canterbury were created with the `[` subsetting operator (highlighted in gray, right)." width="100%"><p class="caption">
FIGURE 4.1: Illustration of spatial subsetting with red triangles representing 101 high points in New Zealand, clustered near the central Canterbuy region (left). The points in Canterbury were created with the <code>[</code> subsetting operator (highlighted in gray, right).
</p>
</div>
<p>Like attribute subsetting, the command <code>x[y, ]</code> (equivalent to <code>nz_height[canterbury, ]</code>) subsets features of a <em>target</em> <code>x</code> using the contents of a <em>source</em> object <code>y</code>.
Instead of <code>y</code> being a vector of class <code>logical</code> or <code>integer</code>, however, for spatial subsetting both <code>x</code> and <code>y</code> must be geographic objects.
Specifically, objects used for spatial subsetting in this way must have the class <code>sf</code> or <code>sfc</code>: both <code>nz</code> and <code>nz_height</code> are geographic vector data frames and have the class <code>sf</code>, and the result of the operation returns another <code>sf</code> object representing the features in the target <code>nz_height</code> object that intersect with (in this case high points that are located within) the <code>canterbury</code> region.</p>
<p>Various <em>topological relations</em> can be used for spatial subsetting which determine the type of spatial relationship that features in the target object must have with the subsetting object to be selected.
These include <em>touches</em>, <em>crosses</em> or <em>within</em>, as we will see shortly in Section <a href="spatial-operations.html#topological-relations">4.2.2</a>.
The default setting <code>st_intersects</code> is a ‘catch all’ topological relation that will return features in the target that <em>touch</em>, <em>cross</em> or are <em>within</em> the source ‘subsetting’ object.
Alternative spatial operators can be specified with the <code>op =</code> argument, as demonstrated in the following command which returns the opposite of <code><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_intersects()</a></code>, points that do not intersect with Canterbury (see Section <a href="spatial-operations.html#topological-relations">4.2.2</a>).</p>
<div class="sourceCode" id="cb113"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">nz_height</span><span class="op">[</span><span class="va">canterbury</span>, , op <span class="op">=</span> <span class="va">st_disjoint</span><span class="op">]</span></span></code></pre></div>
<div class="rmdnote">
Note the empty argument — denoted with <code>, ,</code> — in the preceding code chunk is included to highlight <code>op</code>, the third argument in <code>[</code> for <code>sf</code> objects.
One can use this to change the subsetting operation in many ways.
<code>nz_height[canterbury, 2, op = st_disjoint]</code>, for example, returns the same rows but only includes the second attribute column (see <code><a href="https://r-spatial.github.io/sf/reference/sf.html">sf::`[.sf`</a></code> and the <code><a href="https://r-spatial.github.io/sf/reference/sf.html">?sf</a></code> for details).
</div>
<p>For many applications, this is all you’ll need to know about spatial subsetting for vector data: it just works.
If you are impatient to learn about more topological relations, beyond <code><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_intersects()</a></code> and <code><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_disjoint()</a></code>, skip to the next section (<a href="spatial-operations.html#topological-relations">4.2.2</a>).
If you’re interested in the details, including other ways of subsetting, read on.</p>
<p>Another way of doing spatial subsetting uses objects returned by topological operators.
These objects can be useful in their own right, for example when exploring the graph network of relationships between contiguous regions, but they can also be used for subsetting, as demonstrated in the code chunk below.</p>
<div class="sourceCode" id="cb114"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">sel_sgbp</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_intersects</a></span><span class="op">(</span>x <span class="op">=</span> <span class="va">nz_height</span>, y <span class="op">=</span> <span class="va">canterbury</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/class.html">class</a></span><span class="op">(</span><span class="va">sel_sgbp</span><span class="op">)</span></span>
<span><span class="co">#> [1] "sgbp" "list"</span></span>
<span><span class="va">sel_sgbp</span></span>
<span><span class="co">#> Sparse geometry binary predicate list of length 101, where the</span></span>
<span><span class="co">#> predicate was `intersects'</span></span>
<span><span class="co">#> first 10 elements:</span></span>
<span><span class="co">#> 1: (empty)</span></span>
<span><span class="co">#> 2: (empty)</span></span>
<span><span class="co">#> 3: (empty)</span></span>
<span><span class="co">#> 4: (empty)</span></span>
<span><span class="co">#> 5: 1</span></span>
<span><span class="co">#> 6: 1</span></span>
<span><span class="va">....</span></span>
<span><span class="va">sel_logical</span> <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/lengths.html">lengths</a></span><span class="op">(</span><span class="va">sel_sgbp</span><span class="op">)</span> <span class="op">></span> <span class="fl">0</span></span>
<span><span class="va">canterbury_height2</span> <span class="op">=</span> <span class="va">nz_height</span><span class="op">[</span><span class="va">sel_logical</span>, <span class="op">]</span></span></code></pre></div>
<p>The above code chunk creates an object of class <code>sgbp</code> (a sparse geometry binary predicate, a list of length <code>x</code> in the spatial operation) and then converts it into a logical vector <code>sel_logical</code> (containing only <code>TRUE</code> and <code>FALSE</code> values, something that can also be used by <strong>dplyr</strong>’s filter function).
The function <code><a href="https://rdrr.io/r/base/lengths.html">lengths()</a></code> identifies which features in <code>nz_height</code> intersect with <em>any</em> objects in <code>y</code>.
In this case 1 is the greatest possible value but for more complex operations one could use the method to subset only features that intersect with, for example, 2 or more features from the source object.</p>
<div class="rmdnote">
Note: another way to return a logical output is by setting <code>sparse = FALSE</code> (meaning ‘return a dense matrix not a sparse one’) in operators such as <code><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_intersects()</a></code>. The command <code>st_intersects(x = nz_height, y = canterbury, sparse = FALSE)[, 1]</code>, for example, would return an output identical to <code>sel_logical</code>.
Note: the solution involving <code>sgbp</code> objects is more generalisable though, as it works for many-to-many operations and has lower memory requirements.
</div>
<p>The same result can be also achieved with the <strong>sf</strong> function <code><a href="https://r-spatial.github.io/sf/reference/st_join.html">st_filter()</a></code> which was <a href="https://github.com/r-spatial/sf/issues/1148">created</a> to increase compatibility between <code>sf</code> objects and <strong>dplyr</strong> data manipulation code:</p>
<div class="sourceCode" id="cb115"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">canterbury_height3</span> <span class="op">=</span> <span class="va">nz_height</span> <span class="op">|></span></span>
<span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_join.html">st_filter</a></span><span class="op">(</span>y <span class="op">=</span> <span class="va">canterbury</span>, .predicate <span class="op">=</span> <span class="va">st_intersects</span><span class="op">)</span></span></code></pre></div>
<p>At this point, there are three identical (in all but row names) versions of <code>canterbury_height</code>, one created using the <code>[</code> operator, one created via an intermediary selection object, and another using <strong>sf</strong>’s convenience function <code><a href="https://r-spatial.github.io/sf/reference/st_join.html">st_filter()</a></code>.
<!-- RL: commented out for now as old. Todo: if we ever update that vignette uncomment the next line. -->
<!-- To explore spatial subsetting in more detail, see the supplementary vignettes on `subsetting` and [`tidyverse-pitfalls`](https://geocompr.github.io/geocompkg/articles/) on the [geocompkg website](https://geocompr.github.io/geocompkg/articles/). -->
The next section explores different types of spatial relation, also known as binary predicates, that can be used to identify whether or not two features are spatially related or not.</p>
</div>
<div id="topological-relations" class="section level3" number="4.2.2">
<h3>
<span class="header-section-number">4.2.2</span> Topological relations<a class="anchor" aria-label="anchor" href="#topological-relations"><i class="fas fa-link"></i></a>
</h3>
<p>Topological relations describe the spatial relationships between objects.
“Binary topological relationships”, to give them their full name, are logical statements (in that the answer can only be <code>TRUE</code> or <code>FALSE</code>) about the spatial relationships between two objects defined by ordered sets of points (typically forming points, lines and polygons) in two or more dimensions <span class="citation">(<a href="references.html#ref-egenhofer_mathematical_1990">Egenhofer and Herring 1990</a>)</span>.
That may sound rather abstract and, indeed, the definition and classification of topological relations is based on mathematical foundations first published in book form in 1966 <span class="citation">(<a href="references.html#ref-spanier_algebraic_1995">Spanier 1995</a>)</span>, with the field of algebraic topology continuing into the 21<sup>st</sup> century <span class="citation">(<a href="references.html#ref-dieck_algebraic_2008">Dieck 2008</a>)</span>.</p>
<p>Despite their mathematical origins, topological relations can be understood intuitively with reference to visualizations of commonly used functions that test for common types of spatial relationships.
Figure <a href="spatial-operations.html#fig:relations">4.2</a> shows a variety of geometry pairs and their associated relations.
The third and fourth pairs in Figure <a href="spatial-operations.html#fig:relations">4.2</a> (from left to right and then down) demonstrate that, for some relations, order is important.
While the relations <em>equals</em>, <em>intersects</em>, <em>crosses</em>, <em>touches</em> and <em>overlaps</em> are symmetrical, meaning that if <code>function(x, y)</code> is true, <code>function(y, x)</code> will also be true, relations in which the order of the geometries are important such as <em>contains</em> and <em>within</em> are not.
Notice that each geometry pair has a “DE-9IM” string such as FF2F11212, described in the next section.
</p>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:relations"></span>
<img src="figures/relations-1.png" alt="Topological relations between vector geometries, inspired by Figures 1 and 2 in Egenhofer and Herring (1990). The relations for which the function(x, y) is true are printed for each geometry pair, with x represented in pink and y represented in blue. The nature of the spatial relationship for each pair is described by the Dimensionally Extended 9-Intersection Model string." width="100%"><p class="caption">
FIGURE 4.2: Topological relations between vector geometries, inspired by Figures 1 and 2 in Egenhofer and Herring (1990). The relations for which the function(x, y) is true are printed for each geometry pair, with x represented in pink and y represented in blue. The nature of the spatial relationship for each pair is described by the Dimensionally Extended 9-Intersection Model string.
</p>
</div>
<p>In <code>sf</code>, functions testing for different types of topological relations are called ‘binary predicates’, as described in the vignette <em>Manipulating Simple Feature Geometries</em>, which can be viewed with the command <a href="https://r-spatial.github.io/sf/articles/sf3.html"><code>vignette("sf3")</code></a>, and in the help page <a href="https://r-spatial.github.io/sf/reference/geos_binary_ops.html"><code>?geos_binary_pred</code></a>.
To see how topological relations work in practice, let’s create a simple reproducible example, building on the relations illustrated in Figure <a href="spatial-operations.html#fig:relations">4.2</a> and consolidating knowledge of how vector geometries are represented from a previous chapter (Section <a href="spatial-class.html#geometry">2.2.4</a>).
Note that to create tabular data representing coordinates (x and y) of the polygon vertices, we use the base R function <code><a href="https://rdrr.io/r/base/cbind.html">cbind()</a></code> to create a matrix representing coordinates points, a <code>POLYGON</code>, and finally an <code>sfc</code> object, as described in Chapter <a href="spatial-class.html#spatial-class">2</a>:</p>
<div class="sourceCode" id="cb116"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">polygon_matrix</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>
<span> x <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">0</span>, <span class="fl">0</span>, <span class="fl">1</span>, <span class="fl">1</span>, <span class="fl">0</span><span class="op">)</span>,</span>
<span> y <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">0</span>, <span class="fl">1</span>, <span class="fl">1</span>, <span class="fl">0.5</span>, <span class="fl">0</span><span class="op">)</span></span>
<span><span class="op">)</span></span>
<span><span class="va">polygon_sfc</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/sfc.html">st_sfc</a></span><span class="op">(</span><span class="fu"><a href="https://r-spatial.github.io/sf/reference/st.html">st_polygon</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/list.html">list</a></span><span class="op">(</span><span class="va">polygon_matrix</span><span class="op">)</span><span class="op">)</span><span class="op">)</span></span></code></pre></div>
<p>We will create additional geometries to demonstrate spatial relations with the following commands which, when plotted on top of the polygon created above, relate in space to one another, as shown in Figure <a href="spatial-operations.html#fig:relation-objects">4.3</a>.
Note the use of the function <code><a href="https://r-spatial.github.io/sf/reference/st_as_sf.html">st_as_sf()</a></code> and the argument <code>coords</code> to efficiently convert from a data frame containing columns representing coordinates to an <code>sf</code> object containing points:</p>
<div class="sourceCode" id="cb117"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">line_sfc</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/sfc.html">st_sfc</a></span><span class="op">(</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="fu"><a href="https://rdrr.io/r/base/cbind.html">cbind</a></span><span class="op">(</span></span>
<span> x <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">0.4</span>, <span class="fl">1</span><span class="op">)</span>,</span>
<span> y <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">0.2</span>, <span class="fl">0.5</span><span class="op">)</span></span>
<span><span class="op">)</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="co"># create points</span></span>
<span><span class="va">point_df</span> <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/data.frame.html">data.frame</a></span><span class="op">(</span></span>
<span> x <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">0.2</span>, <span class="fl">0.7</span>, <span class="fl">0.4</span><span class="op">)</span>,</span>
<span> y <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">0.1</span>, <span class="fl">0.2</span>, <span class="fl">0.8</span><span class="op">)</span></span>
<span><span class="op">)</span></span>
<span><span class="va">point_sf</span> <span class="op">=</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="va">point_df</span>, coords <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="st">"x"</span>, <span class="st">"y"</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:relation-objects"></span>
<img src="figures/relation-objects-1.png" alt="Points, line and polygon objects arranged to illustrate topological relations." width="50%"><p class="caption">
FIGURE 4.3: Points, line and polygon objects arranged to illustrate topological relations.
</p>
</div>
<p>A simple query is: which of the points in <code>point_sf</code> intersect in some way with polygon <code>polygon_sfc</code>?
The question can be answered by inspection (points 1 and 3 are touching and within the polygon, respectively).
This question can be answered with the spatial predicate <code><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_intersects()</a></code> as follows:</p>
<div class="sourceCode" id="cb118"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_intersects</a></span><span class="op">(</span><span class="va">point_sf</span>, <span class="va">polygon_sfc</span><span class="op">)</span></span>
<span><span class="co">#> Sparse geometry binary predicate... `intersects'</span></span>
<span><span class="co">#> 1: 1</span></span>
<span><span class="co">#> 2: (empty)</span></span>
<span><span class="co">#> 3: 1</span></span></code></pre></div>
<p>The result should match your intuition:
positive (<code>1</code>) results are returned for the first and third point, and a negative result (represented by an empty vector) for the second are outside the polygon’s border.
What may be unexpected is that the result comes in the form of a list of vectors.
This <em>sparse matrix</em> output only registers a relation if one exists, reducing the memory requirements of topological operations on multi-feature objects.
As we saw in the previous section, a <em>dense matrix</em> consisting of <code>TRUE</code> or <code>FALSE</code> values is returned when <code>sparse = FALSE</code>.</p>
<div class="sourceCode" id="cb119"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_intersects</a></span><span class="op">(</span><span class="va">point_sf</span>, <span class="va">polygon_sfc</span>, sparse <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span></span>
<span><span class="co">#> [,1]</span></span>
<span><span class="co">#> [1,] TRUE</span></span>
<span><span class="co">#> [2,] FALSE</span></span>
<span><span class="co">#> [3,] TRUE</span></span></code></pre></div>
<p>In the above output each row represents a feature in the target (argument <code>x</code>) object and each column represents a feature in the selecting object (<code>y</code>).
In this case, there is only one feature in the <code>y</code> object <code>polygon_sfc</code> so the result, which can be used for subsetting as we saw in Section <a href="spatial-operations.html#spatial-subsetting">4.2.1</a>, has only one column.</p>
<p><code><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_intersects()</a></code> returns <code>TRUE</code> even in cases where the features just touch: <em>intersects</em> is a ‘catch-all’ topological operation which identifies many types of spatial relation, as illustrated in Figure <a href="spatial-operations.html#fig:relations">4.2</a>.
More restrictive questions include which points lie within the polygon, and which features are on or contain a shared boundary with <code>y</code>?
These can be answered as follows (results not shown):</p>
<div class="sourceCode" id="cb120"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_within</a></span><span class="op">(</span><span class="va">point_sf</span>, <span class="va">polygon_sfc</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_touches</a></span><span class="op">(</span><span class="va">point_sf</span>, <span class="va">polygon_sfc</span><span class="op">)</span></span></code></pre></div>
<p>Note that although the first point <em>touches</em> the boundary polygon, it is not within it; the third point is within the polygon but does not touch any part of its border.
The opposite of <code><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_intersects()</a></code> is <code><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_disjoint()</a></code>, which returns only objects that do not spatially relate in any way to the selecting object (note <code>[, 1]</code> converts the result into a vector).</p>
<div class="sourceCode" id="cb121"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_disjoint</a></span><span class="op">(</span><span class="va">point_sf</span>, <span class="va">polygon_sfc</span>, sparse <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span><span class="op">[</span>, <span class="fl">1</span><span class="op">]</span></span>
<span><span class="co">#> [1] FALSE TRUE FALSE</span></span></code></pre></div>
<p>The function <code><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_is_within_distance()</a></code> detects features that <em>almost touch</em> the selection object, which has an additional <code>dist</code> argument.
It can be used to set how close target objects need to be before they are selected.
The ‘is within distance’ binary spatial predicate is demonstrated in the code chunk below, the results of which show that every point is within 0.2 units of the polygon.</p>
<div class="sourceCode" id="cb122"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_is_within_distance</a></span><span class="op">(</span><span class="va">point_sf</span>, <span class="va">polygon_sfc</span>, dist <span class="op">=</span> <span class="fl">0.2</span>, sparse <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span><span class="op">[</span>, <span class="fl">1</span><span class="op">]</span></span>
<span><span class="co">#> [1] TRUE TRUE TRUE</span></span></code></pre></div>
<p>Note that although point 2 is more than 0.2 units of distance from the nearest vertex of <code>polygon_sfc</code>, it is still selected when the distance is set to 0.2.
This is because distance is measured to the nearest edge, in this case the part of the polygon that lies directly above point 2 in Figure <a href="spatial-operations.html#fig:relation-objects">4.3</a>.
(You can verify the actual distance between point 2 and the polygon is 0.13 with the command <code>st_distance(point_sf, polygon_sfc)</code>.)</p>
<div class="rmdnote">
Functions for calculating topological relations use spatial indices to largely speed up spatial query performance.
They achieve that using the Sort-Tile-Recursive (STR) algorithm.
The <code>st_join</code> function, mentioned in the next section, also uses the spatial indexing.
You can learn more at <a href="https://www.r-spatial.org/r/2017/06/22/spatial-index.html" class="uri">https://www.r-spatial.org/r/2017/06/22/spatial-index.html</a>.
</div>
</div>
<div id="distance-relations" class="section level3" number="4.2.3">
<h3>
<span class="header-section-number">4.2.3</span> Distance relations<a class="anchor" aria-label="anchor" href="#distance-relations"><i class="fas fa-link"></i></a>
</h3>
<p>While the topological relations presented in the previous section are binary (a feature either intersects with another or does not) distance relations are continuous.
The distance between two <code>sf</code> objects is calculated with <code><a href="https://r-spatial.github.io/sf/reference/geos_measures.html">st_distance()</a></code>, which is also used behind the scenes in Section <a href="spatial-operations.html#non-overlapping-joins">4.2.6</a> for distance-based joins.
This is illustrated in the code chunk below, which finds the distance between the highest point in New Zealand and the geographic centroid of the Canterbury region, created in Section <a href="spatial-operations.html#spatial-subsetting">4.2.1</a>:
</p>
<div class="sourceCode" id="cb123"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">nz_highest</span> <span class="op">=</span> <span class="va">nz_height</span> <span class="op">|></span> <span class="fu"><a href="https://dplyr.tidyverse.org/reference/slice.html">slice_max</a></span><span class="op">(</span>n <span class="op">=</span> <span class="fl">1</span>, order_by <span class="op">=</span> <span class="va">elevation</span><span class="op">)</span></span>
<span><span class="va">canterbury_centroid</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/geos_unary.html">st_centroid</a></span><span class="op">(</span><span class="va">canterbury</span><span class="op">)</span></span>
<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">nz_highest</span>, <span class="va">canterbury_centroid</span><span class="op">)</span></span>
<span><span class="co">#> Units: [m]</span></span>
<span><span class="co">#> [,1]</span></span>
<span><span class="co">#> [1,] 115540</span></span></code></pre></div>
<p>There are two potentially surprising things about the result:</p>
<ul>
<li>It has <code>units</code>, telling us the distance is 100,000 meters, not 100,000 inches, or any other measure of distance</li>
<li>It is returned as a matrix, even though the result only contains a single value</li>
</ul>
<p>This second feature hints at another useful feature of <code><a href="https://r-spatial.github.io/sf/reference/geos_measures.html">st_distance()</a></code>, its ability to return <em>distance matrices</em> between all combinations of features in objects <code>x</code> and <code>y</code>.
This is illustrated in the command below, which finds the distances between the first three features in <code>nz_height</code> and the Otago and Canterbury regions of New Zealand represented by the object <code>co</code>.</p>
<div class="sourceCode" id="cb124"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">co</span> <span class="op">=</span> <span class="fu"><a href="https://dplyr.tidyverse.org/reference/filter.html">filter</a></span><span class="op">(</span><span class="va">nz</span>, <span class="fu"><a href="https://rdrr.io/r/base/grep.html">grepl</a></span><span class="op">(</span><span class="st">"Canter|Otag"</span>, <span class="va">Name</span><span class="op">)</span><span class="op">)</span></span>
<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">nz_height</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">3</span>, <span class="op">]</span>, <span class="va">co</span><span class="op">)</span></span>
<span><span class="co">#> Units: [m]</span></span>
<span><span class="co">#> [,1] [,2]</span></span>
<span><span class="co">#> [1,] 123537 15498</span></span>
<span><span class="co">#> [2,] 94283 0</span></span>
<span><span class="co">#> [3,] 93019 0</span></span></code></pre></div>
<p>Note that the distance between the second and third features in <code>nz_height</code> and the second feature in <code>co</code> is zero.
This demonstrates the fact that distances between points and polygons refer to the distance to <em>any part of the polygon</em>:
The second and third points in <code>nz_height</code> are <em>in</em> Otago, which can be verified by plotting them (result not shown):</p>
<div class="sourceCode" id="cb125"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rspatial.github.io/terra/reference/plot.html">plot</a></span><span class="op">(</span><span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_geometry.html">st_geometry</a></span><span class="op">(</span><span class="va">co</span><span class="op">)</span><span class="op">[</span><span class="fl">2</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="fu"><a href="https://r-spatial.github.io/sf/reference/st_geometry.html">st_geometry</a></span><span class="op">(</span><span class="va">nz_height</span><span class="op">)</span><span class="op">[</span><span class="fl">2</span><span class="op">:</span><span class="fl">3</span><span class="op">]</span>, add <span class="op">=</span> <span class="cn">TRUE</span><span class="op">)</span></span></code></pre></div>
</div>
<div id="DE-9IM-strings" class="section level3" number="4.2.4">
<h3>
<span class="header-section-number">4.2.4</span> DE-9IM strings<a class="anchor" aria-label="anchor" href="#DE-9IM-strings"><i class="fas fa-link"></i></a>
</h3>
<p>Underlying the binary predicates demonstrated in the previous section is the Dimensionally Extended 9-Intersection Model (DE-9IM).
As the cryptic name suggests, this is not an easy topic to understand, but it is worth knowing about because it underlies many spatial operations and enables the creation of custom spatial predicates.
The model was originally labelled “DE + 9IM” by its inventors, referring to the “dimension of the intersections of boundaries, interiors, and exteriors of two features” <span class="citation">(<a href="references.html#ref-clementini_comparison_1995">Clementini and Di Felice 1995</a>)</span>, but is now referred to as DE-9IM <span class="citation">(<a href="references.html#ref-shen_classification_2018">Shen, Chen, and Liu 2018</a>)</span>.
DE-9IM is applicable to 2-dimensional objects (points, lines and polygons) in Euclidean space, meaning that the model (and software implementing it such as GEOS) assumes you are working with data in a projected coordinate reference system, described in Chapter <a href="reproj-geo-data.html#reproj-geo-data">7</a>.</p>
<p>To demonstrate how DE-9IM strings work, let’s take a look at the various ways that the first geometry pair in Figure <a href="spatial-operations.html#fig:relations">4.2</a> relate.
Figure <a href="spatial-operations.html#fig:de9imgg">4.4</a> illustrates the 9 intersection model (9IM) which shows the intersections between every combination of each object’s interior, boundary and exterior: when each component of the first object <code>x</code> is arranged as columns and each component of <code>y</code> is arranged as rows, a facetted graphic is created with the intersections between each element highlighted.</p>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:de9imgg"></span>
<img src="figures/de9imgg-1.png" alt="Illustration of how the Dimensionally Extended 9 Intersection Model (DE-9IM) works. Colors not in the legend represent the overlap between different components. The thick lines highlight 2 dimensional intesections, e.g., between the boundary of object x and the interior of object y, shown in the middle top facet." width="100%"><p class="caption">
FIGURE 4.4: Illustration of how the Dimensionally Extended 9 Intersection Model (DE-9IM) works. Colors not in the legend represent the overlap between different components. The thick lines highlight 2 dimensional intesections, e.g., between the boundary of object x and the interior of object y, shown in the middle top facet.
</p>
</div>
<p>DE-9IM strings are derived from the dimension of each type of relation.
In this case the red intersections in Figure <a href="spatial-operations.html#fig:de9imgg">4.4</a> have dimensions of 0 (points), 1 (lines), and 2 (polygons), as shown in Table <a href="spatial-operations.html#tab:de9emtable">4.1</a>.</p>
<div class="inline-table"><table class="table table-sm">
<caption>
<span id="tab:de9emtable">TABLE 4.1: </span>Table showing relations between interiors, boundaries and exteriors of geometries x and y.</caption>
<thead><tr class="header">
<th align="left"></th>
<th align="left">Interior (x)</th>
<th align="left">Boundary (x)</th>
<th align="left">Exterior (x)</th>
</tr></thead>
<tbody>
<tr class="odd">
<td align="left">Interior (y)</td>
<td align="left">2</td>
<td align="left">1</td>
<td align="left">2</td>
</tr>
<tr class="even">
<td align="left">Boundary (y)</td>
<td align="left">1</td>
<td align="left">1</td>
<td align="left">1</td>
</tr>
<tr class="odd">
<td align="left">Exterior (y)</td>
<td align="left">2</td>
<td align="left">1</td>
<td align="left">2</td>
</tr>
</tbody>
</table></div>
<p>Flattening this matrix ‘row-wise’ (meaning concatenating the first row, then the second, then the third) results in the string <code>212111212</code>.
Another example will serve to demonstrate the system:
the relation shown in Figure <a href="spatial-operations.html#fig:relations">4.2</a> (the third polygon pair in the third column and 1st row) can be defined in the DE-9IM system as follows:</p>
<ul>
<li>The intersections between the <em>interior</em> of the larger object <code>x</code> and the interior, boundary and exterior of <code>y</code> have dimensions of 2, 1 and 2 respectively</li>
<li>The intersections between the <em>boundary</em> of the larger object <code>x</code> and the interior, boundary and exterior of <code>y</code> have dimensions of F, F and 1 respectively, where ‘F’ means ‘false’, the objects are disjoint</li>
<li>The intersections between the <em>exterior</em> of <code>x</code> and the interior, boundary and exterior of <code>y</code> have dimensions of F, F and 2 respectively: the exterior of the larger object does not touch the interior or boundary of <code>y</code>, but the exterior of the smaller and larger objects cover the same area</li>
</ul>
<p>These three components, when concatenated, create the string <code>212</code>, <code>FF1</code>, and <code>FF2</code>.
This is the same as the result obtained from the function <code><a href="https://r-spatial.github.io/sf/reference/st_relate.html">st_relate()</a></code> (see the source code of this chapter to see how other geometries in Figure <a href="spatial-operations.html#fig:relations">4.2</a> were created):</p>
<div class="sourceCode" id="cb126"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">xy2sfc</span> <span class="op">=</span> <span class="kw">function</span><span class="op">(</span><span class="va">x</span>, <span class="va">y</span><span class="op">)</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/sfc.html">st_sfc</a></span><span class="op">(</span><span class="fu"><a href="https://r-spatial.github.io/sf/reference/st.html">st_polygon</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/list.html">list</a></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">x</span>, <span class="va">y</span><span class="op">)</span><span class="op">)</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="va">x</span> <span class="op">=</span> <span class="fu">xy2sfc</span><span class="op">(</span>x <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">0</span>, <span class="fl">0</span>, <span class="fl">1</span>, <span class="fl">1</span>, <span class="fl">0</span><span class="op">)</span>, y <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">0</span>, <span class="fl">1</span>, <span class="fl">1</span>, <span class="fl">0.5</span>, <span class="fl">0</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="va">y</span> <span class="op">=</span> <span class="fu">xy2sfc</span><span class="op">(</span>x <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">0.7</span>, <span class="fl">0.7</span>, <span class="fl">0.9</span>, <span class="fl">0.7</span><span class="op">)</span>, y <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">0.8</span>, <span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">0.8</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_relate.html">st_relate</a></span><span class="op">(</span><span class="va">x</span>, <span class="va">y</span><span class="op">)</span></span>
<span><span class="co">#> [,1] </span></span>
<span><span class="co">#> [1,] "212FF1FF2"</span></span></code></pre></div>
<p>Understanding DE-9IM strings allows new binary spatial predicates to be developed.
The help page <code><a href="https://r-spatial.github.io/sf/reference/st_relate.html">?st_relate</a></code> contains function definitions for ‘queen’ and ‘rook’ relations in which polygons share a border or only a point, respectively.
‘Queen’ relations mean that ‘boundary-boundary’ relations (the cell in the second column and the second row in Table <a href="spatial-operations.html#tab:de9emtable">4.1</a>, or the 5th element of the DE-9IM string) must not be empty, corresponding to the pattern <code>F***T****</code>, while for ‘rook’ relations the same element must be 1 (meaning a linear intersection).
These are implemented as follows:</p>
<div class="sourceCode" id="cb127"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">st_queen</span> <span class="op">=</span> <span class="kw">function</span><span class="op">(</span><span class="va">x</span>, <span class="va">y</span><span class="op">)</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_relate.html">st_relate</a></span><span class="op">(</span><span class="va">x</span>, <span class="va">y</span>, pattern <span class="op">=</span> <span class="st">"F***T****"</span><span class="op">)</span></span>
<span><span class="va">st_rook</span> <span class="op">=</span> <span class="kw">function</span><span class="op">(</span><span class="va">x</span>, <span class="va">y</span><span class="op">)</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_relate.html">st_relate</a></span><span class="op">(</span><span class="va">x</span>, <span class="va">y</span>, pattern <span class="op">=</span> <span class="st">"F***1****"</span><span class="op">)</span></span></code></pre></div>
<p>Building on the object <code>x</code> created previously, we can use the newly created functions to find out which elements in the grid are a ‘queen’ and ‘rook’ in relation to the middle square of the grid as follows:</p>
<div class="sourceCode" id="cb128"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">grid</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_make_grid.html">st_make_grid</a></span><span class="op">(</span><span class="va">x</span>, n <span class="op">=</span> <span class="fl">3</span><span class="op">)</span></span>
<span><span class="va">grid_sf</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/sf.html">st_sf</a></span><span class="op">(</span><span class="va">grid</span><span class="op">)</span></span>
<span><span class="va">grid_sf</span><span class="op">$</span><span class="va">queens</span> <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/lengths.html">lengths</a></span><span class="op">(</span><span class="fu">st_queen</span><span class="op">(</span><span class="va">grid</span>, <span class="va">grid</span><span class="op">[</span><span class="fl">5</span><span class="op">]</span><span class="op">)</span><span class="op">)</span> <span class="op">></span> <span class="fl">0</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">grid</span>, col <span class="op">=</span> <span class="va">grid_sf</span><span class="op">$</span><span class="va">queens</span><span class="op">)</span></span>
<span><span class="va">grid_sf</span><span class="op">$</span><span class="va">rooks</span> <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/lengths.html">lengths</a></span><span class="op">(</span><span class="fu">st_rook</span><span class="op">(</span><span class="va">grid</span>, <span class="va">grid</span><span class="op">[</span><span class="fl">5</span><span class="op">]</span><span class="op">)</span><span class="op">)</span> <span class="op">></span> <span class="fl">0</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">grid</span>, col <span class="op">=</span> <span class="va">grid_sf</span><span class="op">$</span><span class="va">rooks</span><span class="op">)</span></span></code></pre></div>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:queens"></span>
<img src="figures/queens-1.png" alt="Demonstration of custom binary spatial predicates for finding 'queen' (left) and 'rook' (right) relations to the central square in a grid with 9 geometries." width="100%"><p class="caption">
FIGURE 4.5: Demonstration of custom binary spatial predicates for finding ‘queen’ (left) and ‘rook’ (right) relations to the central square in a grid with 9 geometries.
</p>
</div>
<!-- Another of a custom binary spatial predicate is 'overlapping lines' which detects lines that overlap for some or all of another line's geometry. -->
<!-- This can be implemented as follows, with the pattern signifying that the intersection between the two line interiors must be a line: -->
</div>
<div id="spatial-joining" class="section level3" number="4.2.5">
<h3>
<span class="header-section-number">4.2.5</span> Spatial joining<a class="anchor" aria-label="anchor" href="#spatial-joining"><i class="fas fa-link"></i></a>
</h3>
<p>Joining two non-spatial datasets relies on a shared ‘key’ variable, as described in Section <a href="attr.html#vector-attribute-joining">3.2.4</a>.
Spatial data joining applies the same concept, but instead relies on spatial relations, described in the previous section.
As with attribute data, joining adds new columns to the target object (the argument <code>x</code> in joining functions), from a source object (<code>y</code>).
</p>
<p>The process is illustrated by the following example: imagine you have ten points randomly distributed across the Earth’s surface and you ask, for the points that are on land, which countries are they in?
Implementing this idea in a <a href="https://github.com/geocompx/geocompr/blob/main/code/04-spatial-join.R">reproducible example</a> will build your geographic data handling skills and show how spatial joins work.
The starting point is to create points that are randomly scattered over the Earth’s surface.</p>
<div class="sourceCode" id="cb129"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/base/Random.html">set.seed</a></span><span class="op">(</span><span class="fl">2018</span><span class="op">)</span> <span class="co"># set seed for reproducibility</span></span>
<span><span class="op">(</span><span class="va">bb</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_bbox.html">st_bbox</a></span><span class="op">(</span><span class="va">world</span><span class="op">)</span><span class="op">)</span> <span class="co"># the world's bounds</span></span>
<span><span class="co">#> xmin ymin xmax ymax </span></span>
<span><span class="co">#> -180.0 -89.9 180.0 83.6</span></span>
<span><span class="va">random_df</span> <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/data.frame.html">data.frame</a></span><span class="op">(</span></span>
<span> x <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Uniform.html">runif</a></span><span class="op">(</span>n <span class="op">=</span> <span class="fl">10</span>, min <span class="op">=</span> <span class="va">bb</span><span class="op">[</span><span class="fl">1</span><span class="op">]</span>, max <span class="op">=</span> <span class="va">bb</span><span class="op">[</span><span class="fl">3</span><span class="op">]</span><span class="op">)</span>,</span>
<span> y <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/Uniform.html">runif</a></span><span class="op">(</span>n <span class="op">=</span> <span class="fl">10</span>, min <span class="op">=</span> <span class="va">bb</span><span class="op">[</span><span class="fl">2</span><span class="op">]</span>, max <span class="op">=</span> <span class="va">bb</span><span class="op">[</span><span class="fl">4</span><span class="op">]</span><span class="op">)</span></span>
<span><span class="op">)</span></span>
<span><span class="va">random_points</span> <span class="op">=</span> <span class="va">random_df</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>coords <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="st">"x"</span>, <span class="st">"y"</span><span class="op">)</span>, crs <span class="op">=</span> <span class="st">"EPSG:4326"</span><span class="op">)</span> <span class="co"># set coordinates and CRS</span></span></code></pre></div>
<p>The scenario illustrated in Figure <a href="spatial-operations.html#fig:spatial-join">4.6</a> shows that the <code>random_points</code> object (top left) lacks attribute data, while the <code>world</code> (top right) has attributes, including country names shown for a sample of countries in the legend.
Spatial joins are implemented with <code><a href="https://r-spatial.github.io/sf/reference/st_join.html">st_join()</a></code>, as illustrated in the code chunk below.
The output is the <code>random_joined</code> object which is illustrated in Figure <a href="spatial-operations.html#fig:spatial-join">4.6</a> (bottom left).
Before creating the joined dataset, we use spatial subsetting to create <code>world_random</code>, which contains only countries that contain random points, to verify the number of country names returned in the joined dataset should be four (Figure <a href="spatial-operations.html#fig:spatial-join">4.6</a>, top right panel).</p>
<div class="sourceCode" id="cb130"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">world_random</span> <span class="op">=</span> <span class="va">world</span><span class="op">[</span><span class="va">random_points</span>, <span class="op">]</span></span>
<span><span class="fu"><a href="https://rspatial.github.io/terra/reference/dimensions.html">nrow</a></span><span class="op">(</span><span class="va">world_random</span><span class="op">)</span></span>
<span><span class="co">#> [1] 4</span></span>
<span><span class="va">random_joined</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_join.html">st_join</a></span><span class="op">(</span><span class="va">random_points</span>, <span class="va">world</span><span class="op">[</span><span class="st">"name_long"</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:spatial-join"></span>
<img src="figures/spatial-join-1.png" alt="Illustration of a spatial join. A new attribute variable is added to random points (top left) from source world object (top right) resulting in the data represented in the final panel." width="100%"><p class="caption">
FIGURE 4.6: Illustration of a spatial join. A new attribute variable is added to random points (top left) from source world object (top right) resulting in the data represented in the final panel.
</p>
</div>
<p>By default, <code><a href="https://r-spatial.github.io/sf/reference/st_join.html">st_join()</a></code> performs a left join, meaning that the result is an object containing all rows from <code>x</code> including rows with no match in <code>y</code> (see Section <a href="attr.html#vector-attribute-joining">3.2.4</a>), but it can also do inner joins by setting the argument <code>left = FALSE</code>.
Like spatial subsetting, the default topological operator used by <code><a href="https://r-spatial.github.io/sf/reference/st_join.html">st_join()</a></code> is <code><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_intersects()</a></code>, which can be changed by setting the <code>join</code> argument (see <code><a href="https://r-spatial.github.io/sf/reference/st_join.html">?st_join</a></code> for details).
The example above demonstrates the addition of a column from a polygon layer to a point layer, but the approach works regardless of geometry types.
In such cases, for example when <code>x</code> contains polygons, each of which match multiple objects in <code>y</code>, spatial joins will result in duplicate features by creating a new row for each match in <code>y</code>.</p>
</div>
<div id="non-overlapping-joins" class="section level3" number="4.2.6">
<h3>
<span class="header-section-number">4.2.6</span> Distance-based joins<a class="anchor" aria-label="anchor" href="#non-overlapping-joins"><i class="fas fa-link"></i></a>
</h3>
<p>Sometimes two geographic datasets do not intersect but still have a strong geographic relationship due to their proximity.
The datasets <code>cycle_hire</code> and <code>cycle_hire_osm</code>, already attached in the <strong>spData</strong> package, provide a good example.
Plotting them shows that they are often closely related but they do not touch, as shown in Figure <a href="spatial-operations.html#fig:cycle-hire">4.7</a>, a base version of which is created with the following code below:
</p>
<div class="sourceCode" id="cb131"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rspatial.github.io/terra/reference/plot.html">plot</a></span><span class="op">(</span><span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_geometry.html">st_geometry</a></span><span class="op">(</span><span class="va">cycle_hire</span><span class="op">)</span>, col <span class="op">=</span> <span class="st">"blue"</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="fu"><a href="https://r-spatial.github.io/sf/reference/st_geometry.html">st_geometry</a></span><span class="op">(</span><span class="va">cycle_hire_osm</span><span class="op">)</span>, add <span class="op">=</span> <span class="cn">TRUE</span>, pch <span class="op">=</span> <span class="fl">3</span>, col <span class="op">=</span> <span class="st">"red"</span><span class="op">)</span></span></code></pre></div>
<p>We can check if any points are the same using <code><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_intersects()</a></code> as shown below:</p>
<div class="sourceCode" id="cb132"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/base/any.html">any</a></span><span class="op">(</span><span class="fu"><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_intersects</a></span><span class="op">(</span><span class="va">cycle_hire</span>, <span class="va">cycle_hire_osm</span>, sparse <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="co">#> [1] FALSE</span></span></code></pre></div>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:cycle-hire"></span>
<div class="leaflet html-widget html-fill-item" id="htmlwidget-9580c84d0c32f0a22a8c" style="width:100%;height:389.34px;"></div>
<script type="application/json" data-for="htmlwidget-9580c84d0c32f0a22a8c">{"x":{"options":{"crs":{"crsClass":"L.CRS.EPSG3857","code":null,"proj4def":null,"projectedBounds":null,"options":{}}},"calls":[{"method":"addCircles","args":[[51.52916347,51.49960695,51.52128377,51.53005939,51.49313,51.51811784,51.53430039,51.52834133,51.5073853,51.50597426,51.52395143,51.52168078,51.51991453,51.52994371,51.51772703,51.52635795,51.5216612,51.51477076,51.52505093,51.52773634,51.53007835,51.5222641,51.51943538,51.51908011,51.5288338,51.52728093,51.51382102,51.52351808,51.513735,51.52915444,51.52953709,51.52469624,51.5341235,51.50173726,51.49159394,51.4973875,51.5263778,51.52127071,51.52001715,51.53099181,51.52026,51.51073687,51.522511,51.507131,51.52334476,51.51248445,51.50706909,51.52867339,51.52671796,51.52295439,51.5099923,51.52174785,51.51707521,51.52058381,51.52334672,51.52452699,51.53176825,51.52644828,51.49738251,51.4907579,51.53089041,51.50946212,51.52522753,51.51795029,51.51882555,51.52059681,51.5262363,51.53136059,51.5154186,51.52352001,51.52572618,51.48591714,51.53219984,51.52559505,51.52341837,51.52486887,51.50069361,51.52025302,51.514274,51.50963938,51.51593725,51.50064702,51.48947903,51.51646835,51.51858757,51.5262503,51.53301907,51.49368637,51.49889832,51.53440868,51.49506109,51.5208417,51.53095071,51.49792478,51.52554222,51.51457763,51.49043573,51.51155322,51.51340693,51.50472376,51.51159481,51.51552971,51.51410514,51.52600832,51.49812559,51.51563144,51.53304322,51.5100172,51.51580998,51.49646288,51.52451738,51.51423368,51.51449962,51.49288067,51.49582705,51.52589324,51.51573534,51.51891348,51.52111369,51.52836014,51.49654462,51.50069491,51.51782144,51.51700801,51.49536226,51.5118973,51.50950627,51.53300545,51.52364804,51.50136494,51.504904,51.52326004,51.51196176,51.49087475,51.51227622,51.49488108,51.52096262,51.51530805,51.49372451,51.4968865,51.48894022,51.50074359,51.48836528,51.51494305,51.49211134,51.48478899,51.49705603,51.51213691,51.49217002,51.51183419,51.50379168,51.49586666,51.49443626,51.50039792,51.49085368,51.51461995,51.50663341,51.49234577,51.51760685,51.4931848,51.515607,51.517932,51.50185512,51.49395092,51.50040123,51.51474612,51.52784273,51.4916156,51.49121192,51.50486,51.512529,51.52174384,51.51791921,51.49616092,51.48985626,51.5129118,51.49941247,51.52202903,51.52071513,51.48805753,51.51733558,51.49247977,51.5165179,51.53166681,51.48997562,51.50311799,51.51251523,51.50581776,51.50462759,51.50724437,51.50368837,51.50556905,51.51048489,51.51492456,51.5225965,51.51236389,51.51821864,51.52659961,51.518144,51.518154,51.50135267,51.52505151,51.49358391,51.51681444,51.49464523,51.50658458,51.50274025,51.52683806,51.51906932,51.49094565,51.51615461,51.48971651,51.49016361,51.49066456,51.49481649,51.50275704,51.49148474,51.51476963,51.50935342,51.50844614,51.52891573,51.50742485,51.50654321,51.50669284,51.49396755,51.51501025,51.50777049,51.53450449,51.49571828,51.51838043,51.5084448,51.5233534,51.52443845,51.50545935,51.52200801,51.49096258,51.51611887,51.4853572,51.52285301,51.530052,51.50645179,51.49859784,51.489932,51.518908,51.5046364,51.53404294,51.53051587,51.52557531,51.51362054,51.52248185,51.50143293,51.50494561,51.5136846,51.5134891,51.49874469,51.51422502,51.52644342,51.51595344,51.50102668,51.49206037,51.49320445,51.50082346,51.4863434,51.53583617,51.50144456,51.50613324,51.49671237,51.52004497,51.50930161,51.50315739,51.5034938,51.51862243,51.490083,51.49580589,51.51906446,51.49914063,51.5272947,51.52367314,51.49236962,51.50923022,51.51678023,51.48483991,51.49888404,51.49815779,51.488226,51.50029631,51.49363156,51.49612799,51.50227992,51.49398524,51.50501351,51.51475963,51.466907,51.50295379,51.51211869,51.48677988,51.51816295,51.50990837,51.49907558,51.50963123,51.4908679,51.51918144,51.53088935,51.51734403,51.50088934,51.52536703,51.49768448,51.49679128,51.51004801,51.52702563,51.49357351,51.49785559,51.526293,51.523196,51.49652013,51.50908747,51.53266186,51.53114,51.53095,51.53589283,51.51641749,51.52085887,51.49334336,51.50860544,51.505044,51.51196803,51.504942,51.51108452,51.51175646,51.5333196,51.51444134,51.50810309,51.53692216,51.528246,51.48802358,51.50013942,51.51217033,51.51196,51.5017154373867,51.529423,51.486965,51.49459148,51.506767,51.486575,51.494412,51.520994,51.51643491,51.49782999,51.49675303,51.50391972,51.536264,51.5291212008901,51.519656,51.530344,51.5109192966489,51.5173721,51.50024195,51.520205,51.49775,51.515208,51.49980661,51.50402793,51.50194596,51.49188409,51.50535447,51.49559291,51.51431171,51.50686435,51.51953043,51.50935171,51.51310333,51.496481,51.51352755,51.49369988,51.51070161,51.51013066,51.521776,51.51066202,51.49942855,51.51733427,51.524826,51.492462,51.51809,51.51348,51.502319,51.52261762,51.517703,51.528187,51.52289229,51.51689296,51.50204238,51.49418566,51.512303,51.51222,51.51994326,51.493146,51.519968,51.49337264,51.488105,51.485821,51.504043,51.504044,51.49456127,51.490491,51.533379,51.493072,51.51397065,51.499917,51.48902,51.534474,51.496957,51.524564,51.488852,51.51824,51.488124,51.5338,51.483145,51.495656,51.510101,51.515256,51.52568,51.52388,51.528936,51.493381,51.50623,51.511088,51.515975,51.504719,51.505697,51.508447,51.487679,51.49447,51.538071,51.542138,51.504749,51.535179,51.516196,51.516,51.541603,51.534776,51.51417,51.53558,51.534464,51.523538,51.521564,51.513757,51.530535,51.533283,51.529452,51.496137,51.498125,51.499041,51.489096,51.49109,51.521889,51.527152,51.5128711,51.487129,51.509843,51.51328,51.528828,51.531127,51.525645,51.511811,51.506946,51.513074,51.508622,51.522507,51.524677,51.50196,51.528169,51.52512,51.527058,51.519265,51.522561,51.502635,51.520893,51.496454,51.520398,51.517475,51.528692,51.519362,51.517842,51.509303,51.511066,51.532091,51.5142228,51.516204,51.500088,51.5112,51.532513,51.528224,51.518811,51.526041,51.534137,51.525941,51.51549,51.511654,51.504714,51.503143,51.503802,51.507326,51.486892,51.508896,51.530326,51.50357,51.528222,51.531864,51.537349,51.51793,51.508981,51.531091,51.528302,51.506613,51.514115,51.509591,51.526153,51.539957,51.517428,51.509474,51.498386,51.49605,51.521564,51.503447,51.511542,51.535678,51.513548,51.502661,51.51601,51.493978,51.501391,51.511624,51.518369,51.51746,51.499286,51.509943,51.521905,51.509158,51.51616,51.53213,51.503083,51.506256,51.539099,51.485587,51.53356,51.497304,51.528869,51.527607,51.511246,51.48692917,51.497622,51.51244,51.490645,51.50964,51.52959,51.487196,51.53256,51.506093,51.5171,51.531066,51.513875,51.493267,51.472817,51.473471,51.494499,51.485743,51.481747,51.514767,51.472993,51.477839,51.538792,51.520331,51.504199,51.468814,51.491093,51.46512358,51.48256792,51.50646524,51.46925984,51.50403821,51.47817208,51.4768851,51.48102131,51.47107905,51.47625965,51.4737636,51.47518024,51.46086446,51.50748124,51.46193072,51.4729184,51.47761941,51.48438657,51.4687905,51.46881971,51.45995384,51.47073264,51.4795017,51.47053858,51.48959104,51.49434708,51.48606206,51.46706414,51.45787019,51.46663393,51.48357068,51.47286577,51.49824168,51.46916161,51.5190427,51.48373225,51.50173215,51.49886563,51.50035306,51.46904022,51.48180515,51.51687069,51.48089844,51.51148696,51.47084722,51.48267821,51.48294452,51.46841875,51.4996806,51.47729232,51.46437067,51.49087074,51.51632095,51.48498496,51.51323001,51.474376,51.46108367,51.49610093,51.5015946,51.49422354,51.47614939,51.46718562,51.475089,51.4646884,51.46517078,51.53546778,51.46348914,51.47787084,51.46866929,51.46231278,51.47453545,51.47768469,51.47303687,51.48810829,51.46506424,51.45475251,51.46067005,51.48814438,51.49760804,51.46095151,51.45922541,51.47047503,51.47946386,51.54211855,51.464786,51.53638435,51.48728535,51.53639219,51.53908372,51.5366541,51.47696496,51.47287627,51.52868155,51.51505991,51.45682071,51.45971528,51.50215353,51.49021762,51.46161068,51.46321128,51.47439218,51.48335692,51.51854104,51.54100708,51.47311696,51.51563007,51.51542791,51.53658514,51.53571683,51.53642464,51.48724429,51.53603947,51.52458353,51.46822047,51.45799126,51.47732253,51.47505096,51.47569809,51.46199911,51.47893931,51.49208492,51.47727637,51.50630441,51.50542628,51.46230566,51.46489445,51.47993289,51.47514228,51.48176572,51.51092871,51.5129814,51.51510818,51.46079243,51.46745485,51.47816972,51.4795738,51.48796408,51.53727795,51.53932857,51.4710956,51.51612862,51.45816465,51.49263658,51.5129006,51.48512191,51.47515398,51.48795853,51.51787005,51.52456169,51.52526975,51.48321729,51.50070305,51.52059714,51.46239255,51.46760141,51.45752945,51.45705988,51.46134382,51.473611,51.491026,51.509224,51.47250956,51.511891,51.470131,51.496664,51.460333,51.4619230679],[-0.109970527,-0.197574246,-0.084605692,-0.120973687,-0.156876,-0.144228881,-0.1680743,-0.170134484,-0.09644075100000001,-0.092754157,-0.122502346,-0.130431727,-0.136039674,-0.123616824,-0.127854211,-0.125979294,-0.109006325,-0.12221963,-0.131161087,-0.135273468,-0.13884627,-0.114079481,-0.119123345,-0.124678402,-0.132250369,-0.11829517,-0.107927706,-0.143613641,-0.193487,-0.093421615,-0.08335332300000001,-0.084439283,-0.129386874,-0.184980612,-0.192369256,-0.197245586,-0.07813092100000001,-0.0755789,-0.08391116799999999,-0.093903825,-0.157183945,-0.144165239,-0.162298,-0.06691,-0.183846408,-0.099141408,-0.145904427,-0.08745937600000001,-0.104298194,-0.094934859,-0.143495266,-0.09447507199999999,-0.086685542,-0.154701411,-0.120202614,-0.079248081,-0.114329032,-0.172190727,-0.089446947,-0.106323685,-0.089782579,-0.124749274,-0.13518856,-0.108657431,-0.108028472,-0.116688468,-0.134407652,-0.117069978,-0.098850915,-0.108340165,-0.08848618799999999,-0.124469948,-0.105480698,-0.144083893,-0.124121774,-0.099489485,-0.102091246,-0.141327271,-0.111257,-0.131510949,-0.111778348,-0.078600401,-0.115156562,-0.079684557,-0.132053392,-0.123509611,-0.139174593,-0.111014912,-0.100440521,-0.109025404,-0.08581448899999999,-0.09734016199999999,-0.078505384,-0.183834706,-0.138231303,-0.158264483,-0.122806861,-0.0929401,-0.076793375,-0.192538767,-0.07712132200000001,-0.190240716,-0.147301667,-0.096317627,-0.132102166,-0.132328837,-0.172528678,-0.157275636,-0.105270275,-0.183289032,-0.158963647,-0.07353765399999999,-0.141423695,-0.114934001,-0.13547809,-0.090847761,-0.093080779,-0.156166631,-0.078869751,-0.104724625,-0.150905245,-0.094524319,-0.096496865,-0.09388536,-0.185296516,-0.137043852,-0.07545948199999999,-0.136792671,-0.074754872,-0.191462381,-0.06797,-0.104708922,-0.097441687,-0.153319609,-0.157436972,-0.117974901,-0.085634242,-0.147203711,-0.198286569,-0.161203828,-0.111435796,-0.202759212,-0.129361842,-0.11614642,-0.138364847,-0.110683213,-0.168917077,-0.201554966,-0.101536865,-0.174292825,-0.11282408,-0.191933711,-0.092921165,-0.193068385,-0.196170309,-0.137841333,-0.131773845,-0.141334487,-0.121328408,-0.167894973,-0.183118788,-0.183716959,-0.159237081,-0.147624377,-0.195455928,-0.165164288,-0.108068155,-0.186753859,-0.173715911,-0.113001,-0.115163,-0.08111889999999999,-0.188098863,-0.140947636,-0.141923621,-0.153645496,-0.152317537,-0.165842551,-0.14521173,-0.140741432,-0.175810943,-0.178433004,-0.164393768,-0.109914711,-0.132845681,-0.153520935,-0.133201961,-0.100186337,-0.091773776,-0.106237501,-0.098497684,-0.111606696,-0.082989638,-0.06607803700000001,-0.161113413,-0.06954201,-0.100791005,-0.112432615,-0.06275,-0.062697,-0.153194766,-0.166304359,-0.165101392,-0.151926305,-0.158105512,-0.199004026,-0.149569201,-0.130504336,-0.088285377,-0.181190899,-0.08242239899999999,-0.170194408,-0.19039362,-0.166485083,-0.13045856,-0.155349725,-0.090220911,-0.188129731,-0.196422,-0.131961389,-0.115480888,-0.134621209,-0.123179697,-0.103137426,-0.17873226,-0.112753217,-0.130699733,-0.106992706,-0.110889274,-0.073438925,-0.067176443,-0.175116099,-0.138019439,-0.10569204,-0.151359288,-0.139625122,-0.128585022,-0.142207481,-0.099994052,-0.168314,-0.170279555,-0.096191134,-0.162727,-0.079249,-0.116542278,-0.08637971699999999,-0.106408455,-0.179592915,-0.116764211,-0.154907218,-0.178656971,-0.123247648,-0.135580879,-0.191351186,-0.103132904,-0.08066008299999999,-0.109256828,-0.169249375,-0.180246101,-0.132224622,-0.144132875,-0.089740764,-0.122492418,-0.156285395,-0.110699309,-0.114686385,-0.20528437,-0.09217644699999999,-0.084985356,-0.191496313,-0.07962099,-0.176645823,-0.162418,-0.127575233,-0.059642081,-0.112031483,-0.174653609,-0.128377673,-0.147478734,-0.151296092,-0.175488803,-0.138089062,-0.165471605,-0.209494128,-0.135635511,-0.092762704,-0.190603326,-0.106000855,-0.074189225,-0.136928582,-0.172729559,-0.148105415,-0.216573,-0.158456089,-0.16209757,-0.115853961,-0.135025698,-0.187842717,-0.08566631600000001,-0.119047563,-0.116911864,-0.140485596,-0.176770502,-0.138072691,-0.08315935200000001,-0.153463612,-0.141943703,-0.093913472,-0.138846453,-0.08854277100000001,-0.139956043,-0.081608045,-0.07395500000000001,-0.083067,-0.101384068,-0.129697889,-0.099981142,-0.086016,-0.085603,-0.160854428,-0.179135079,-0.089887855,-0.194757949,-0.193764092,-0.115851,-0.120718759,-0.115533,-0.197524944,-0.119643424,-0.111781191,-0.087587447,-0.12602103,-0.150181444,-0.10102611,-0.166878535,-0.113936001,-0.150481272,-0.142783033,-0.1798541843891,-0.097122,-0.116625,-0.134234258,-0.123702,-0.117286,-0.173881,-0.139016,-0.124332175,-0.135440826,-0.138733562,-0.11342628,-0.133952,-0.171185284853,-0.132339,-0.100168,-0.1511263847351,-0.1642075,-0.15934065,-0.174593,-0.10988,-0.117863,-0.176415994,-0.11386435,-0.194392952,-0.125674815,-0.113656543,-0.179077626,-0.200838199,-0.150666888,-0.13577731,-0.14744969,-0.13121385,-0.192404,-0.130110822,-0.121394101,-0.121723604,-0.155757901,-0.068856,-0.142345694,-0.179702476,-0.103604248,-0.176268,-0.159919,-0.163609,-0.17977,-0.200742,-0.071653961,-0.154106,-0.075375,-0.171681991,-0.158249929,-0.184400221,-0.18267094,-0.159988,-0.160785,-0.170704337,-0.099828,-0.169774,-0.09968067,-0.110121,-0.149004,-0.105312,-0.104778,-0.15393398,-0.149186,-0.139159,-0.129925,-0.09294031,-0.174554,-0.17524,-0.122203,-0.173894,-0.116279,-0.105593,-0.11655,-0.120903,-0.118677,-0.113134,-0.114605,-0.211358,-0.058641,-0.055312,-0.06507599999999999,-0.055894,-0.007542,-0.02296,-0.057159,-0.053177,-0.063531,-0.07054199999999999,-0.055167,-0.021582,-0.014409,-0.144664,-0.145393,-0.057544,-0.03338,-0.029138,-0.038775,-0.138853,-0.071881,-0.052099,-0.08248999999999999,-0.07634100000000001,-0.030556,-0.022694,-0.020467,-0.025492,-0.028155,-0.027616,-0.019355,-0.011457,-0.020157,-0.009205,-0.018716,-0.04667,-0.058005,-0.0389866,-0.009001,-0.02377,-0.047784,-0.013258,-0.048017,-0.06954299999999999,-0.025626,-0.058681,-0.064094,-0.06500599999999999,-0.041378,-0.03562,-0.016251,-0.018703,-0.015578,-0.025296,-0.021345,-0.054883,-0.022702,-0.051394,-0.009506000000000001,-0.026768,-0.07585500000000001,-0.059091,-0.074431,-0.090075,-0.025996,-0.053558,-0.06142,-0.055656,-0.155525,-0.211316,-0.014438,-0.033085,-0.037471,-0.011662,-0.047218,-0.037366,-0.036017,-0.013475,-0.179668,-0.014293,-0.008428,-0.215808,-0.145827,-0.170983,-0.012413,-0.042744,-0.020068,-0.069743,-0.066035,-0.147154,-0.067937,-0.00699,-0.075901,-0.144466,-0.142844,-0.033828,-0.204666,-0.102208,-0.145246,-0.107987,-0.002275,-0.107913,-0.104193,-0.039264,-0.016233,-0.056667,-0.062546,-0.005659,-0.021596,-0.0985,-0.127554,-0.205991,-0.205921,-0.043371,-0.12335,-0.009152,-0.117619,-0.063386,-0.224103,-0.18697,-0.08298999999999999,-0.017676,-0.218337,-0.141728,-0.18119,-0.09315,-0.022793,-0.047548,-0.057133,-0.09305099999999999,-0.102996299,-0.125978,-0.19096,-0.014582,-0.08497,-0.0801,-0.179369,-0.16862,-0.2242237,-0.18377,-0.11934,-0.117774,-0.21985,-0.199783,-0.20782,-0.228188,-0.223616,-0.124642,-0.225787,-0.133972,-0.116493,-0.138535,-0.163667,-0.210941,-0.210279,-0.216493,-0.157788279,-0.172078187,-0.208486599,-0.141812513,-0.217400093,-0.144690541,-0.215895601,-0.209973497,-0.207842908,-0.193254007,-0.197010096,-0.167160736,-0.187427294,-0.205535908,-0.180791784,-0.132102704,-0.149551631,-0.20481514,-0.158230901,-0.184318843,-0.190184054,-0.126994068,-0.141770709,-0.163041605,-0.209378594,-0.215804559,-0.214428378,-0.193502076,-0.174691623,-0.169821175,-0.202038682,-0.148059277,-0.117495865,-0.174485792,-0.204764421,-0.223852256,-0.100292412,-0.137424571,-0.217515071,-0.19627483,-0.18027465,-0.213872396,-0.183853573,-0.218190203,-0.17070367,-0.117661574,-0.219346128,-0.199135704,-0.221791552,-0.16478637,-0.174619404,-0.206029743,-0.202608612,-0.167919869,-0.211593602,-0.155442787,-0.191722864,-0.208158259,-0.222293381,-0.236769936,-0.1232585,-0.152248582,-0.201968,-0.173656546,-0.18038939,-0.11619105,-0.182126248,-0.126874471,-0.146544642,-0.211468596,-0.170210533,-0.170329317,-0.214749808,-0.22660621,-0.163750945,-0.195197203,-0.198735357,-0.222456468,-0.21145598,-0.20066766,-0.180884959,-0.152130083,-0.195777222,-0.028941601,-0.215618902,-0.102757578,-0.217995921,-0.112721065,-0.070329419,-0.07023031,-0.174347066,-0.176267008,-0.06555032099999999,-0.10534448,-0.202802098,-0.212145939,-0.083632928,-0.215087092,-0.21614583,-0.215550761,-0.163347594,-0.216305546,-0.034903714,-0.14326094,-0.137235175,-0.049067243,-0.02356501,-0.07588568599999999,-0.060291813,-0.054162264,-0.205279052,-0.026262677,-0.058631453,-0.190346493,-0.184806157,-0.138748723,-0.150908371,-0.20587627,-0.206240805,-0.208485293,-0.229116862,-0.189210466,-0.087262995,-0.150817316,-0.175407201,-0.17302926,-0.19411695,-0.187278987,-0.185273723,-0.214594781,-0.219486603,-0.208565479,-0.212607684,-0.172293499,-0.18243547,-0.17903854,-0.161765173,-0.079201849,-0.07428467499999999,-0.157850096,-0.120909408,-0.20600248,-0.234094148,-0.214762686,-0.174971902,-0.159169801,-0.187404506,-0.201005397,-0.165668686,-0.163795009,-0.211860644,-0.129698963,-0.032566533,-0.16829214,-0.20682737,-0.192165613,-0.200806304,-0.159322467,-0.191803,-0.209121,-0.216016,-0.122831913,-0.107349,-0.20464,-0.223868,-0.167029,-0.165297856693],10,null,null,{"interactive":true,"className":"","stroke":true,"color":"#03F","weight":5,"opacity":0.5,"fill":true,"fillColor":"#03F","fillOpacity":0.2},null,null,null,{"interactive":false,"permanent":false,"direction":"auto","opacity":1,"offset":[0,0],"textsize":"10px","textOnly":false,"className":"","sticky":true},null,null]},{"method":"addCircles","args":[[51.52912521362305,51.53401565551758,51.52729034423828,51.52582931518555,51.53001403808594,51.52594375610352,51.52664947509766,51.52085876464844,51.51682281494141,51.51473617553711,51.49416351318359,51.50061798095703,51.51780319213867,51.51329803466797,51.51405334472656,51.51224899291992,51.4882698059082,51.50081253051758,51.52448272705078,51.50056076049805,51.50337982177734,51.49507141113281,51.50747680664062,51.51416015625,51.51640701293945,51.51606750488281,51.50171279907227,51.50201034545898,51.49529647827148,51.51899337768555,51.51996994018555,51.51194000244141,51.52170944213867,51.51566314697266,51.51248550415039,51.497802734375,51.5099983215332,51.49079895019531,51.49155807495117,51.51231002807617,51.49150466918945,51.51769638061523,51.51214218139648,51.4901008605957,51.51478576660156,51.49358367919922,51.49903106689453,51.49145126342773,51.50217819213867,51.5250129699707,51.52245712280273,51.52287673950195,51.52532196044922,51.52196884155273,51.5164794921875,51.5145263671875,51.51681137084961,51.53588104248047,51.51811218261719,51.52334213256836,51.52565383911133,51.51673126220703,51.52835464477539,51.52250289916992,51.53087615966797,51.51727676391602,51.51192092895508,51.53435897827148,51.51531219482422,51.52731704711914,51.52647399902344,51.52053451538086,51.53308486938477,51.52260208129883,51.53691482543945,51.51469421386719,51.52555847167969,51.52023696899414,51.51406478881836,51.52341842651367,51.5207405090332,51.52516174316406,51.51856994628906,51.52337646484375,51.50148010253906,51.52470016479492,51.49089813232422,51.51391220092773,51.51790237426758,51.5232048034668,51.51870346069336,51.51491165161133,51.52326965332031,51.52390670776367,51.53091430664062,51.53261947631836,51.52450561523438,51.50960540771484,51.5159797668457,51.50600814819336,51.50090026855469,51.50944519042969,51.51948928833008,51.51821517944336,51.53089904785156,51.52178573608398,51.52351760864258,51.5269660949707,51.50862121582031,51.50478363037109,51.50655364990234,51.50936508178711,51.51110458374023,51.50312423706055,51.50276184082031,51.5063591003418,51.50270462036133,51.50302886962891,51.52287292480469,51.50719833374023,51.50673675537109,51.53334808349609,51.53291702270508,51.51734924316406,51.52548599243164,51.52022933959961,51.5236930847168,51.52767562866211,51.53001022338867,51.52880477905273,51.52627182006836,51.50925064086914,51.51859283447266,51.53307342529297,51.5191764831543,51.52675247192383,51.51288986206055,51.51224899291992,51.51071548461914,51.5118522644043,51.50993728637695,51.50992584228516,51.51190948486328,51.5057373046875,51.50535202026367,51.51589202880859,51.51157760620117,51.51702117919922,51.5191650390625,51.52096939086914,51.52167510986328,51.53197479248047,51.49435043334961,51.52944183349609,51.51446151733398,51.51610946655273,51.51904678344727,51.52439880371094,51.52506637573242,51.5257682800293,51.52671051025391,51.52104187011719,51.49912261962891,51.49487686157227,51.4920654296875,51.49595642089844,51.49212265014648,51.53086471557617,51.50924682617188,51.52347564697266,51.51543426513672,51.50062561035156,51.49895477294922,51.50366973876953,51.49112319946289,51.52059555053711,51.49082183837891,51.49552536010742,51.50117111206055,51.49394226074219,51.51353073120117,51.49674224853516,51.49799728393555,51.50186538696289,51.49391937255859,51.49413681030273,51.49236297607422,51.507080078125,51.49700927734375,51.49364471435547,51.49655914306641,51.49946212768555,51.49881744384766,51.50128936767578,51.49464797973633,51.49319076538086,51.49240493774414,51.49787902832031,51.49010848999023,51.48990249633789,51.49111938476562,51.48974227905273,51.48804092407227,51.49075317382812,51.50144577026367,51.51197814941406,51.51037216186523,51.49218368530273,51.51503753662109,51.50650405883789,51.51491928100586,51.5050163269043,51.51174926757812,51.51435470581055,51.506591796875,51.50767517089844,51.50740432739258,51.50800323486328,51.50844955444336,51.50907516479492,51.50968933105469,51.51758193969727,51.51374435424805,51.51470565795898,51.51447296142578,51.53133010864258,51.52838516235352,51.52645874023438,51.5267333984375,51.53044891357422,51.52826690673828,51.52888107299805,51.5291633605957,51.5278205871582,51.53181838989258,51.53451538085938,51.53160095214844,51.53448486328125,51.51836013793945,51.51002502441406,51.52993011474609,51.52622222900391,51.50460433959961,51.49643325805664,51.51578903198242,51.51375961303711,51.51592254638672,51.50947189331055,51.49486923217773,51.49581527709961,51.49396896362305,51.49613571166992,51.52635192871094,51.48826599121094,51.49322128295898,51.49234771728516,51.49362945556641,51.48994445800781,51.50579833984375,51.48537826538086,51.48810958862305,51.48983383178711,51.4981575012207,51.51644515991211,51.53055953979492,51.48485565185547,51.51791381835938,51.51488494873047,51.51557922363281,51.51791763305664,51.51560592651367,51.51712799072266,51.51890182495117,51.51425933837891,51.51215744018555,51.52226638793945,51.49868392944336,51.49961090087891,51.50031661987305,51.49734878540039,51.4966926574707,51.50032806396484,51.49629974365234,51.50531768798828,51.50063705444336,51.49807739257812,51.51834106445312,51.49370193481445,51.49331665039062,51.4958381652832,51.50140762329102,51.51897048950195,51.50468826293945,51.52130126953125,51.49744033813477,51.49075317382812,51.48680877685547,51.48625946044922,51.49372863769531,51.48944473266602,51.49286270141602,51.48483657836914,51.49043655395508,51.49571228027344,51.48894882202148,51.48603057861328,51.4908561706543,51.50009155273438,51.51220703125,51.50178909301758,51.51728439331055,51.50294494628906,51.51157760620117,51.51779556274414,51.51667022705078,51.49860000610352,51.49674987792969,51.51644897460938,51.51970672607422,51.49796676635742,51.50197982788086,51.52048873901367,51.49689483642578,51.51947021484375,51.52097702026367,51.51348495483398,51.51299285888672,51.52950668334961,51.51092147827148,51.51813125610352,51.50935363769531,51.51070404052734,51.50686645507812,51.49713134765625,51.4997673034668,51.50234985351562,51.53605651855469,51.5134391784668,51.51768112182617,51.51806259155273,51.51361083984375,51.51009750366211,51.51062393188477,51.49942779541016,51.49175262451172,51.51727676391602,51.5125617980957,51.5200080871582,51.50026702880859,51.52265167236328,51.51991271972656,51.52025604248047,51.51988983154297,51.51507186889648,51.50411224365234,51.52243041992188,51.52361297607422,51.5246696472168,51.52186584472656,51.50592422485352,51.49340438842773,51.49310684204102,51.53213119506836,51.53573226928711,51.51632308959961,51.53445053100586,51.53379821777344,51.49440765380859,51.52646636962891,51.53155136108398,51.53030395507812,51.53110504150391,51.51013946533203,51.5091667175293,51.50912857055664,51.51172637939453,51.51735305786133,51.50736999511719,51.50645065307617,51.50137329101562,51.50635147094727,51.53342437744141,51.51915740966797,51.5282096862793,51.50611114501953,51.50505065917969,51.50486373901367,51.50396728515625,51.50378799438477,51.51714706420898,51.51609420776367,51.51236724853516,51.52350616455078,51.51633834838867,51.51122665405273,51.51162338256836,51.51406478881836,51.51754379272461,51.52091598510742,51.52569198608398,51.52255630493164,51.53409194946289,51.52899169921875,51.53515625,51.53250503540039,51.53331756591797,51.53050994873047,51.5294303894043,51.52883529663086,51.52602005004883,51.49992370605469,51.52867889404297,51.51124572753906,51.51420593261719,51.51108169555664,51.51920318603516,51.5088996887207,51.5390510559082,51.50929641723633,51.48754119873047,51.5162467956543,51.53112411499023,51.53093719482422,51.50707626342773,51.49605560302734,51.52891540527344,51.52882385253906,51.52891540527344,51.52885055541992,51.52882766723633,51.52885818481445,51.52874755859375,51.52875518798828,51.49691772460938,51.50173568725586,51.45926666259766,51.46195983886719,51.46232986450195,51.4643669128418,51.46463394165039,51.46493530273438,51.51509857177734,51.49330902099609,51.49103164672852,51.54203796386719,51.53428268432617,51.53427886962891,51.53406143188477,51.53406143188477,51.4881706237793,51.51298141479492,51.53634262084961,51.49882507324219,51.50459289550781,51.52611923217773,51.5133056640625,51.5135383605957,51.50959014892578,51.51366806030273,51.52515029907227,51.5159912109375,51.46982192993164,51.46978759765625,51.46990585327148,51.46994018554688,51.48889923095703,51.52953720092773,51.51374435424805,51.52818298339844,51.50558471679688,51.50534057617188,51.5048713684082,51.50386810302734,51.51820755004883,51.51516342163086,51.52561187744141,51.52750015258789,51.52717208862305,51.49206924438477,51.47304916381836,51.51152038574219,51.52157974243164,51.51946258544922,51.53804016113281,51.49589157104492,51.50403213500977,51.50404739379883,51.50413131713867,51.52257537841797,51.50473785400391,51.50476455688477,51.50480651855469,51.50482940673828,51.52818298339844,51.47513198852539,51.53466033935547,51.52820587158203,51.54682540893555,51.54629898071289,51.50030136108398,51.53590393066406,51.52753829956055,51.540283203125,51.53598022460938,51.51189804077148,51.49393081665039,51.48147964477539,51.54103469848633,51.54169464111328,51.48582458496094,51.46757125854492,51.53876876831055,51.51645278930664,51.50397491455078,51.47622680664062,51.47941589355469,51.47994232177734,51.48358535766602,51.48390197753906,51.50931167602539,51.50931930541992,51.50932312011719,51.50933837890625,51.50934219360352,51.50983428955078,51.50983810424805,51.50997161865234,51.50997543334961,51.54087066650391,51.54308319091797,51.53842926025391,51.52247619628906,51.48233413696289,51.51785659790039,51.51787567138672,51.51789855957031,51.51802825927734,51.49266815185547,51.50943374633789,51.50895309448242],[-0.09338779747486115,-0.129309207201004,-0.1182352006435394,-0.09083600342273712,-0.1210571974515915,-0.1038272008299828,-0.1123251020908356,-0.08988530188798904,-0.1582169979810715,-0.1221619993448257,-0.1825312972068787,-0.09451589733362198,-0.09637469798326492,-0.07670540362596512,-0.07358500361442566,-0.06940989941358566,-0.1356287002563477,-0.0898251011967659,-0.1588800996541977,-0.07856839895248413,-0.07957950234413147,-0.08594369888305664,-0.09640560299158096,-0.08061859756708145,-0.07962480187416077,-0.08212079852819443,-0.1848890036344528,-0.1843793988227844,-0.1852709054946899,-0.1247294023633003,-0.1358640938997269,-0.1207367032766342,-0.1304273009300232,-0.1322115063667297,-0.1331615000963211,-0.08164320141077042,-0.157212495803833,-0.196125403046608,-0.186663806438446,-0.1597979962825775,-0.1924726963043213,-0.1280096024274826,-0.1620880961418152,-0.1904651969671249,-0.1652061939239502,-0.1906657963991165,-0.08559329807758331,-0.09012889862060547,-0.07422350347042084,-0.1662103980779648,-0.1548451036214828,-0.171658992767334,-0.153420701622963,-0.1513393968343735,-0.1644168943166733,-0.1582493036985397,-0.1518975049257278,-0.160705104470253,-0.1441397964954376,-0.1838506013154984,-0.1439844071865082,-0.1755494028329849,-0.1700911968946457,-0.1622716933488846,-0.1768050044775009,-0.1758725941181183,-0.174384206533432,-0.1681527942419052,-0.1470558047294617,-0.1746665984392166,-0.1721556931734085,-0.1547646969556808,-0.1726416945457458,-0.1610399037599564,-0.150157704949379,-0.1480903029441833,-0.1795047968626022,-0.1570902019739151,-0.1472848057746887,-0.175151601433754,-0.1450770050287247,-0.135173499584198,-0.1766694933176041,-0.1241701990365982,-0.1107200980186462,-0.08494079858064651,-0.1394933015108109,-0.09281720221042633,-0.1084283962845802,-0.1047310009598732,-0.108025997877121,-0.06618840247392654,-0.1204117983579636,-0.1224915981292725,-0.09387639909982681,-0.09995549917221069,-0.07924109697341919,-0.07465779781341553,-0.1692546010017395,-0.09271299839019775,-0.08347310125827789,-0.1243832036852837,-0.1191532015800476,-0.06269069761037827,-0.07854320108890533,-0.1092348992824554,-0.1083744987845421,-0.08861500024795532,-0.1937263011932373,-0.1925584971904755,-0.1989907026290894,-0.1963624060153961,-0.1974851042032242,-0.1535256057977676,-0.1494038999080658,-0.1701322048902512,-0.1552689969539642,-0.191404402256012,-0.09995950013399124,-0.1061898022890091,-0.1032897979021072,-0.1117890030145645,-0.1367686986923218,-0.1381060928106308,-0.1382132023572922,-0.1412868946790695,-0.1284389048814774,-0.1353798061609268,-0.1387432068586349,-0.1322190016508102,-0.1342364996671677,-0.151116207242012,-0.1319689005613327,-0.1391572952270508,-0.1405548006296158,-0.1305487006902695,-0.1536446958780289,-0.1575558930635452,-0.1441929936408997,-0.142809197306633,-0.1388102024793625,-0.1879072934389114,-0.1370330005884171,-0.09979409724473953,-0.105653703212738,-0.09304329752922058,-0.09291200339794159,-0.09384860098361969,-0.08823960274457932,-0.08566469699144363,-0.09442149847745895,-0.1054814979434013,-0.09292580187320709,-0.08337190002202988,-0.08765500038862228,-0.1285894960165024,-0.05971070006489754,-0.1381545066833496,-0.1311517059803009,-0.0884820967912674,-0.0781090036034584,-0.07875949889421463,-0.1120520979166031,-0.1179369017481804,-0.1321800053119659,-0.1353285014629364,-0.1383174955844879,-0.08993770182132721,-0.08472809940576553,-0.1435666978359222,-0.09880249947309494,-0.1019288972020149,-0.1001908034086227,-0.09848310053348541,-0.09712529927492142,-0.116697296500206,-0.1813008040189743,-0.1791877001523972,-0.1802363991737366,-0.1787330061197281,-0.1913225948810577,-0.1388126015663147,-0.1438243985176086,-0.1592289954423904,-0.1476075947284698,-0.1546085029840469,-0.1476241052150726,-0.1459303945302963,-0.1688434034585953,-0.1649394929409027,-0.1509203016757965,-0.152325302362442,-0.165447399020195,-0.1532440930604935,-0.1580667048692703,-0.1680227965116501,-0.1784097999334335,-0.1838334947824478,-0.1625753045082092,-0.1625996977090836,-0.1736893951892853,-0.1702129989862442,-0.1668481975793839,-0.1662604063749313,-0.1784280985593796,-0.09751179814338684,-0.08290299773216248,-0.1015529036521912,-0.1126210987567902,-0.1232742965221405,-0.1162168979644775,-0.1727204024791718,-0.1197232007980347,-0.1184678003191948,-0.1317393034696579,-0.131008505821228,-0.1346541047096252,-0.1259603053331375,-0.1318870931863785,-0.1297453045845032,-0.1314696967601776,-0.1213387995958328,-0.1354451030492783,-0.1375370025634766,-0.1414363980293274,-0.1170132011175156,-0.1010444983839989,-0.1091234982013702,-0.1042855009436607,-0.1063415035605431,-0.1047158986330032,-0.1154187023639679,-0.1099506989121437,-0.1080510020256042,-0.1142947971820831,-0.1090492978692055,-0.1098138988018036,-0.10696080327034,-0.0730184018611908,-0.1434559971094131,-0.1235487014055252,-0.1234614998102188,-0.09175200015306473,-0.1016196012496948,-0.1052607968449593,-0.1078827977180481,-0.1117516979575157,-0.1190444976091385,-0.1305709928274155,-0.1275409013032913,-0.1368478983640671,-0.1408818960189819,-0.1259883940219879,-0.1292545050382614,-0.1441729068756104,-0.1413037925958633,-0.1400548070669174,-0.132818803191185,-0.136564701795578,-0.1421748995780945,-0.1405968070030212,-0.1420411020517349,-0.1319441050291061,-0.1791722029447556,-0.1674931049346924,-0.138083204627037,-0.1879784017801285,-0.1881255954504013,-0.1831178069114685,-0.1836283057928085,-0.190236896276474,-0.08666019886732101,-0.1560654938220978,-0.2009900957345963,-0.201490193605423,-0.11407820135355,-0.1039230972528458,-0.1975031048059464,-0.1930058002471924,-0.1972882002592087,-0.2052401006221771,-0.195428803563118,-0.1058785989880562,-0.112193301320076,-0.2025800049304962,-0.2094330042600632,-0.1006769984960556,-0.198381707072258,-0.1947298943996429,-0.1918330043554306,-0.1915650963783264,-0.07889190316200256,-0.1232506036758423,-0.08454930037260056,-0.08951500058174133,-0.1062221974134445,-0.1158493012189865,-0.1222822964191437,-0.1110500991344452,-0.1150930970907211,-0.1149597987532616,-0.1106337010860443,-0.1227603033185005,-0.1110092028975487,-0.1114194020628929,-0.1243795976042747,-0.1168795973062515,-0.1139174029231071,-0.1504797041416168,-0.1797240972518921,-0.1644082069396973,-0.1585202068090439,-0.0769961029291153,-0.09002619981765747,-0.1795178949832916,-0.09621690213680267,-0.09396609663963318,-0.1242700964212418,-0.1323571056127548,-0.1350177973508835,-0.194334402680397,-0.09731210023164749,-0.1613789945840836,-0.1357111930847168,-0.1389185041189194,-0.1300491988658905,-0.1311278939247131,-0.09709309786558151,-0.1512179970741272,-0.1349456012248993,-0.1474936008453369,-0.1421763002872467,-0.1506551057100296,-0.1922765970230103,-0.1763094961643219,-0.2007306963205338,-0.133782297372818,-0.1797408014535904,-0.1540451943874359,-0.1635331958532333,-0.1934694051742554,-0.1556915044784546,-0.1216448023915291,-0.1797050982713699,-0.1256998926401138,-0.1035635992884636,-0.115013100206852,-0.09218680113554001,-0.09272850304841995,-0.07158850133419037,-0.1697838008403778,-0.1745657026767731,-0.1705185025930405,-0.1181337013840675,-0.2174990028142929,-0.04176019877195358,-0.07494830340147018,-0.03569810092449188,-0.04663209989666939,-0.07068169862031937,-0.09982819855213165,-0.09984319657087326,-0.06135139986872673,-0.06271539628505707,-0.09847539663314819,-0.1220870018005371,-0.118647001683712,-0.1735181957483292,-0.0286301001906395,-0.06618860363960266,-0.04279280081391335,-0.04798439890146255,-0.2114519029855728,-0.2242292016744614,-0.2196248024702072,-0.2059940993785858,-0.1231873035430908,-0.1457414031028748,-0.1430086046457291,-0.2059323042631149,-0.2184094041585922,-0.09320160001516342,-0.1477314978837967,-0.03735850006341934,-0.1146802976727486,-0.1159529015421867,-0.1153946965932846,-0.1138684973120689,-0.1130378022789955,-0.1834370046854019,-0.1873559057712555,-0.1909584999084473,-0.0305780004709959,-0.02917139977216721,-0.09299620240926743,-0.06898610293865204,-0.1111695989966393,-0.1079486981034279,-0.05132989957928658,-0.05524060130119324,-0.05480609834194183,-0.03731809929013252,-0.05588379874825478,-0.03340740129351616,-0.03296750038862228,-0.02815829962491989,-0.02544780075550079,-0.02760040014982224,-0.0477617010474205,-0.0358303003013134,-0.1745879054069519,-0.08740100264549255,-0.01404820010066032,-0.03362049907445908,-0.05736390128731728,-0.02137940004467964,-0.01240990031510592,-0.1416262984275818,-0.02587329968810081,-0.1168603971600533,-0.1552148014307022,-0.08595810085535049,-0.08555299788713455,-0.06689970195293427,-0.1042110025882721,-0.01334539987146854,-0.01331450045108795,-0.01330539956688881,-0.01330920029431581,-0.01334969978779554,-0.01334539987146854,-0.01334269996732473,-0.01337889954447746,-0.1738771051168442,-0.1003087982535362,-0.1808453053236008,-0.1808360069990158,-0.1753010004758835,-0.1747123003005981,-0.1738660931587219,-0.1728828996419907,-0.1054434031248093,-0.2191217988729477,-0.216540202498436,-0.02883020043373108,-0.08636900037527084,-0.08640450239181519,-0.08638259768486023,-0.08634699881076813,-0.2223017960786819,-0.06408549845218658,-0.1025628000497818,-0.1373721957206726,-0.1165795028209686,-0.04696319997310638,-0.04797070100903511,-0.1167192980647087,-0.08465509861707687,-0.1176247969269753,-0.015591099858284,-0.1208008974790573,-0.1408015936613083,-0.1407676041126251,-0.1404854953289032,-0.140521302819252,-0.105547197163105,-0.08009859919548035,-0.02040489949285984,-0.07548379898071289,-0.1117931008338928,-0.1136531010270119,-0.1129046976566315,-0.1134240031242371,-0.1165020987391472,-0.05844509974122047,-0.06955330073833466,-0.05706809833645821,-0.05791860073804855,-0.2291229963302612,-0.2147257030010223,-0.05670920014381409,-0.02237650007009506,-0.0744313970208168,-0.1446426063776016,-0.1728768050670624,-0.2174569964408875,-0.2173631936311722,-0.2174052000045776,-0.04101530089974403,-0.06754639744758606,-0.06752920150756836,-0.06778910011053085,-0.06777189671993256,-0.06979779899120331,-0.1592990010976791,-0.124966099858284,-0.06943450123071671,-0.01454740017652512,-0.01002360042184591,-0.1590677946805954,-0.1560537070035934,-0.1348813027143478,-0.02166059985756874,-0.02660050056874752,-0.1071562990546227,-0.1274670958518982,-0.1380670964717865,-0.1431670933961868,-0.1390734016895294,-0.1488102972507477,-0.2066828012466431,-0.1384492963552475,-0.1183784976601601,-0.01320760045200586,-0.1932822018861771,-0.1957414001226425,-0.1941318064928055,-0.2020470052957535,-0.1975594013929367,-0.02592330053448677,-0.02581189945340157,-0.02573350071907043,-0.02573220059275627,-0.02582200057804585,-0.02373870089650154,-0.0237748995423317,-0.02368880063295364,-0.02372509986162186,-0.01074429973959923,-0.007984300144016743,-0.01189500000327826,-0.04172470048069954,-0.1362718045711517,-0.04322149977087975,-0.04317649826407433,-0.04341059923171997,-0.04321610182523727,-0.0923537015914917,-0.002419099910184741,-0.006909300107508898],10,null,null,{"interactive":true,"className":"","stroke":true,"color":"red","weight":5,"opacity":0.5,"fill":true,"fillColor":"red","fillOpacity":0.2},null,null,null,{"interactive":false,"permanent":false,"direction":"auto","opacity":1,"offset":[0,0],"textsize":"10px","textOnly":false,"className":"","sticky":true},null,null]}],"limits":{"lat":[51.45475251,51.54682540893555],"lng":[-0.236769936,-0.002275]}},"evals":[],"jsHooks":[]}</script><p class="caption">
FIGURE 4.7: The spatial distribution of cycle hire points in London based on official data (blue) and OpenStreetMap data (red).
</p>
</div>
<p>Imagine that we need to join the <code>capacity</code> variable in <code>cycle_hire_osm</code> onto the official ‘target’ data contained in <code>cycle_hire</code>.
This is when a non-overlapping join is needed.
The simplest method is to use the binary predicate <code><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_is_within_distance()</a></code>, as demonstrated below using a threshold distance of 20 m.
One can set the threshold distance in metric units also for unprojected data (e.g., lon/lat CRSs such as WGS84), if the spherical geometry engine (S2) is enabled, as it is in <strong>sf</strong> by default (see Section <a href="spatial-class.html#s2">2.2.9</a>).</p>
<div class="sourceCode" id="cb133"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">sel</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/geos_binary_pred.html">st_is_within_distance</a></span><span class="op">(</span><span class="va">cycle_hire</span>, <span class="va">cycle_hire_osm</span>, </span>
<span> dist <span class="op">=</span> <span class="fu">units</span><span class="fu">::</span><span class="fu"><a href="https://r-quantities.github.io/units/reference/units.html">set_units</a></span><span class="op">(</span><span class="fl">20</span>, <span class="st">"m"</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rspatial.github.io/terra/reference/summary.html">summary</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/lengths.html">lengths</a></span><span class="op">(</span><span class="va">sel</span><span class="op">)</span> <span class="op">></span> <span class="fl">0</span><span class="op">)</span></span>
<span><span class="co">#> Mode FALSE TRUE </span></span>
<span><span class="co">#> logical 304 438</span></span></code></pre></div>
<p>This shows that there are 438 points in the target object <code>cycle_hire</code> within the threshold distance of <code>cycle_hire_osm</code>.
How to retrieve the <em>values</em> associated with the respective <code>cycle_hire_osm</code> points?
The solution is again with <code><a href="https://r-spatial.github.io/sf/reference/st_join.html">st_join()</a></code>, but with an additional <code>dist</code> argument (set to 20 m below):</p>
<div class="sourceCode" id="cb134"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">z</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_join.html">st_join</a></span><span class="op">(</span><span class="va">cycle_hire</span>, <span class="va">cycle_hire_osm</span>, <span class="va">st_is_within_distance</span>, </span>
<span> dist <span class="op">=</span> <span class="fu">units</span><span class="fu">::</span><span class="fu"><a href="https://r-quantities.github.io/units/reference/units.html">set_units</a></span><span class="op">(</span><span class="fl">20</span>, <span class="st">"m"</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rspatial.github.io/terra/reference/dimensions.html">nrow</a></span><span class="op">(</span><span class="va">cycle_hire</span><span class="op">)</span></span>
<span><span class="co">#> [1] 742</span></span>
<span><span class="fu"><a href="https://rspatial.github.io/terra/reference/dimensions.html">nrow</a></span><span class="op">(</span><span class="va">z</span><span class="op">)</span></span>
<span><span class="co">#> [1] 762</span></span></code></pre></div>
<p>Note that the number of rows in the joined result is greater than the target.
This is because some cycle hire stations in <code>cycle_hire</code> have multiple matches in <code>cycle_hire_osm</code>.
To aggregate the values for the overlapping points and return the mean, we can use the aggregation methods learned in Chapter <a href="attr.html#attr">3</a>, resulting in an object with the same number of rows as the target.</p>
<div class="sourceCode" id="cb135"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">z</span> <span class="op">=</span> <span class="va">z</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/summarise.html">summarize</a></span><span class="op">(</span>capacity <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/summarize-generics.html">mean</a></span><span class="op">(</span><span class="va">capacity</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="fu"><a href="https://rspatial.github.io/terra/reference/dimensions.html">nrow</a></span><span class="op">(</span><span class="va">z</span><span class="op">)</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">cycle_hire</span><span class="op">)</span></span>
<span><span class="co">#> [1] TRUE</span></span></code></pre></div>
<p>The capacity of nearby stations can be verified by comparing a plot of the capacity of the source <code>cycle_hire_osm</code> data with the results in this new object (plots not shown):</p>
<div class="sourceCode" id="cb136"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rspatial.github.io/terra/reference/plot.html">plot</a></span><span class="op">(</span><span class="va">cycle_hire_osm</span><span class="op">[</span><span class="st">"capacity"</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">z</span><span class="op">[</span><span class="st">"capacity"</span><span class="op">]</span><span class="op">)</span></span></code></pre></div>
<p>The result of this join has used a spatial operation to change the attribute data associated with simple features; the geometry associated with each feature has remained unchanged.</p>
</div>
<div id="spatial-aggr" class="section level3" number="4.2.7">
<h3>
<span class="header-section-number">4.2.7</span> Spatial aggregation<a class="anchor" aria-label="anchor" href="#spatial-aggr"><i class="fas fa-link"></i></a>
</h3>
<p>As with attribute data aggregation, spatial data aggregation <em>condenses</em> data: aggregated outputs have fewer rows than non-aggregated inputs.
Statistical <em>aggregating functions</em>, such as mean average or sum, summarise multiple values of a variable, and return a single value per <em>grouping variable</em>.
Section <a href="attr.html#vector-attribute-aggregation">3.2.3</a> demonstrated how <code><a href="https://rspatial.github.io/terra/reference/aggregate.html">aggregate()</a></code> and <code>group_by() |> summarize()</code> condense data based on attribute variables, this section shows how the same functions work with spatial objects.
</p>
<p>Returning to the example of New Zealand, imagine you want to find out the average height of high points in each region: it is the geometry of the source (<code>y</code> or <code>nz</code> in this case) that defines how values in the target object (<code>x</code> or <code>nz_height</code>) are grouped.
This can be done in a single line of code with base R’s <code><a href="https://rspatial.github.io/terra/reference/aggregate.html">aggregate()</a></code> method.</p>
<div class="sourceCode" id="cb137"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">nz_agg</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/aggregate.html">aggregate</a></span><span class="op">(</span>x <span class="op">=</span> <span class="va">nz_height</span>, by <span class="op">=</span> <span class="va">nz</span>, FUN <span class="op">=</span> <span class="va">mean</span><span class="op">)</span></span></code></pre></div>
<p>The result of the previous command is an <code>sf</code> object with the same geometry as the (spatial) aggregating object (<code>nz</code>), which you can verify with the command <code>identical(st_geometry(nz), st_geometry(nz_agg))</code>.
The result of the previous operation is illustrated in Figure <a href="spatial-operations.html#fig:spatial-aggregation">4.8</a>, which shows the average value of features in <code>nz_height</code> within each of New Zealand’s 16 regions.
The same result can also be generated by piping the output from <code><a href="https://r-spatial.github.io/sf/reference/st_join.html">st_join()</a></code> into the ‘tidy’ functions <code><a href="https://dplyr.tidyverse.org/reference/group_by.html">group_by()</a></code> and <code><a href="https://dplyr.tidyverse.org/reference/summarise.html">summarize()</a></code> as follows:</p>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:spatial-aggregation"></span>
<img src="figures/spatial-aggregation-1.png" alt="Average height of the top 101 high points across the regions of New Zealand." width="50%"><p class="caption">
FIGURE 4.8: Average height of the top 101 high points across the regions of New Zealand.
</p>
</div>
<div class="sourceCode" id="cb138"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">nz_agg2</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/st_join.html">st_join</a></span><span class="op">(</span>x <span class="op">=</span> <span class="va">nz</span>, y <span class="op">=</span> <span class="va">nz_height</span><span class="op">)</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">Name</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>elevation <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/summarize-generics.html">mean</a></span><span class="op">(</span><span class="va">elevation</span>, na.rm <span class="op">=</span> <span class="cn">TRUE</span><span class="op">)</span><span class="op">)</span></span></code></pre></div>
<p>The resulting <code>nz_agg</code> objects have the same geometry as the aggregating object <code>nz</code> but with a new column summarizing the values of <code>x</code> in each region using the function <code><a href="https://rspatial.github.io/terra/reference/summarize-generics.html">mean()</a></code>.
Other functions could be used instead of <code><a href="https://rspatial.github.io/terra/reference/summarize-generics.html">mean()</a></code> here, including <code><a href="https://rspatial.github.io/terra/reference/summarize-generics.html">median()</a></code>, <code><a href="https://rdrr.io/r/stats/sd.html">sd()</a></code> and other functions that return a single value per group.
Note: one difference between the <code><a href="https://rspatial.github.io/terra/reference/aggregate.html">aggregate()</a></code> and <code>group_by() |> summarize()</code> approaches is that the former results in <code>NA</code> values for unmatching region names while the latter preserves region names.
The ‘tidy’ approach is thus more flexible in terms of aggregating functions and the column names of the results.
Aggregating operations that also create new geometries are covered in Section <a href="geometry-operations.html#geometry-unions">5.2.7</a>.</p>
</div>
<div id="incongruent" class="section level3" number="4.2.8">
<h3>
<span class="header-section-number">4.2.8</span> Joining incongruent layers<a class="anchor" aria-label="anchor" href="#incongruent"><i class="fas fa-link"></i></a>
</h3>
<p>Spatial congruence is an important concept related to spatial aggregation.
An <em>aggregating object</em> (which we will refer to as <code>y</code>) is <em>congruent</em> with the target object (<code>x</code>) if the two objects have shared borders.
Often this is the case for administrative boundary data, whereby larger units — such as Middle Layer Super Output Areas (<a href="https://www.ons.gov.uk/methodology/geography/ukgeographies/censusgeography">MSOAs</a>) in the UK or districts in many other European countries — are composed of many smaller units.</p>
<p><em>Incongruent</em> aggregating objects, by contrast, do not share common borders with the target <span class="citation">(<a href="references.html#ref-qiu_development_2012">Qiu, Zhang, and Zhou 2012</a>)</span>.
This is problematic for spatial aggregation (and other spatial operations) illustrated in Figure <a href="spatial-operations.html#fig:areal-example">4.9</a>: aggregating the centroid of each sub-zone will not return accurate results.
Areal interpolation overcomes this issue by transferring values from one set of areal units to another, using a range of algorithms including simple area weighted approaches and more sophisticated approaches such as ‘pycnophylactic’ methods <span class="citation">(<a href="references.html#ref-tobler_smooth_1979">Waldo R. Tobler 1979</a>)</span>.</p>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:areal-example"></span>
<img src="figures/areal-example-1.png" alt="Illustration of congruent (left) and incongruent (right) areal units with respect to larger aggregating zones (translucent red borders)." width="100%"><p class="caption">
FIGURE 4.9: Illustration of congruent (left) and incongruent (right) areal units with respect to larger aggregating zones (translucent red borders).
</p>
</div>
<p>The <strong>spData</strong> package contains a dataset named <code>incongruent</code> (colored polygons with black borders in the right panel of Figure <a href="spatial-operations.html#fig:areal-example">4.9</a>) and a dataset named <code>aggregating_zones</code> (the two polygons with the translucent blue border in the right panel of Figure <a href="spatial-operations.html#fig:areal-example">4.9</a>).
Let us assume that the <code>value</code> column of <code>incongruent</code> refers to the total regional income in million Euros.
How can we transfer the values of the underlying nine spatial polygons into the two polygons of <code>aggregating_zones</code>?</p>
<p>The simplest useful method for this is <em>area weighted</em> spatial interpolation, which transfers values from the <code>incongruent</code> object to a new column in <code>aggregating_zones</code> in proportion with the area of overlap: the larger the spatial intersection between input and output features, the larger the corresponding value.
This is implemented in <code><a href="https://r-spatial.github.io/sf/reference/interpolate_aw.html">st_interpolate_aw()</a></code>, as demonstrated in the code chunk below.</p>
<div class="sourceCode" id="cb139"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">iv</span> <span class="op">=</span> <span class="va">incongruent</span><span class="op">[</span><span class="st">"value"</span><span class="op">]</span> <span class="co"># keep only the values to be transferred</span></span>
<span><span class="va">agg_aw</span> <span class="op">=</span> <span class="fu"><a href="https://r-spatial.github.io/sf/reference/interpolate_aw.html">st_interpolate_aw</a></span><span class="op">(</span><span class="va">iv</span>, <span class="va">aggregating_zones</span>, extensive <span class="op">=</span> <span class="cn">TRUE</span><span class="op">)</span></span>
<span><span class="co">#> Warning in st_interpolate_aw.sf(iv, aggregating_zones, extensive = TRUE):</span></span>
<span><span class="co">#> st_interpolate_aw assumes attributes are constant or uniform over areas of x</span></span>
<span><span class="va">agg_aw</span><span class="op">$</span><span class="va">value</span></span>
<span><span class="co">#> [1] 19.6 25.7</span></span></code></pre></div>
<p>In our case it is meaningful to sum up the values of the intersections falling into the aggregating zones since total income is a so-called spatially extensive variable (which increases with area), assuming income is evenly distributed across the smaller zones (hence the warning message above).
This would be different for spatially <a href="https://geodacenter.github.io/workbook/3b_rates/lab3b.html#spatially-extensive-and-spatially-intensive-variables">intensive</a> variables such as <em>average</em> income or percentages, which do not increase as the area increases.
<code><a href="https://r-spatial.github.io/sf/reference/interpolate_aw.html">st_interpolate_aw()</a></code> works equally with spatially intensive variables: set the <code>extensive</code> parameter to <code>FALSE</code> and it will use an average rather than a sum function when doing the aggregation.</p>
</div>
</div>
<div id="spatial-ras" class="section level2" number="4.3">
<h2>
<span class="header-section-number">4.3</span> Spatial operations on raster data<a class="anchor" aria-label="anchor" href="#spatial-ras"><i class="fas fa-link"></i></a>
</h2>
<p>This section builds on Section <a href="attr.html#manipulating-raster-objects">3.3</a>, which highlights various basic methods for manipulating raster datasets, to demonstrate more advanced and explicitly spatial raster operations, and uses the objects <code>elev</code> and <code>grain</code> manually created in Section <a href="attr.html#manipulating-raster-objects">3.3</a>.
For the reader’s convenience, these datasets can be also found in the <strong>spData</strong> package.</p>
<div class="sourceCode" id="cb140"><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">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></code></pre></div>
<div id="spatial-raster-subsetting" class="section level3" number="4.3.1">
<h3>
<span class="header-section-number">4.3.1</span> Spatial subsetting<a class="anchor" aria-label="anchor" href="#spatial-raster-subsetting"><i class="fas fa-link"></i></a>
</h3>
<p>The previous chapter (Section <a href="attr.html#manipulating-raster-objects">3.3</a>) demonstrated how to retrieve values associated with specific cell IDs or row and column combinations.
Raster objects can also be extracted by location (coordinates) and other spatial objects.
To use coordinates for subsetting, one can ‘translate’ the coordinates into a cell ID with the <strong>terra</strong> function <code><a href="https://rspatial.github.io/terra/reference/xyCellFrom.html">cellFromXY()</a></code>.
An alternative is to use <code><a href="https://rspatial.github.io/terra/reference/extract.html">terra::extract()</a></code> (be careful, there is also a function called <code><a href="https://rspatial.github.io/terra/reference/extract.html">extract()</a></code> in the <strong>tidyverse</strong>) to extract values.
Both methods are demonstrated below to find the value of the cell that covers a point located at coordinates of 0.1, 0.1.
</p>
<div class="sourceCode" id="cb141"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">id</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/xyCellFrom.html">cellFromXY</a></span><span class="op">(</span><span class="va">elev</span>, xy <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/matrix.html">matrix</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="fl">0.1</span>, <span class="fl">0.1</span><span class="op">)</span>, ncol <span class="op">=</span> <span class="fl">2</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="va">elev</span><span class="op">[</span><span class="va">id</span><span class="op">]</span></span>
<span><span class="co"># the same as</span></span>
<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">elev</span>, <span class="fu"><a href="https://rdrr.io/r/base/matrix.html">matrix</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="fl">0.1</span>, <span class="fl">0.1</span><span class="op">)</span>, ncol <span class="op">=</span> <span class="fl">2</span><span class="op">)</span><span class="op">)</span></span></code></pre></div>
<p>Raster objects can also be subset with another raster object, as demonstrated in the code chunk below:</p>
<div class="sourceCode" id="cb142"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">clip</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>xmin <span class="op">=</span> <span class="fl">0.9</span>, xmax <span class="op">=</span> <span class="fl">1.8</span>, ymin <span class="op">=</span> <span class="op">-</span><span class="fl">0.45</span>, ymax <span class="op">=</span> <span class="fl">0.45</span>,</span>
<span> resolution <span class="op">=</span> <span class="fl">0.3</span>, vals <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/rep.html">rep</a></span><span class="op">(</span><span class="fl">1</span>, <span class="fl">9</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="va">elev</span><span class="op">[</span><span class="va">clip</span><span class="op">]</span></span>
<span><span class="co"># we can also use extract</span></span>
<span><span class="co"># terra::extract(elev, ext(clip))</span></span></code></pre></div>
<p>This amounts to retrieving the values of the first raster object (in this case <code>elev</code>) that fall within the extent of a second raster (here: <code>clip</code>), as illustrated in Figure <a href="spatial-operations.html#fig:raster-subset">4.10</a>.</p>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:raster-subset"></span>
<img src="figures/04_raster_subset.png" alt="Original raster (left). Raster mask (middle). Output of masking a raster (right)." width="100%"><p class="caption">
FIGURE 4.10: Original raster (left). Raster mask (middle). Output of masking a raster (right).
</p>
</div>
<p>The example above returned the values of specific cells, but in many cases spatial outputs from subsetting operations on raster datasets are needed.
This can be done by setting the <code>drop</code> argument of the <code>[</code> operator to <code>FALSE</code>.
The code below returns the first two cells of <code>elev</code>, i.e., the first two cells of the top row, as a raster object (only the first 2 lines of the output is shown):</p>
<div class="sourceCode" id="cb143"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">elev</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">2</span>, drop <span class="op">=</span> <span class="cn">FALSE</span><span class="op">]</span> <span class="co"># spatial subsetting with cell IDs</span></span>
<span><span class="co">#> class : SpatRaster </span></span>
<span><span class="co">#> dimensions : 1, 2, 1 (nrow, ncol, nlyr)</span></span>
<span><span class="co">#> ...</span></span></code></pre></div>
<p>Another common use case of spatial subsetting is when a raster with <code>logical</code> (or <code>NA</code>) values is used to mask another raster with the same extent and resolution, as illustrated in Figure <a href="spatial-operations.html#fig:raster-subset">4.10</a>.
In this case, the <code>[</code> and <code><a href="https://rspatial.github.io/terra/reference/mask.html">mask()</a></code> functions can be used (results not shown).</p>
<div class="sourceCode" id="cb144"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="co"># create raster mask</span></span>
<span><span class="va">rmask</span> <span class="op">=</span> <span class="va">elev</span></span>
<span><span class="fu"><a href="https://rspatial.github.io/terra/reference/values.html">values</a></span><span class="op">(</span><span class="va">rmask</span><span class="op">)</span> <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/sample.html">sample</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="cn">NA</span>, <span class="cn">TRUE</span><span class="op">)</span>, <span class="fl">36</span>, replace <span class="op">=</span> <span class="cn">TRUE</span><span class="op">)</span></span></code></pre></div>
<p>In the code chunk above, we have created a mask object called <code>rmask</code> with values randomly assigned to <code>NA</code> and <code>TRUE</code>.
Next, we want to keep those values of <code>elev</code> which are <code>TRUE</code> in <code>rmask</code>.
In other words, we want to mask <code>elev</code> with <code>rmask</code>.</p>
<div class="sourceCode" id="cb145"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="co"># spatial subsetting</span></span>
<span><span class="va">elev</span><span class="op">[</span><span class="va">rmask</span>, drop <span class="op">=</span> <span class="cn">FALSE</span><span class="op">]</span> <span class="co"># with [ operator</span></span>
<span><span class="co"># we can also use mask</span></span>
<span><span class="co"># mask(elev, rmask)</span></span></code></pre></div>
<p>The above approach can be also used to replace some values (e.g., expected to be wrong) with NA.</p>
<div class="sourceCode" id="cb146"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">elev</span><span class="op">[</span><span class="va">elev</span> <span class="op"><</span> <span class="fl">20</span><span class="op">]</span> <span class="op">=</span> <span class="cn">NA</span></span></code></pre></div>
<p>These operations are in fact Boolean local operations since we compare cell-wise two rasters.
The next subsection explores these and related operations in more detail.</p>
</div>
<div id="map-algebra" class="section level3" number="4.3.2">
<h3>
<span class="header-section-number">4.3.2</span> Map algebra<a class="anchor" aria-label="anchor" href="#map-algebra"><i class="fas fa-link"></i></a>
</h3>
<p>
The term ‘map algebra’ was coined in the late 1970s to describe a “set of conventions, capabilities, and techniques” for the analysis of geographic raster <em>and</em> (although less prominently) vector data <span class="citation">(<a href="references.html#ref-tomlin_map_1994">Tomlin 1994</a>)</span>.
In this context, we define map algebra more narrowly, as operations that modify or summarize raster cell values, with reference to surrounding cells, zones, or statistical functions that apply to every cell.</p>
<p>Map algebra operations tend to be fast, because raster datasets only implicitly store coordinates, hence the <a href="https://geozoneblog.wordpress.com/2013/04/19/raster-vs-vector/">old adage</a> “raster is faster but vector is corrector”.
The location of cells in raster datasets can be calculated by using its matrix position and the resolution and origin of the dataset (stored in the header).
For the processing, however, the geographic position of a cell is barely relevant as long as we make sure that the cell position is still the same after the processing.
Additionally, if two or more raster datasets share the same extent, projection and resolution, one could treat them as matrices for the processing.</p>
<p>This is the way that map algebra works with the <strong>terra</strong> package.
First, the headers of the raster datasets are queried and (in cases where map algebra operations work on more than one dataset) checked to ensure the datasets are compatible.
Second, map algebra retains the so-called one-to-one locational correspondence, meaning that cells cannot move.
This differs from matrix algebra, in which values change position, for example when multiplying or dividing matrices.</p>
<p>Map algebra (or cartographic modeling with raster data) divides raster operations into four subclasses <span class="citation">(<a href="references.html#ref-tomlin_geographic_1990">Tomlin 1990</a>)</span>, with each working on one or several grids simultaneously:</p>
<ol style="list-style-type: decimal">
<li>
<em>Local</em> or per-cell operations</li>
<li>
<em>Focal</em> or neighborhood operations.
Most often the output cell value is the result of a 3 x 3 input cell block</li>
<li>
<em>Zonal</em> operations are similar to focal operations, but the surrounding pixel grid on which new values are computed can have irregular sizes and shapes</li>
<li>
<em>Global</em> or per-raster operations.
That means the output cell derives its value potentially from one or several entire rasters</li>
</ol>
<p>This typology classifies map algebra operations by the number of cells used for each pixel processing step and the type of the output.
For the sake of completeness, we should mention that raster operations can also be classified by discipline such as terrain, hydrological analysis, or image classification.
The following sections explain how each type of map algebra operations can be used, with reference to worked examples.</p>
</div>
<div id="local-operations" class="section level3" number="4.3.3">
<h3>
<span class="header-section-number">4.3.3</span> Local operations<a class="anchor" aria-label="anchor" href="#local-operations"><i class="fas fa-link"></i></a>
</h3>
<p>
<strong>Local</strong> operations comprise all cell-by-cell operations in one or several layers.
This includes adding or subtracting values from a raster, squaring and multiplying rasters.
Raster algebra also allows logical operations such as finding all raster cells that are greater than a specific value (5 in our example below).
The <strong>terra</strong> package supports all these operations and more, as demonstrated below (Figure <a href="spatial-operations.html#fig:04-local-operations">4.11</a>):</p>
<div class="sourceCode" id="cb147"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">elev</span> <span class="op">+</span> <span class="va">elev</span></span>
<span><span class="va">elev</span><span class="op">^</span><span class="fl">2</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/Log.html">log</a></span><span class="op">(</span><span class="va">elev</span><span class="op">)</span></span>
<span><span class="va">elev</span> <span class="op">></span> <span class="fl">5</span></span></code></pre></div>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:04-local-operations"></span>
<img src="figures/04-local-operations.png" alt="Examples of different local operations of the elev raster object: adding two rasters, squaring, applying logarithmic transformation, and performing a logical operation." width="100%"><p class="caption">
FIGURE 4.11: Examples of different local operations of the elev raster object: adding two rasters, squaring, applying logarithmic transformation, and performing a logical operation.
</p>
</div>
<p>Another good example of local operations is the classification of intervals of numeric values into groups such as grouping a digital elevation model into low (class 1), middle (class 2) and high elevations (class 3).
Using the <code><a href="https://rspatial.github.io/terra/reference/classify.html">classify()</a></code> command, we need first to construct a reclassification matrix, where the first column corresponds to the lower and the second column to the upper end of the class.
The third column represents the new value for the specified ranges in column one and two.</p>
<div class="sourceCode" id="cb148"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">rcl</span> <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/matrix.html">matrix</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="fl">0</span>, <span class="fl">12</span>, <span class="fl">1</span>, <span class="fl">12</span>, <span class="fl">24</span>, <span class="fl">2</span>, <span class="fl">24</span>, <span class="fl">36</span>, <span class="fl">3</span><span class="op">)</span>, ncol <span class="op">=</span> <span class="fl">3</span>, byrow <span class="op">=</span> <span class="cn">TRUE</span><span class="op">)</span></span>
<span><span class="va">rcl</span></span>
<span><span class="co">#> [,1] [,2] [,3]</span></span>
<span><span class="co">#> [1,] 0 12 1</span></span>
<span><span class="co">#> [2,] 12 24 2</span></span>
<span><span class="co">#> [3,] 24 36 3</span></span></code></pre></div>
<p>Here, we assign the raster values in the ranges 0–12, 12–24 and 24–36 are <em>reclassified</em> to take values 1, 2 and 3, respectively.</p>
<div class="sourceCode" id="cb149"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">recl</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/classify.html">classify</a></span><span class="op">(</span><span class="va">elev</span>, rcl <span class="op">=</span> <span class="va">rcl</span><span class="op">)</span></span></code></pre></div>
<p>The <code><a href="https://rspatial.github.io/terra/reference/classify.html">classify()</a></code> function can be also used when we want to reduce the number of classes in our categorical rasters.
We will perform several additional reclassifications in Chapter <a href="location.html#location">14</a>.</p>
<p>Apart from applying arithmetic operators directly, one can also use the <code><a href="https://rspatial.github.io/terra/reference/app.html">app()</a></code>, <code><a href="https://rspatial.github.io/terra/reference/tapp.html">tapp()</a></code> and <code><a href="https://rspatial.github.io/terra/reference/lapp.html">lapp()</a></code> functions.
They are more efficient, hence, they are preferable in the presence of large raster datasets.
Additionally, they allow you to save an output file directly.
The <code><a href="https://rspatial.github.io/terra/reference/app.html">app()</a></code> function applies a function to each cell of a raster and is used to summarize (e.g., calculating the sum) the values of multiple layers into one layer.
<code><a href="https://rspatial.github.io/terra/reference/tapp.html">tapp()</a></code> is an extension of <code><a href="https://rspatial.github.io/terra/reference/app.html">app()</a></code>, allowing us to select a subset of layers (see the <code>index</code> argument) for which we want to perform a certain operation.
Finally, the <code><a href="https://rspatial.github.io/terra/reference/lapp.html">lapp()</a></code> function allows us to apply a function to each cell using layers as arguments – an application of <code><a href="https://rspatial.github.io/terra/reference/lapp.html">lapp()</a></code> is presented below.</p>
<p>The calculation of the normalized difference vegetation index (NDVI) is a well-known local (pixel-by-pixel) raster operation.
It returns a raster with values between -1 and 1; positive values indicate the presence of living plants (mostly > 0.2).
NDVI is calculated from red and near-infrared (NIR) bands of remotely sensed imagery, typically from satellite systems such as Landsat or Sentinel.
Vegetation absorbs light heavily in the visible light spectrum, and especially in the red channel, while reflecting NIR light. Here’s the NVDI formula:</p>
<p><span class="math display">\[
\begin{split}
NDVI&= \frac{\text{NIR} - \text{Red}}{\text{NIR} + \text{Red}}\\
\end{split}
\]</span></p>
<p>Let’s calculate NDVI for the multispectral satellite file of the Zion National Park.</p>
<div class="sourceCode" id="cb150"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">multi_raster_file</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/landsat.tif"</span>, package <span class="op">=</span> <span class="st">"spDataLarge"</span><span class="op">)</span></span>
<span><span class="va">multi_rast</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="va">multi_raster_file</span><span class="op">)</span></span></code></pre></div>
<p>Our raster object has four satellite bands from the Landsat 8 satellite — blue, green, red, and near-infrared (NIR).
Importantly, Landsat level-2 products are stored as integers to save disk space, and thus we need to convert them to floating-point numbers before doing any calculations.
For that purpose, we need to apply a scaling factor (0.0000275) and add an offset (-0.2) to the original values.<a class="footnote-ref" tabindex="0" data-toggle="popover" data-content='<p>You can read more about it at <a href="https://www.usgs.gov/faqs/how-do-i-use-a-scale-factor-landsat-level-2-science-products" class="uri">https://www.usgs.gov/faqs/how-do-i-use-a-scale-factor-landsat-level-2-science-products</a>.</p>'><sup>22</sup></a></p>
<div class="sourceCode" id="cb151"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">multi_rast</span> <span class="op">=</span> <span class="op">(</span><span class="va">multi_rast</span> <span class="op">*</span> <span class="fl">0.0000275</span><span class="op">)</span> <span class="op">-</span> <span class="fl">0.2</span></span></code></pre></div>
<p>The proper values now should be in a range between 0 and 1.
This is not the case here, probably due to the presence of clouds and other atmospheric effects, thus we need to replace below 0 to 0.</p>
<div class="sourceCode" id="cb152"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">multi_rast</span><span class="op">[</span><span class="va">multi_rast</span> <span class="op"><</span> <span class="fl">0</span><span class="op">]</span> <span class="op">=</span> <span class="fl">0</span></span></code></pre></div>
<p>The next step should be to implement the NDVI formula into an R function:</p>
<div class="sourceCode" id="cb153"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">ndvi_fun</span> <span class="op">=</span> <span class="kw">function</span><span class="op">(</span><span class="va">nir</span>, <span class="va">red</span><span class="op">)</span><span class="op">{</span></span>
<span> <span class="op">(</span><span class="va">nir</span> <span class="op">-</span> <span class="va">red</span><span class="op">)</span> <span class="op">/</span> <span class="op">(</span><span class="va">nir</span> <span class="op">+</span> <span class="va">red</span><span class="op">)</span></span>
<span><span class="op">}</span></span></code></pre></div>
<p>This function accepts two numerical arguments, <code>nir</code> and <code>red</code>, and returns a numerical vector with NDVI values.
It can be used as the <code>fun</code> argument of <code><a href="https://rspatial.github.io/terra/reference/lapp.html">lapp()</a></code>.
We just need to remember that our function expects two bands (not four from the original raster), and they need to be in the NIR, red order.
That is why we subset the input raster with <code>multi_rast[[c(4, 3)]]</code> before doing any calculations.</p>
<div class="sourceCode" id="cb154"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">ndvi_rast</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/lapp.html">lapp</a></span><span class="op">(</span><span class="va">multi_rast</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">4</span>, <span class="fl">3</span><span class="op">)</span><span class="op">]</span><span class="op">]</span>, fun <span class="op">=</span> <span class="va">ndvi_fun</span><span class="op">)</span></span></code></pre></div>
<p>The result, shown on the right panel in Figure <a href="spatial-operations.html#fig:04-ndvi">4.12</a>, can be compared to the RGB image of the same area (left panel of the same Figure).
It allows us to see that the largest NDVI values are connected to northern areas of dense forest, while the lowest values are related to the lake in the north and snowy mountain ridges.</p>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:04-ndvi"></span>
<img src="figures/04-ndvi.png" alt="RGB image (left) and NDVI values (right) calculated for the example satellite file of the Zion National Park" width="100%"><p class="caption">
FIGURE 4.12: RGB image (left) and NDVI values (right) calculated for the example satellite file of the Zion National Park
</p>
</div>
<p>Predictive mapping is another interesting application of local raster operations.
The response variable corresponds to measured or observed points in space, for example, species richness, the presence of landslides, tree disease or crop yield.
Consequently, we can easily retrieve space- or airborne predictor variables from various rasters (elevation, pH, precipitation, temperature, land cover, soil class, etc.).
Subsequently, we model our response as a function of our predictors using <code><a href="https://rdrr.io/r/stats/lm.html">lm()</a></code>, <code><a href="https://rdrr.io/r/stats/glm.html">glm()</a></code>, <code>gam()</code> or a machine-learning technique.
Spatial predictions on raster objects can therefore be made by applying estimated coefficients to the predictor raster values, and summing the output raster values (see Chapter <a href="eco.html#eco">15</a>).</p>
</div>
<div id="focal-operations" class="section level3" number="4.3.4">
<h3>
<span class="header-section-number">4.3.4</span> Focal operations<a class="anchor" aria-label="anchor" href="#focal-operations"><i class="fas fa-link"></i></a>
</h3>
<p>
While local functions operate on one cell, though possibly from multiple layers, <strong>focal</strong> operations take into account a central (focal) cell and its neighbors.
The neighborhood (also named kernel, filter or moving window) under consideration is typically of size 3-by-3 cells (that is the central cell and its eight surrounding neighbors), but can take on any other size or (not necessarily rectangular) shape as defined by the user.
A focal operation applies an aggregation function to all cells within the specified neighborhood, uses the corresponding output as the new value for the central cell, and moves on to the next central cell (Figure <a href="spatial-operations.html#fig:focal-example">4.13</a>).
Other names for this operation are spatial filtering and convolution <span class="citation">(<a href="references.html#ref-burrough_principles_2015">Burrough, McDonnell, and Lloyd 2015</a>)</span>.</p>
<p>In R, we can use the <code><a href="https://rspatial.github.io/terra/reference/focal.html">focal()</a></code> function to perform spatial filtering.
We define the shape of the moving window with a <code>matrix</code> whose values correspond to weights (see <code>w</code> parameter in the code chunk below).
Secondly, the <code>fun</code> parameter lets us specify the function we wish to apply to this neighborhood.
Here, we choose the minimum, but any other summary function, including <code><a href="https://rdrr.io/r/base/sum.html">sum()</a></code>, <code><a href="https://rspatial.github.io/terra/reference/summarize-generics.html">mean()</a></code>, or <code><a href="https://rdrr.io/r/stats/cor.html">var()</a></code> can be used.</p>
<div class="sourceCode" id="cb155"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">r_focal</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/focal.html">focal</a></span><span class="op">(</span><span class="va">elev</span>, w <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/matrix.html">matrix</a></span><span class="op">(</span><span class="fl">1</span>, nrow <span class="op">=</span> <span class="fl">3</span>, ncol <span class="op">=</span> <span class="fl">3</span><span class="op">)</span>, fun <span class="op">=</span> <span class="va">min</span><span class="op">)</span></span></code></pre></div>
<p>This function also accepts additional arguments, for example, should it remove NAs in the process (<code>na.rm = TRUE</code>) or not (<code>na.rm = FALSE</code>).</p>
<div class="figure" style="text-align: center">
<span style="display:block;" id="fig:focal-example"></span>
<img src="figures/04_focal_example.png" alt="Input raster (left) and resulting output raster (right) due to a focal operation - finding the minimum value in 3-by-3 moving windows." width="100%"><p class="caption">
FIGURE 4.13: Input raster (left) and resulting output raster (right) due to a focal operation - finding the minimum value in 3-by-3 moving windows.
</p>
</div>
<p>We can quickly check if the output meets our expectations.
In our example, the minimum value has to be always the upper left corner of the moving window (remember we have created the input raster by row-wise incrementing the cell values by one starting at the upper left corner).
In this example, the weighting matrix consists only of 1s, meaning each cell has the same weight on the output, but this can be changed.</p>
<p>Focal functions or filters play a dominant role in image processing.
Low-pass or smoothing filters use the mean function to remove extremes.
In the case of categorical data, we can replace the mean with the mode, which is the most common value.
By contrast, high-pass filters accentuate features.
The line detection Laplace and Sobel filters might serve as an example here.
Check the <code><a href="https://rspatial.github.io/terra/reference/focal.html">focal()</a></code> help page for how to use them in R (this will also be used in the exercises at the end of this chapter).</p>
<p>Terrain processing, the calculation of topographic characteristics such as slope, aspect and flow directions, relies on focal functions.
<code><a href="https://rspatial.github.io/terra/reference/terrain.html">terrain()</a></code> can be used to calculate these metrics, although some terrain algorithms, including the Zevenbergen and Thorne method to compute slope, are not implemented in this <strong>terra</strong> function.
Many other algorithms — including curvatures, contributing areas and wetness indices — are implemented in open source desktop geographic information system (GIS) software.
Chapter <a href="gis.html#gis">10</a> shows how to access such GIS functionality from within R.</p>
</div>
<div id="zonal-operations" class="section level3" number="4.3.5">
<h3>
<span class="header-section-number">4.3.5</span> Zonal operations<a class="anchor" aria-label="anchor" href="#zonal-operations"><i class="fas fa-link"></i></a>
</h3>
<p>
Just like focal operations, <em>zonal</em> operations apply an aggregation function to multiple raster cells.
However, a second raster, usually with categorical values, defines the <em>zonal filters</em> (or ‘zones’) in the case of zonal operations, as opposed to a neighborhood window in the case of focal operations presented in the previous section.
Consequently, raster cells defining the zonal filter do not necessarily have to be neighbors.
Our grain size raster is a good example, as illustrated in the right panel of Figure <a href="attr.html#fig:cont-raster">3.2</a>: different grain sizes are spread irregularly throughout the raster.
Finally, the result of a zonal operation is a summary table grouped by zone which is why this operation is also known as <em>zonal statistics</em> in the GIS world.
This is in contrast to focal operations which return a raster object by default.</p>
<p>The following code chunk uses the <code><a href="https://rspatial.github.io/terra/reference/zonal.html">zonal()</a></code> function to calculate the mean elevation associated with each grain size class.</p>
<div class="sourceCode" id="cb156"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">z</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/zonal.html">zonal</a></span><span class="op">(</span><span class="va">elev</span>, <span class="va">grain</span>, fun <span class="op">=</span> <span class="st">"mean"</span><span class="op">)</span></span>
<span><span class="va">z</span></span>
<span><span class="co">#> grain elev</span></span>
<span><span class="co">#> 1 clay 14.8</span></span>
<span><span class="co">#> 2 silt 21.2</span></span>
<span><span class="co">#> 3 sand 18.7</span></span></code></pre></div>
<p>This returns the statistics for each category, here the mean altitude for each grain size class.
Note that it is also possible to get a raster with calculated statistics for each zone by setting the <code>as.raster</code> argument to <code>TRUE</code>.</p>
</div>
<div id="global-operations-and-distances" class="section level3" number="4.3.6">
<h3>
<span class="header-section-number">4.3.6</span> Global operations and distances<a class="anchor" aria-label="anchor" href="#global-operations-and-distances"><i class="fas fa-link"></i></a>
</h3>
<p><em>Global</em> operations are a special case of zonal operations with the entire raster dataset representing a single zone.
The most common global operations are descriptive statistics for the entire raster dataset such as the minimum or maximum – we already discussed those in Section <a href="attr.html#summarizing-raster-objects">3.3.2</a>.</p>
<p>Aside from that, global operations are also useful for the computation of distance and weight rasters.
In the first case, one can calculate the distance from each cell to a specific target cell.
For example, one might want to compute the distance to the nearest coast (see also <code><a href="https://rspatial.github.io/terra/reference/distance.html">terra::distance()</a></code>).
We might also want to consider topography, that means, we are not only interested in the pure distance but would like also to avoid the crossing of mountain ranges when going to the coast.
To do so, we can weight the distance with elevation so that each additional altitudinal meter ‘prolongs’ the Euclidean distance (in Exercises 8 and 9 at the end of this chapter you will do exactly that).
Visibility and viewshed computations also belong to the family of global operations (in the exercises of Chapter <a href="gis.html#gis">10</a>, you will compute a viewshed raster).</p>
</div>
<div id="map-algebra-counterparts-in-vector-processing" class="section level3" number="4.3.7">
<h3>
<span class="header-section-number">4.3.7</span> Map algebra counterparts in vector processing<a class="anchor" aria-label="anchor" href="#map-algebra-counterparts-in-vector-processing"><i class="fas fa-link"></i></a>
</h3>
<p>Many map algebra operations have a counterpart in vector processing <span class="citation">(<a href="references.html#ref-liu_essential_2009">Liu and Mason 2009</a>)</span>.
Computing a distance raster (global operation) while only considering a maximum distance (logical focal operation) is the equivalent to a vector buffer operation (Section <a href="geometry-operations.html#clipping">5.2.5</a>).
Reclassifying raster data (either local or zonal function depending on the input) is equivalent to dissolving vector data (Section <a href="spatial-operations.html#spatial-joining">4.2.5</a>).
Overlaying two rasters (local operation), where one contains <code>NULL</code> or <code>NA</code> values representing a mask, is similar to vector clipping (Section <a href="geometry-operations.html#clipping">5.2.5</a>).
Quite similar to spatial clipping is intersecting two layers (Section <a href="spatial-operations.html#spatial-subsetting">4.2.1</a>).
The difference is that these two layers (vector or raster) simply share an overlapping area (see Figure <a href="geometry-operations.html#fig:venn-clip">5.8</a> for an example).
However, be careful with the wording.
Sometimes the same words have slightly different meanings for raster and vector data models.
While aggregating polygon geometries means dissolving boundaries, for raster data geometries it means increasing cell sizes and thereby reducing spatial resolution.
Zonal operations dissolve the cells of one raster in accordance with the zones (categories) of another raster dataset using an aggregating function.</p>
</div>
<div id="merging-rasters" class="section level3" number="4.3.8">
<h3>
<span class="header-section-number">4.3.8</span> Merging rasters<a class="anchor" aria-label="anchor" href="#merging-rasters"><i class="fas fa-link"></i></a>
</h3>
<p>
Suppose we would like to compute the NDVI (see Section <a href="spatial-operations.html#local-operations">4.3.3</a>), and additionally want to compute terrain attributes from elevation data for observations within a study area.
Such computations rely on remotely sensed information.
The corresponding imagery is often divided into scenes covering a specific spatial extent, and frequently, a study area covers more than one scene.
Then, we would need to merge the scenes covered by our study area.
In the easiest case, we can just merge these scenes, that is put them side by side.
This is possible, for example, with digital elevation data.
In the following code chunk we first download the SRTM elevation data for Austria and Switzerland (for the country codes, see the <strong>geodata</strong> function <code>country_codes()</code>).
In a second step, we merge the two rasters into one.</p>
<div class="sourceCode" id="cb157"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">aut</span> <span class="op">=</span> <span class="fu">geodata</span><span class="fu">::</span><span class="fu"><a href="https://rdrr.io/pkg/geodata/man/elevation.html">elevation_30s</a></span><span class="op">(</span>country <span class="op">=</span> <span class="st">"AUT"</span>, path <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/tempfile.html">tempdir</a></span><span class="op">(</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="va">ch</span> <span class="op">=</span> <span class="fu">geodata</span><span class="fu">::</span><span class="fu"><a href="https://rdrr.io/pkg/geodata/man/elevation.html">elevation_30s</a></span><span class="op">(</span>country <span class="op">=</span> <span class="st">"CHE"</span>, path <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/tempfile.html">tempdir</a></span><span class="op">(</span><span class="op">)</span><span class="op">)</span></span>
<span><span class="va">aut_ch</span> <span class="op">=</span> <span class="fu"><a href="https://rspatial.github.io/terra/reference/merge.html">merge</a></span><span class="op">(</span><span class="va">aut</span>, <span class="va">ch</span><span class="op">)</span></span></code></pre></div>
<p><strong>terra</strong>’s <code><a href="https://rspatial.github.io/terra/reference/merge.html">merge()</a></code> command combines two images, and in case they overlap, it uses the value of the first raster.</p>
<p>The merging approach is of little use when the overlapping values do not correspond to each other.
This is frequently the case when you want to combine spectral imagery from scenes that were taken on different dates.
The <code><a href="https://rspatial.github.io/terra/reference/merge.html">merge()</a></code> command will still work but you will see a clear border in the resulting image.
On the other hand, the <code><a href="https://rspatial.github.io/terra/reference/mosaic.html">mosaic()</a></code> command lets you define a function for the overlapping area.
For instance, we could compute the mean value – this might smooth the clear border in the merged result but it will most likely not make it disappear.
For a more detailed introduction to remote sensing with R, see <span class="citation">Wegmann, Leutner, and Dech (<a href="references.html#ref-wegmann_remote_2016">2016</a>)</span>.</p>
</div>
</div>
<div id="exercises-2" class="section level2" number="4.4">
<h2>
<span class="header-section-number">4.4</span> Exercises<a class="anchor" aria-label="anchor" href="#exercises-2"><i class="fas fa-link"></i></a>
</h2>
<div class="sourceCode" id="cb158"><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://dplyr.tidyverse.org">dplyr</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></code></pre></div>
<p>E1. It was established in Section <a href="spatial-operations.html#spatial-vec">4.2</a> that Canterbury was the region of New Zealand containing most of the 101 highest points in the country.
How many of these high points does the Canterbury region contain?</p>
<p><strong>Bonus:</strong> plot the result using the <code><a href="https://rspatial.github.io/terra/reference/plot.html">plot()</a></code> function to show all of New Zealand, <code>canterbury</code> region highlighted in yellow, high points in Canterbury represented by red crosses (hint: <code>pch = 7</code>) and high points in other parts of New Zealand represented by blue circles. See the help page <code><a href="https://rspatial.github.io/terra/reference/lines.html">?points</a></code> for details with an illustration of different <code>pch</code> values.</p>
<p>E2. Which region has the second highest number of <code>nz_height</code> points, and how many does it have?</p>
<p>E3. Generalizing the question to all regions: how many of New Zealand’s 16 regions contain points which belong to the top 101 highest points in the country? Which regions?</p>
<ul>
<li>Bonus: create a table listing these regions in order of the number of points and their name.</li>
</ul>
<p>E4. Test your knowledge of spatial predicates by finding out and plotting how US states relate to each other and other spatial objects.</p>
<p>The starting point of this exercise is to create an object representing Colorado state in the USA. Do this with the command
<code>colorado = us_states[us_states$NAME == "Colorado",]</code> (base R) or with with the <code><a href="https://dplyr.tidyverse.org/reference/filter.html">filter()</a></code> function (tidyverse) and plot the resulting object in the context of US states.</p>
<ul>
<li>Create a new object representing all the states that geographically intersect with Colorado and plot the result (hint: the most concise way to do this is with the subsetting method <code>[</code>).</li>
<li>Create another object representing all the objects that touch (have a shared boundary with) Colorado and plot the result (hint: remember you can use the argument <code>op = st_intersects</code> and other spatial relations during spatial subsetting operations in base R).</li>
<li>Bonus: create a straight line from the centroid of the District of Columbia near the East coast to the centroid of California near the West coast of the USA (hint: functions <code><a href="https://r-spatial.github.io/sf/reference/geos_unary.html">st_centroid()</a></code>, <code><a href="https://r-spatial.github.io/sf/reference/geos_combine.html">st_union()</a></code> and <code><a href="https://r-spatial.github.io/sf/reference/st_cast.html">st_cast()</a></code> described in Chapter 5 may help) and identify which states this long East-West line crosses.</li>
</ul>
<p>E5. Use <code>dem = rast(system.file("raster/dem.tif", package = "spDataLarge"))</code>, and reclassify the elevation in three classes: low (<300), medium and high (>500).
Secondly, read the NDVI raster (<code>ndvi = rast(system.file("raster/ndvi.tif", package = "spDataLarge"))</code>) and compute the mean NDVI and the mean elevation for each altitudinal class.</p>
<p>E6. Apply a line detection filter to <code>rast(system.file("ex/logo.tif", package = "terra"))</code>.
Plot the result.
Hint: Read <code>?terra::focal()</code>.</p>
<p>E7. Calculate the Normalized Difference Water Index (NDWI; <code>(green - nir)/(green + nir)</code>) of a Landsat image.
Use the Landsat image provided by the <strong>spDataLarge</strong> package (<code>system.file("raster/landsat.tif", package = "spDataLarge")</code>).
Also, calculate a correlation between NDVI and NDWI for this area (hint: you can use the <code><a href="https://rspatial.github.io/terra/reference/layerCor.html">layerCor()</a></code> function).</p>
<p>E8. A StackOverflow <a href="https://stackoverflow.com/questions/35555709/global-raster-of-geographic-distances">post</a> shows how to compute distances to the nearest coastline using <code><a href="https://rspatial.github.io/terra/reference/distance.html">raster::distance()</a></code>.
Try to do something similar but with <code><a href="https://rspatial.github.io/terra/reference/distance.html">terra::distance()</a></code>: retrieve a digital elevation model of Spain, and compute a raster which represents distances to the coast across the country (hint: use <code><a href="https://rdrr.io/pkg/geodata/man/elevation.html">geodata::elevation_30s()</a></code>).
Convert the resulting distances from meters to kilometers.
Note: it may be wise to increase the cell size of the input raster to reduce compute time during this operation (<code><a href="https://rspatial.github.io/terra/reference/aggregate.html">aggregate()</a></code>).</p>
<p>E9. Try to modify the approach used in the above exercise by weighting the distance raster with the elevation raster; every 100 altitudinal meters should increase the distance to the coast by 10 km.
Next, compute and visualize the difference between the raster created using the Euclidean distance (E7) and the raster weighted by elevation.</p>
</div>
</div>
<div class="chapter-nav">
<div class="prev"><a href="attr.html"><span class="header-section-number">3</span> Attribute data operations</a></div>
<div class="next"><a href="geometry-operations.html"><span class="header-section-number">5</span> Geometry operations</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="#spatial-operations"><span class="header-section-number">4</span> Spatial data operations</a></li>
<li><a class="nav-link" href="#prerequisites-2">Prerequisites</a></li>
<li><a class="nav-link" href="#introduction-1"><span class="header-section-number">4.1</span> Introduction</a></li>
<li>
<a class="nav-link" href="#spatial-vec"><span class="header-section-number">4.2</span> Spatial operations on vector data</a><ul class="nav navbar-nav">
<li><a class="nav-link" href="#spatial-subsetting"><span class="header-section-number">4.2.1</span> Spatial subsetting</a></li>
<li><a class="nav-link" href="#topological-relations"><span class="header-section-number">4.2.2</span> Topological relations</a></li>
<li><a class="nav-link" href="#distance-relations"><span class="header-section-number">4.2.3</span> Distance relations</a></li>
<li><a class="nav-link" href="#DE-9IM-strings"><span class="header-section-number">4.2.4</span> DE-9IM strings</a></li>
<li><a class="nav-link" href="#spatial-joining"><span class="header-section-number">4.2.5</span> Spatial joining</a></li>
<li><a class="nav-link" href="#non-overlapping-joins"><span class="header-section-number">4.2.6</span> Distance-based joins</a></li>
<li><a class="nav-link" href="#spatial-aggr"><span class="header-section-number">4.2.7</span> Spatial aggregation</a></li>
<li><a class="nav-link" href="#incongruent"><span class="header-section-number">4.2.8</span> Joining incongruent layers</a></li>
</ul>
</li>
<li>
<a class="nav-link" href="#spatial-ras"><span class="header-section-number">4.3</span> Spatial operations on raster data</a><ul class="nav navbar-nav">
<li><a class="nav-link" href="#spatial-raster-subsetting"><span class="header-section-number">4.3.1</span> Spatial subsetting</a></li>
<li><a class="nav-link" href="#map-algebra"><span class="header-section-number">4.3.2</span> Map algebra</a></li>
<li><a class="nav-link" href="#local-operations"><span class="header-section-number">4.3.3</span> Local operations</a></li>
<li><a class="nav-link" href="#focal-operations"><span class="header-section-number">4.3.4</span> Focal operations</a></li>
<li><a class="nav-link" href="#zonal-operations"><span class="header-section-number">4.3.5</span> Zonal operations</a></li>
<li><a class="nav-link" href="#global-operations-and-distances"><span class="header-section-number">4.3.6</span> Global operations and distances</a></li>
<li><a class="nav-link" href="#map-algebra-counterparts-in-vector-processing"><span class="header-section-number">4.3.7</span> Map algebra counterparts in vector processing</a></li>
<li><a class="nav-link" href="#merging-rasters"><span class="header-section-number">4.3.8</span> Merging rasters</a></li>
</ul>
</li>
<li><a class="nav-link" href="#exercises-2"><span class="header-section-number">4.4</span> Exercises</a></li>
</ul>