-
Notifications
You must be signed in to change notification settings - Fork 1
/
rv.qmd
117 lines (102 loc) · 2.54 KB
/
rv.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
---
title: "Raster or vector?"
author: "Edzer Pebesma"
date: "August 30, 2022"
---
What are the differences between the following two datasets, `r` and `v`?
```{r}
library(sf)
library(stars)
bb = st_bbox(c(xmin = 0, xmax = 10, ymin = 0, ymax = 10))
r = st_as_stars(bb, nx = 10, ny = 10)
r$values = rnorm(100)
plot(r, breaks = "equal", axes = TRUE, col = sf.colors(), reset = FALSE)
box()
```
```{r}
v = st_as_sf(r)
plot(v, axes = TRUE)
```
Let's see how the objects are printed:
```{r}
r
v
```
Let's see how they behave when we intersect a point inside a cell:
```{r}
pt = st_as_sfc("POINT(0.5 0.5)")
x = st_intersects(r, pt, as_points = FALSE)
t(x)
st_intersects(pt, v)
plot(v, axes = TRUE, reset = FALSE)
plot(v[91,], col = 'grey', add = TRUE)
plot(pt, add = TRUE, col = 'red', pch = 16, cex = 1.5)
text(.5, .8, "91")
```
... with a point on a cell border:
```{r}
pt = st_as_sfc("POINT(1 0.5)")
x = st_intersects(r, pt, as_points = FALSE)
t(x)
y = st_intersects(pt, v)
plot(v, axes = TRUE, reset = FALSE)
plot(v[91:92,], col = 'grey', add = TRUE)
plot(pt, add = TRUE, col = 'red', pch = 16, cex = 1.5)
text(c(.5, 1.5), c(.5, .5), c("91", "92"))
```
... with a point on a cell corner:
```{r}
pt = st_as_sfc("POINT(1 1)")
x = st_intersects(r, pt, as_points = FALSE)
t(x)
y = st_intersects(pt, v)
y
plot(v, axes = TRUE, reset = FALSE)
plot(v[c(81,82,91,92),], col = 'grey', add = TRUE)
plot(pt, add = TRUE, col = 'red', pch = 16, cex = 1.5)
text(c(.5, 1.5, .5, 1.5), c(1.5, 1.5, .5, .5), c("81", "82", "91", "92"))
```
... with s2:
```{r}
st_crs(r) = 4326
st_crs(v) = 4326
pt = st_as_sfc("POINT(1 0.5)")
st_crs(pt) = 4326
y = st_intersects(pt, v)
y
plot(v, axes = TRUE, reset = FALSE)
plot(v[91:92,], col = 'grey', add = TRUE)
plot(pt, add = TRUE, col = 'red', pch = 16, cex = 1.5)
text(c(.5, 1.5), c(.5, .5), c("91", "92"))
pt = st_as_sfc("POINT(0.5 1)")
st_crs(pt) = 4326
y = st_intersects(pt, v)
y
```
Southern Hemisphere:
```{r}
bb = st_bbox(c(xmin = 0, xmax = 10, ymin = -10, ymax = 0))
r_s = st_as_stars(bb, nx = 10, ny = 10)
r_s$values = r$values
st_crs(r_s) = 4326
plot(r_s, breaks = "equal", axes = TRUE, col = sf.colors(), reset = FALSE)
box()
pt = st_as_sfc("POINT(0.5 -1)")
st_crs(pt) = 4326
v = st_as_sf(r_s)
st_intersects(pt, v)
```
```{r}
plot(v, axes = TRUE, reset = FALSE)
plot(v[1,], col = 'grey', add = TRUE)
plot(pt, add = TRUE, col = 'red', pch = 16, cex = 1.5)
text(c(.5, .5), c(-.5, -1.5), c("1", "11"))
```
Removing CRS assumes $R^2$:
```{r}
st_crs(pt) = NA
st_crs(v) = NA
st_intersects(pt, v)
st_crs(r_s) = NA
st_intersects(r_s, pt, as_points = FALSE) |> t()
```