4

I have a shape that I upload form the internet. Based on its calculated centroid I would like to plot from it a 50-degree line and find the coordinates that he intersect with the outline. Any idea how can I do it?

Thanks.

Script:

library(ggplot2)

df = read.table("E:/cloud1.txt")  #stored at https://ufile.io/zq679
colnames(df)<- c("x","y")

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

ggplot(df, aes(x, y))+ geom_point()+ geom_point(aes(meanX, meanY),colour="green",size=2)
lbusett
  • 5,801
  • 2
  • 24
  • 47
elyraz
  • 473
  • 5
  • 18

1 Answers1

6

Here is a solution with sf :

df <- read.table("~/Bureau/cloud1.txt")  #stored at https://ufile.io/zq679
colnames(df) <- c("x","y")
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)
library(sf)
poly <- st_sf(st_sfc(st_polygon(list(as.matrix(df)))))

# Choose the angle (in degrees)
angle <- 50

# 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
# You might need to install the latest github version :
# devtools::install_github("tidyverse/ggplot2")
library(ggplot2)
ggplot() + geom_sf(data=poly, fill = NA) + 
    geom_sf(data=line, color = "gray80", lwd = 3) + 
    geom_sf(data=intersect_line, 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()

enter image description here


Edit 2018-07-03

The dataset from the original question is now gone.
Here is a fully reproducible example with a heart shape and without using geom_sf for the graph to answer a question in the comments.

# 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)


# Transform your data.frame in a sf polygon (the first and last points
# must have the same coordinates)
library(sf)
#> 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)
angle <- 50

# 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]

library(ggplot2)
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()

Created on 2018-07-03 by the reprex package (v0.2.0).

Gilles San Martin
  • 4,224
  • 1
  • 18
  • 31
  • I have updated all my packages (I use R 3.3.3) and i get this error "Error: could not find function "geom_sf"". Did Imissed something? any idea where is the problem? – elyraz Jan 21 '18 at 19:21
  • Yes, you need the latest version of ggplot2 from github. You can install it with the following command : `devtools::install_github("tidyverse/ggplot2")`. I have updated the answer. – Gilles San Martin Jan 21 '18 at 21:56
  • A very nice and cool code. I was trying to change the angles and check the intersections with all the 360 degree possibilities. Unfortunately, it does not work, any idea how to solve it? – wrek Jan 22 '18 at 06:17
  • You are right ! It works only for angles between 0 and 90°. The reason is probably that my basic trigonometry is to far away... The updated version should work for any angle... – Gilles San Martin Jan 22 '18 at 11:09
  • How can I remove the grid and change the color of the outline? – hs100 Jan 23 '18 at 16:38
  • To change the color of the outline you just need to specify the color in `geom_sf` e.g. `geom_sf(data=poly, fill = NA, color = "lightblue")`. The command to remove the graticule is a little more criptic. You need to add this to the ggplot graph : `coord_sf(datum = NA)` – Gilles San Martin Jan 23 '18 at 18:46
  • I get this error with another figure "Error in CPL_geos_op2(op, st_geometry(x), st_geometry(y)) : Evaluation error: TopologyException: Input geom 0 is invalid: Self-intersection at or near point 908 913 at 908 913.". What does it mean? and how can I solve it? – elyraz Jan 28 '18 at 15:48
  • Difficult to help without a reproducible example. Your figure is probably not an simple open figure (the lines are crossing each-other). You can visualize it with a simple plot. – Gilles San Martin Jan 28 '18 at 17:53
  • Is there a way to bypass the geom_sf in the plot section? – hs100 Jul 02 '18 at 18:33
  • @user2916044 : yes you can bypass the geom_sf. See the updated answer – Gilles San Martin Jul 02 '18 at 22:56
  • Thanks a lot for the rapid response. – elyraz Jul 03 '18 at 06:13