1

Based on this answer How to get the coordinates of an intesected line with an outline - R ,I tried to run a loop using the script below. Any idea why I can not plot all the intersection points and lines? The shape is different than the answer given

Code:

library(ggplot2)
library(sf)
  
t <- seq(0, 2*pi, by=0.1)
df <- data.frame(x = 13*sin(t)^3,
                 y = 4*cos(t)-2*cos(3*t)-5*cos(4*t)-cos(2*t))
df <- rbind(df, df[1,]) # close the polygon
 
meanX <- mean(df$x)
meanY <- mean(df$y)
 
# Transform your data.frame in a sf polygon (the first and last points
# must have the same coordinates)
 
#> Linking to GEOS 3.5.1, GDAL 2.1.3, proj.4 4.9.2
poly <- st_sf(st_sfc(st_polygon(list(as.matrix(df)))))

# Choose the angle (in degrees)
rotAngles <- 5

for(angle in seq(0,359,rotAngles)) {
  
# Find the minimum length for the line segment to be always 
# outside the cloud whatever the choosen angle
maxX <- max(abs(abs(df[,"x"]) - abs(meanX)))
maxY <- max(abs(abs(df[,"y"]) - abs(meanY)))
line_length = sqrt(maxX^2 + maxY^2) + 1

# Find the coordinates of the 2 points to draw a line with 
# the intended angle.
# This is the gray line on the graph below
line <- rbind(c(meanX,meanY), 
              c(meanX + line_length * cos((pi/180)*angle), 
                meanY + line_length * sin((pi/180)*angle)))
# Transform into a sf line object
line <- st_sf(st_sfc(st_linestring(line)))

# Intersect the polygon and line. The result is a two points line
# shown in black on the plot below
intersect_line <- st_intersection(poly, line)

# Extract only the second point of this line.
# This is the intersecting point
intersect_point <- st_coordinates(intersect_line)[2,c("X","Y")] 

# Visualise this with ggplot and without geom_sf

# you need first transform back the lines into data.frame
line <- as.data.frame(st_coordinates(line))[,1:2]
intersect_line <- as.data.frame(st_coordinates(intersect_line))[,1:2]


ggplot() + geom_path(data=df, aes(x = x, y = y)) + 
  geom_line(data=line, aes(x = X, y = Y), color = "gray80", lwd = 3) + 
  geom_line(data=intersect_line, aes(x = X, y = Y), color = "gray20", lwd = 1) + 
  geom_point(aes(meanX, meanY), colour="orangered", size=2) +
  geom_point(aes(intersect_point["X"], intersect_point["Y"]), 
             colour="orangered", size=2) +
  theme_bw()
}
JohnSal
  • 378
  • 3
  • 11
  • What did you expect to happen wrapping the one polygon, one line code into a one polygon, 71 line loop? You've made no provision for tracking iterations LHS <- RHS, which generally would look like LHS[i] <- RHS[i], all the way down through plotting, unless you just want the last. – Chris Jul 27 '21 at 20:46
  • Thanks but where shell I put the [angle]. Furthermore, can I store all the intersection points in a dataframe? – JohnSal Jul 28 '21 at 12:24

1 Answers1

1

First we'll go back to @Gilles polygon shape as it is more consistent with his reasoning and presentation:

# Generate a heart shape

t <- seq(0, 2*pi, by=0.1)
df <- data.frame(x = 16*sin(t)^3, 
                 y = 13*cos(t)-5*cos(2*t)-2*cos(3*t)-cos(4*t))
df <- rbind(df, df[1,]) # close the polygon

meanX <- mean(df$x)
meanY <- mean(df$y)

library(sf)
poly <- st_sf(st_sfc(st_polygon(list(as.matrix(df)))))

These elements don't change and don't need to calculated multiple times inside a loop:

maxX <- max(abs(abs(df[,"x"]) - abs(meanX)))
maxY <- max(abs(abs(df[,"y"]) - abs(meanY)))
line_length = sqrt(maxX^2 + maxY^2) + 1

Then your rotAngle and angle:

rotAngle <- 5
angle <- seq(0, 359, rotAngle)

Focusing attention on the for loop, the first line call has elements that do and don't change. Let's make an empty list to hold our results, made outside the for loop, that will hold 2x2 matrices:

line_lst <- list()

for (j in 1:length(angle)) {
line_lst[[j]] <- matrix(nrow = 2, ncol=2)
line_lst[[j]][1,1] <- meanX
line_lst[[j]][1,2] <- meanY
line_lst[[j]][2,1] <- meanX + line_length * cos((pi/180)*angle[j])
 line_lst[[j]][2,2] <- meanY + line_length * sin((pi/180)*angle[j])
 }

line_lst[[1]]
             [,1]       [,2]
[1,] 1.225402e-06 0.09131118
[2,] 2.425684e+01 0.09131118

line_lst[[72]]
             [,1]        [,2]
[1,] 1.225402e-06  0.09131118
[2,] 2.416454e+01 -2.02281169

Those seem reasonable, and this was mainly what I wanted to show, explicating on the LHS[j] <- RHS[j], with the which iteration we're on in 1:length(angle). And on to linestring, intersection, and points, same make an empty receiver, loop thru:

# here we have mismatch of establishing an `i` counter then
# counting `j`, which look close enough to tired eyes
# this will result in NULL(s)
linestring_lst <- list()
for (i in 1:length(line_lst)) { # this causes future error
linestring_lst[[j]] <- st_sf(st_sfc(st_linestring(line_lst[[j]]))) 
}

# simply keeping our accounting right, using all `i` or all `j`,
# or staying away from things that look alike and using `k` here

for (k in 1:length(line_lst)) { 
linestring_lst[[k]] <- st_sf(st_sfc(st_linestring(line_lst[[k]]))) 
}

intersection_lst <- list()
for (j in 1:length(linestring_lst)) {
intersection_lst[[j]] <- st_intersection(poly, linestring_lst[[j]])
}

intersect_points <- list()
for (j in 1:length(intersection_lst)) {
intersect_points[[j]] <- st_coordinates(intersection_lst[[j]])[2,c('X','Y')]
}

The things to remember here related to for loops, create your receiver objects outside the loop, index both the LHS[j] and RHS[j] ([ for vector-like receivers, [[ for lists). And having done each of these independently, you can put it all in one for loop.

And final step, take the lists to data.frame(s) for use in ggplot.

intersect_pts_df <- as.data.frame(do.call('rbind', intersect_points))
head(intersect_pts_df, n = 3)
         X          Y
1 14.96993 0.09131118
2 15.56797 1.45333163
3 15.87039 2.88968964
Chris
  • 1,647
  • 1
  • 18
  • 25
  • How do you perform the intersection points between the loop lines *(meanX, meanY, meanX + line_length * cos((pi/180)*angle[j]), meanY + line_length * sin((pi/180)*angle[j]))* and the heart shape. – jane Jul 28 '21 at 20:21
  • Sketched out above. – Chris Jul 29 '21 at 00:36
  • Any idea why I get this error: Error in UseMethod("st_intersection") : no applicable method for 'st_intersection' applied to an object of class "function"? – jane Jul 29 '21 at 18:09
  • What is returned from `class(poly)` and `class(intersection_lst)`, but something is a function rather than the result of a function that st_intersection can work with. Where'd it come up for you? And you can get the highlighting you want by putting phrases or words between '`' backticks. – Chris Jul 29 '21 at 20:26
  • class(poly) is a "function" and class(intersection_lst) is a "list". I do not see that you declared poly in your script – jane Jul 29 '21 at 20:58
  • You're right, thanks for checking, forgot to bring that across from @johnsal, and picked up as function from attached pkgs. Updated. – Chris Jul 29 '21 at 22:24
  • Thanks for the help. I get this error: > for (j in 1:length(linestring_lst)) { + intersection_lst[[j]] <- st_intersection(poly, linestring_lst[[j]]) + } *Error in UseMethod("st_geometry") : no applicable method for 'st_geometry' applied to an object of class "NULL"* – JohnSal Jul 30 '21 at 04:07
  • I got that Error in UseMethod("st_geometry") when I had a mismatch of for (j in 1:length(linestring_lst)) { + intersection_lst[[i]], 'j' outside and 'i' inside the loop in my case, though you appear to have j(s) throughout? No, that was the prior, linestring_lst creation that was NULL for the mismatch, then returned Error in st_intersection due to NULL input. Yes, went down that rabbit hole too! – Chris Jul 30 '21 at 12:22
  • I still do not get it. As rotAngle <- 5 there should always be intersection points – JohnSal Jul 31 '21 at 19:08
  • Edited with comments where my bad copy/pasting led you astray. Appears to work now... – Chris Jul 31 '21 at 20:31