Geographic Information Systems Asked on January 19, 2021
I create a SpatialLines from x and y coordinates. I was wondering is there a way to pick 5, 10 or m number points on the network which are equal distance to each other at least (2:m-1) ones have equal distance.
I was thinking maybe I can compute the distance for each point to the previous one and get the cumulative length and use seq(., ., length.out=m)
would give me equal distance but I cannot get x and y coordinates from there.
> Y1
xco yco
1 172868.6 3376152
2 172891.0 3376130
3 172926.8 3376096
4 172949.2 3376074
5 172985.0 3376040
6 173007.3 3376019
7 173029.7 3375997
8 173065.5 3375963
9 173087.9 3375941
10 173123.7 3375907
11 173146.0 3375886
12 173181.9 3375851
13 173204.2 3375830
14 173226.6 3375808
15 173262.4 3375774
16 173320.7 3375718
17 173343.0 3375697
18 173365.4 3375676
19 173401.3 3375641
20 173423.7 3375620
21 173459.6 3375586
22 173482.0 3375564
23 173504.3 3375543
24 173541.2 3375509
25 173563.5 3375489
26 173601.2 3375454
27 173623.4 3375434
28 173644.9 3375415
29 173684.7 3375381
30 173706.7 3375362
31 173747.3 3375328
32 173769.3 3375310
33 173811.0 3375276
34 173832.8 3375259
35 173875.2 3375226
36 173896.9 3375209
37 173918.5 3375192
38 173962.0 3375159
39 173983.6 3375142
40 174027.1 3375109
41 174048.6 3375093
42 174092.1 3375060
43 174113.7 3375043
44 174157.2 3375010
45 174178.8 3374994
46 174222.2 3374960
47 174243.7 3374944
48 174287.1 3374911
49 174308.7 3374894
50 174330.3 3374878
51 174373.6 3374844
52 174395.2 3374828
53 174438.9 3374795
54 174460.4 3374778
55 174504.2 3374745
56 174525.7 3374729
57 174569.4 3374696
58 174590.9 3374680
59 174634.5 3374646
60 174656.0 3374630
61 174700.4 3374597
62 174722.0 3374581
63 174742.4 3374567
64 174790.4 3374536
65 174860.9 3374494
66 174914.5 3374464
67 174932.6 3374455
68 174988.4 3374428
69 175005.4 3374420
70 175063.6 3374395
71 175079.2 3374389
72 175138.5 3374366
73 175200.4 3374346
74 175213.2 3374341
75 175276.0 3374324
76 175340.7 3374307
77 175405.8 3374292
78 175415.9 3374289
79 175481.0 3374274
80 175546.1 3374258
81 175611.2 3374243
82 175621.3 3374240
83 175686.4 3374225
84 175751.5 3374209
85 175761.6 3374207
86 175826.3 3374192
87 175893.1 3374181
88 175960.8 3374172
89 176029.4 3374166
90 176098.8 3374162
91 176168.1 3374161
92 176237.6 3374163
93 176307.1 3374167
94 176375.5 3374175
95 176442.9 3374184
96 176509.2 3374197
97 176517.7 3374198
98 176584.2 3374213
99 176648.6 3374230
100 176712.4 3374248
myLine190=Line(Y1)
myL190<- Lines(list(myLine190), ID = 1)
Spt1 <- SpatialLines(list(myL190))
# Cumulative distance
DD1=sqrt(rowSums(do.call("rbind",lapply(1:(dim(Y1)[1]-1),function(i){coordinates(Y1[i,])-coordinates(Y1[i+1,])}))^2))
cd1=cumsum(c(0,DD1))
grid=seq(min(cd1),max(cd1),length.out = 20)
Y1=read.table(text = ' xco yco
1 172868.6 3376152
2 172891.0 3376130
3 172926.8 3376096
4 172949.2 3376074
5 172985.0 3376040
6 173007.3 3376019
7 173029.7 3375997
8 173065.5 3375963
9 173087.9 3375941
10 173123.7 3375907
11 173146.0 3375886
12 173181.9 3375851
13 173204.2 3375830
14 173226.6 3375808
15 173262.4 3375774
16 173320.7 3375718
17 173343.0 3375697
18 173365.4 3375676
19 173401.3 3375641
20 173423.7 3375620
21 173459.6 3375586
22 173482.0 3375564
23 173504.3 3375543
24 173541.2 3375509
25 173563.5 3375489
26 173601.2 3375454
27 173623.4 3375434
28 173644.9 3375415
29 173684.7 3375381
30 173706.7 3375362
31 173747.3 3375328
32 173769.3 3375310
33 173811.0 3375276
34 173832.8 3375259
35 173875.2 3375226
36 173896.9 3375209
37 173918.5 3375192
38 173962.0 3375159
39 173983.6 3375142
40 174027.1 3375109
41 174048.6 3375093
42 174092.1 3375060
43 174113.7 3375043
44 174157.2 3375010
45 174178.8 3374994
46 174222.2 3374960
47 174243.7 3374944
48 174287.1 3374911
49 174308.7 3374894
50 174330.3 3374878
51 174373.6 3374844
52 174395.2 3374828
53 174438.9 3374795
54 174460.4 3374778
55 174504.2 3374745
56 174525.7 3374729
57 174569.4 3374696
58 174590.9 3374680
59 174634.5 3374646
60 174656.0 3374630
61 174700.4 3374597
62 174722.0 3374581
63 174742.4 3374567
64 174790.4 3374536
65 174860.9 3374494
66 174914.5 3374464
67 174932.6 3374455
68 174988.4 3374428
69 175005.4 3374420
70 175063.6 3374395
71 175079.2 3374389
72 175138.5 3374366
73 175200.4 3374346
74 175213.2 3374341
75 175276.0 3374324
76 175340.7 3374307
77 175405.8 3374292
78 175415.9 3374289
79 175481.0 3374274
80 175546.1 3374258
81 175611.2 3374243
82 175621.3 3374240
83 175686.4 3374225
84 175751.5 3374209
85 175761.6 3374207
86 175826.3 3374192
87 175893.1 3374181
88 175960.8 3374172
89 176029.4 3374166
90 176098.8 3374162
91 176168.1 3374161
92 176237.6 3374163
93 176307.1 3374167
94 176375.5 3374175
95 176442.9 3374184
96 176509.2 3374197
97 176517.7 3374198
98 176584.2 3374213
99 176648.6 3374230
100 176712.4 3374248', header = TRUE)
I would solve your problem as follows:
# packages
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(lwgeom)
#> Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.8.0, PROJ 6.3.1
# data
my_data <- read.table(
text = " xco yco
1 172868.6 3376152
2 172891.0 3376130
3 172926.8 3376096
4 172949.2 3376074
5 172985.0 3376040
6 173007.3 3376019
7 173029.7 3375997
8 173065.5 3375963
9 173087.9 3375941
10 173123.7 3375907
11 173146.0 3375886
12 173181.9 3375851
13 173204.2 3375830
14 173226.6 3375808
15 173262.4 3375774
16 173320.7 3375718
17 173343.0 3375697
18 173365.4 3375676
19 173401.3 3375641
20 173423.7 3375620
21 173459.6 3375586
22 173482.0 3375564
23 173504.3 3375543
24 173541.2 3375509
25 173563.5 3375489
26 173601.2 3375454
27 173623.4 3375434
28 173644.9 3375415
29 173684.7 3375381
30 173706.7 3375362
31 173747.3 3375328
32 173769.3 3375310
33 173811.0 3375276
34 173832.8 3375259
35 173875.2 3375226
36 173896.9 3375209
37 173918.5 3375192
38 173962.0 3375159
39 173983.6 3375142
40 174027.1 3375109
41 174048.6 3375093
42 174092.1 3375060
43 174113.7 3375043
44 174157.2 3375010
45 174178.8 3374994
46 174222.2 3374960
47 174243.7 3374944
48 174287.1 3374911
49 174308.7 3374894
50 174330.3 3374878
51 174373.6 3374844
52 174395.2 3374828
53 174438.9 3374795
54 174460.4 3374778
55 174504.2 3374745
56 174525.7 3374729
57 174569.4 3374696
58 174590.9 3374680
59 174634.5 3374646
60 174656.0 3374630
61 174700.4 3374597
62 174722.0 3374581
63 174742.4 3374567
64 174790.4 3374536
65 174860.9 3374494
66 174914.5 3374464
67 174932.6 3374455
68 174988.4 3374428
69 175005.4 3374420
70 175063.6 3374395
71 175079.2 3374389
72 175138.5 3374366
73 175200.4 3374346
74 175213.2 3374341
75 175276.0 3374324
76 175340.7 3374307
77 175405.8 3374292
78 175415.9 3374289
79 175481.0 3374274
80 175546.1 3374258
81 175611.2 3374243
82 175621.3 3374240
83 175686.4 3374225
84 175751.5 3374209
85 175761.6 3374207
86 175826.3 3374192
87 175893.1 3374181
88 175960.8 3374172
89 176029.4 3374166
90 176098.8 3374162
91 176168.1 3374161
92 176237.6 3374163
93 176307.1 3374167
94 176375.5 3374175
95 176442.9 3374184
96 176509.2 3374197
97 176517.7 3374198
98 176584.2 3374213
99 176648.6 3374230
100 176712.4 3374248",
header = TRUE
)
Convert data into LINESTRING
format
(my_linestring <- st_linestring(x = as.matrix(my_data)))
#> LINESTRING (172868.6 3376152, 172891 3376130, 172926.8 3376096, 172949.2 3376074, 172985 3376040, 173007.3 3376019, 173029.7 3375997, 173065.5 3375963, 173087.9 3375941, 173123.7 3375907, 173146 3375886, 173181.9 3375851, 173204.2 3375830, 173226.6 3375808, 173262.4 3375774, 173320.7 3375718, 173343 3375697, 173365.4 3375676, 173401.3 3375641, 173423.7 3375620, 173459.6 3375586, 173482 3375564, 173504.3 3375543, 173541.2 3375509, 173563.5 3375489, 173601.2 3375454, 173623.4 3375434, 173644.9 3375415, 173684.7 3375381, 173706.7 3375362, 173747.3 3375328, 173769.3 3375310, 173811 3375276, 173832.8 3375259, 173875.2 3375226, 173896.9 3375209, 173918.5 3375192, 173962 3375159, 173983.6 3375142, 174027.1 3375109, 174048.6 3375093, 174092.1 3375060, 174113.7 3375043, 174157.2 3375010, 174178.8 3374994, 174222.2 3374960, 174243.7 3374944, 174287.1 3374911, 174308.7 3374894, 174330.3 3374878, 174373.6 3374844, 174395.2 3374828, 174438.9 3374795, 174460.4 3374778, 174504.2 3374745, 174525.7 3374729, 174569.4 3374696, 174590.9 3374680, 174634.5 3374646, 174656 3374630, 174700.4 3374597, 174722 3374581, 174742.4 3374567, 174790.4 3374536, 174860.9 3374494, 174914.5 3374464, 174932.6 3374455, 174988.4 3374428, 175005.4 3374420, 175063.6 3374395, 175079.2 3374389, 175138.5 3374366, 175200.4 3374346, 175213.2 3374341, 175276 3374324, 175340.7 3374307, 175405.8 3374292, 175415.9 3374289, 175481 3374274, 175546.1 3374258, 175611.2 3374243, 175621.3 3374240, 175686.4 3374225, 175751.5 3374209, 175761.6 3374207, 175826.3 3374192, 175893.1 3374181, 175960.8 3374172, 176029.4 3374166, 176098.8 3374162, 176168.1 3374161, 176237.6 3374163, 176307.1 3374167, 176375.5 3374175, 176442.9 3374184, 176509.2 3374197, 176517.7 3374198, 176584.2 3374213, 176648.6 3374230, 176712.4 3374248)
Let's say I want m = 5
points. The first point will always be at the beginning of the LINESTRING
(i.e. relative distance from the origin is 0), the last point is at the end of the LINESTRING
(i.e. relative distance = 1), the other three points are located at 0.25, 0.5 and 0.75 (relative distance from the origin). Moreover, 5 points imply 4 segments.
first_segment <- st_linesubstring(my_linestring, 0, 0.25)
second_segment <- st_linesubstring(my_linestring, 0.25, 0.5)
third_segment <- st_linesubstring(my_linestring, 0.5, 0.75)
fourth_segment <- st_linesubstring(my_linestring, 0.75, 1)
# Bind the results
my_segments <- st_sfc(first_segment, second_segment, third_segment, fourth_segment)
# extract the endpoints
my_points <- st_boundary(my_segments)
# plot
par(mar = rep(0, 4))
plot(my_linestring, reset = FALSE)
plot(my_points, add = TRUE, pch = 16, cex = 2)
Created on 2020-11-15 by the reprex package (v0.3.0)
The same ideas could be readapted for other values of m
using (for example) a for loop (it depends on your problem). Moreover, you should define the LINESTRING
object setting appropriate values for the CRS
and the precision.
EDIT - Extract Coordinates
# Extract the coordinates of each point
st_coordinates(my_points)
#> X Y L1
#> [1,] 172868.6 3376152 1
#> [2,] 173688.2 3375378 1
#> [3,] 173688.2 3375378 2
#> [4,] 174579.8 3374688 2
#> [5,] 174579.8 3374688 3
#> [6,] 175602.0 3374245 3
#> [7,] 175602.0 3374245 4
#> [8,] 176712.4 3374248 4
# Clearly there are some duplicates (i.e. the first point of the second segment
# is the last point of the first segment and so on). You can remove the
# duplicates as follows:
st_coordinates(st_sfc(unique(st_cast(my_points, "POINT"))))
#> X Y
#> 1 172868.6 3376152
#> 2 173688.2 3375378
#> 3 174579.8 3374688
#> 4 175602.0 3374245
#> 5 176712.4 3374248
Created on 2020-11-16 by the reprex package (v0.3.0)
Please check here and here to understand the key differences between sf and sp approaches.
Correct answer by agila on January 19, 2021
Using rgeos::gInterpolate
(e.g., for m = 5):
library("rgeos")
pts <- gInterpolate(Spt1, seq(0, 1, length.out = 5), normalized = TRUE)
pts
# SpatialPoints:
# x y
# 1 172869 3376152
# 2 173688 3375378
# 3 174580 3374688
# 4 175602 3374245
# 5 176712 3374248
# Coordinate Reference System (CRS) arguments: NA
Plot
plot(Spt1)
plot(pts, add = TRUE)
Answered by rcs on January 19, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP