8

BACKGROUND:

A similar question has been asked before, but not answered; see Creating ellipse around line in 3D space using PostGIS. I am a new user asking this question but hopefully with more details. I am more of an R programmer and not a geospatial expert, so my terminology may seem slightly off!

PROBLEM:

I have a line in 3-dimensional space connecting two antenna points above the earth's surface (actually this is a simplification because I have a stationary receiver and many transmitter locations).

What I need to do is use R to automate the creation of 3D ellipses around each of these 'line of sight' lines. The ellipse represents a fresnel zone that must be inspected to assess blockage due to terrain. Ideally, the output would be a spatial file that can be imported into GIS (.shp) and/or GoogleEarth (.kml).

While I know this is possible--there is a website (www.radiofresnel.com) that generates exactly what I need (a 3D Fresnel Zone with explicit geographic coordinates, provided as a .kml file that is easy to import into Google Earth and assess obstructions)--I have reached a certain point and now I am stuck. I think it should be possible to generate an analogous feature in R.

WHAT I HAVE DONE:

Here is an illustration of what I am trying to achieve: enter image description here

I have successfully implemented steps 0 - 5, and am still struggling with steps 6-8.

REPRODUCIBLE EXAMPLE:

Here is my database, representing the outcome of steps 0 through 5. Each row is one of 13 evenly spaced points at approx. ~1000 m intervals along the 'line-of-sight' line between a transmitter and receiver position. Each point has 3D (X,Y,Z) coordinates. For each of these points, there is also the radius of the short axis of the ellipse at each of these center points.

# LIBRARIES
library(scatterplot3d)
library(rgl)
library(sf)
library(tidyverse)
library(matlib)

# OUTCOME OF STEPS 1-5
x <- c(529931.0, 529091.2, 528251.3, 527411.5, 526571.7, 525731.8, 524892.0, 524052.2, 523212.3, 522372.5, 521532.7, 520692.8, 519853.0)
y <- c(3497607, 3498068, 3498529, 3498990, 3499451, 3499912, 3500374, 3500835, 3501296, 3501757, 3502218, 3502679, 3503140)
z <- c(1602.000, 1612.083, 1622.167, 1632.250, 1642.333, 1652.417, 1662.500, 1672.583, 1682.667, 1692.750, 1702.833, 1712.917, 1723.000)
radius <- c(0.00000, 39.94727, 53.86488, 62.58534, 68.13429, 71.25654, 72.26733, 71.25654, 68.13429, 62.58534, 53.86488, 39.94727,  0.00000)

df <- cbind.data.frame(x,y,z,radius)

Where I am stuck: I only can partially achieve the remaining steps. Specifically, I have only been able to visualize the desired outcome in rgl, but visualization in R/RStudio is not my end goal. I need this shape to exist in some valid geo-spatial data structure, such as a 3D polygon that can be imported in a GIS and/or Google Earth workflow.

Step 6 - render circles using the center point coordinates (xi,yi,zi) and the point specific radius (radiusi)

for (i in 1:length(df$z)) {
  circle3d(c(df$z[i], df$y[i], df$x[i]), df$radius[i])
}

enter image description here

Step 7 - add a series of equally spaced lines (n) perpendicular to the outer edge of the circles

shade3d(turn3d(x = x, y = radius, n=100), col = "green")

enter image description here

Step 8 - output shape as a geo-spatial data structure, such as a 3D polygon, that can be used in a GIS workflow.

# I need to figure out Steps 6 & 7 before I can figure this out. 

UPDATE #1: I think have been able to achieve Step 6 by creating a point cloud that delineates the circumference of the circles with centers x[i], y[i], z[i], and radius[i] along the long axis of the Fresnel Zone. There is probably a better way, but below is the code I wrote to achieve Step 6. I am still trying to figure out how to achieve Steps 7 and 8.

# CUSTOM FUNCTIONS
## Circle Point Function - creates a point cloud that defines the outer boundary of the Fresnel Zone at each point along the long axis 
## CREDIT: Adapted from circle3d function in the matlib v0.9.5 package (Friendly 2021)
circle <- function (center, radius, segs, fill = FALSE, ...) {
   angles <- seq(0, 2 * pi, length = segs)
   x <- center[1] + radius * sin(angles) # x-coordinates of each point along the circumference
   y <- center[2] + radius * cos(angles) # y-coordinates of each point along the circumference
}

# SETUP CIRCLE FRAMEWORK
## Set the number of points wanted to delineate the circumference of each circle (higher number equals higher resolution)
## this is done in terms of quarter arc due to the way the 'circle' function orders coordinates, starting with the max
## this approach is necessary to match up y coordinates (max at 270 degrees, min at 90) to proper z coordinates (max at 0/360 degrees, min at 180) for generating the circles
quarterArc.steps <- 20 # setting this to 20 will create 80 + 1 points (20 per quarter arc)
segs <- (quarterArc.steps*4) + 1
quarterArc <- (segs - 1) /  4

#GENERATE POINTS ALONG THE CIRCUMFERENCE OF EACH CIRCLE
datalist <- c()
for(i in 1:length(df$x)) {
  # Creates y and z coordinates at a specific point 'x' along the long axis, order is ymax to ymax and zmax to zmax 
  dat.circle <- data.frame(
    pts.y = c(circle(c(df$z[i], df$y[i], df$x[i]), df$radius[i], segs = segs)), 
    pts.z = c(circle(c(df$z[i], df$y[i], df$x[i]), df$radius[i], segs = segs)), 
    pts.x = rep(df$x[i], length(seq(from = 1, to = segs))))
  # Allign the y and z coordinates 
  index <- c(seq(from = (quarterArc + 1), to = segs), seq(from = 1, to = quarterArc))
  dat.circleNew <- data.frame(
    Y = dat.circle$pts.y[order(index)],
    Z = dat.circle$pts.z,
    X = dat.circle$pts.x)
  dat.circleNew$CircleID <- i
  datalist[[i]] <- dat.circleNew
}

# STORE OUTCOME AS A DATAFRAME
df.Step6 = do.call(rbind, datalist)
## Convert to a SFC_POINT Geometry
sfc_circumference <- df.Step6 %>%
  select(X, Y, Z) %>% 
  st_as_sf(coords = c("X", "Y", "Z"))
## Visualize
scatterplot3d(df.Step6$X, df.Step6$Y, df.Step6$Z, pch=20)

enter image description here

8
  • What file format is acceptable for your 3d polygon output? What can your desired GIS import?
    – Spacedman
    Jan 6 at 12:20
  • I primarily work with .shp files, but I get the impression that writing a 3D object to .shp may not be possible. I did try outputting what I have so far with st_write(sfc_circumference, "FresZone_try1.shp") and I get this error: "GDAL Error 6: Geometry type of `3D Point' not supported in shapefiles. Type can be overridden with a layer creation option of SHPT=POINT/ARC/POLYGON/MULTIPOINT/POINTZ/ARCZ/POLYGONZ/MULTIPOINTZ/MULTIPATCH." It appears R may read, but not write, .gdb files. ESRI ARC GIS can import .kml and convert to .gdb. Any recommendations? I can post as a separate inquiry.
    – Wu Wei
    Jan 6 at 16:43
  • Your circles that you generated in step 6 will always be perpendicular to the x axis, i.e., the circles will be vertical. Are you sure you want this? Do you not want the circles to be perpendicular to the axis of the ellipsoid that is the fresnel zone?
    – joy
    Jan 10 at 10:40
  • Also, cross sections of an ellipsoid which are not perpendicular to the major axis of the ellipsoid will not be circular. Rather, those cross sections will be ellipses.
    – joy
    Jan 10 at 10:41
  • 1
    @joy, I used 'ScreenToGIF' (screentogif.com) to capture the RGL output on my screen
    – Wu Wei
    Jan 16 at 17:10

2 Answers 2

1
+100

The question asks to solve a number of steps in a bigger problem, which is to produce the Fresnel ellipsoid given the latitudes, longitudes and altitudes of the two end points. This answer considers the original problem as a whole. The following R function fresnelellipsoid_kml creates a kml file of the Fresnel ellipsoid between the two endpoints. The function can also be used to construct any ellipsoid between two endpoints above the surface of the earth.

An ellipsoid is completely determined by the two endpoints, i.e., the poles, and the maximum radius of a cross-sectional circle. The maximum radius of the nth Fresnel zone corresponding to frequency f is given by (cnD/4f)1/2, where c is the speed of light in vacuum (299792458 m) and D is the distance between the two endpoints. The ellipsoid is constructed by considering a number of its cross sectional circles and then appropriately joining equidistant points on the perimeters of the circles. The WGS 84 coordinate system is used for computation, which considers the earth as an oblate spheroid with equatorial radius a = 6378137 m and flattening f = 1/298.257223563.

The arguments of fresnelellipsoid_kml are the following:

  • latlongalt1: a numeric vector with three components: latitude in degrees (southern latitude is negative), longitude in degrees (western longitude is negative) and altitude in meters above the mean sea level for the first endpoint.
  • latlongalt2: a numeric vector with three components: latitude in degrees (southern latitude is negative), longitude in degrees (western longitude is negative) and altitude in meters above the mean sea level for the second endpoint.
  • max_radius: maximum radius (in meters) of a cross-sectional circle of the ellipsoid between the endpoints. This value is used if fresnel = FALSE (see below). Default value is 100.
  • segments: number of line segments used to draw the perimeter of a cross-sectional circle of the ellipsoid between the endpoints. Minimum value is 10, and default value is 100.
  • sections: number of cross-sectional circles used to construct the three dimensional polygon approximating the ellipsoid between the endpoints. Minimum value is 10, and default value is 100.
  • fresnel: logical argument, indicates whether to draw a Fresnel ellipsoid (TRUE) or a custom ellipsoid with a given value of max_radius (FALSE). Default is FALSE.
  • fresnel_freq: the frequency corresponding to the Fresnel zone in Hz. Required if fresnel = TRUE.
  • fresnel_num: The Frensel number corresponding to the Fresnel zone. Required if fresnel = TRUE.
  • filename: the name of the kml file to be written in the current working directory.
  • overwrite: whether to overwrite a kml file present in the current directory has the same name as filename. Default is TRUE.
  • ...: optional arguments used in the preparation of the kml file. These are described below.
    • indentsymbol: the indent character while writing the kml file, either '\t' (tab) or ' ' (double spaces). Default taken is '\t'.
    • lineofsight_color: a numeric vector of four components specifying the color of the straight line joining the two endpoints. The four components are between 0 and 1, representing opacity (0 is transparent and 1 is completely opaque), and the intensities of blue, green and red. Default values are opacity = 0.7, blue = 1, green = 0 and red = 0.
    • lineofsight_width: the width (in pixels) of the straight line joining the two endpoints. Default value taken is 5.
    • polygon_color: a numeric vector of four components specifying the color used for drawing the polygonal faces of the three dimensional polygon approximating the ellipsoid. The four components are between 0 and 1, representing opacity, and the intensities of blue, green and red. Default values are opacity = 0.3, blue = 0, green = 1 and red = 1, i.e., the default style of the three dimensional polygon is translucent yellow.
    • p1_name: name of the first endpoint. Default is Point 1.
    • p2_name: name of the second endpoint. Default is Point 2.
fresnelellipsoid_kml = function(latlongalt1, latlongalt2, max_radius = 100, segments = 100, sections = 100,
                                fresnel = FALSE, fresnel_freq, fresnel_num,
                                filename, overwrite = TRUE, ...){
  ##### Functions for transformations between coordinate systems #####
  
  latlong2spherical = function(latlongalt){
    lat = latlongalt[1]
    long = latlongalt[2]
    altitude = latlongalt[3]
    
    lat = (lat / 180) * pi
    long = (long / 180) * pi
    
    a = 6378137
    b = a
    
    f = 1 / 298.257223563
    c = a * (1 - f)
    
    phi = pi - long
    theta = (pi / 2) - lat
    
    r = sqrt((a * cos(phi) * sin(theta))^2 + (b * sin(phi) * sin(theta))^2 + (c * cos(theta))^2) + altitude
    
    spherical = c(r, phi, theta)
    
    return(spherical)
  }
  
  spherical2latlong = function(spherical){
    r = spherical[1]
    phi = spherical[2]
    theta = spherical[3]
    
    a = 6378137
    b = a
    
    f = 1 / 298.257223563
    c = a * (1 - f)
    
    altitude = r - sqrt((a * cos(phi) * sin(theta))^2 + (b * sin(phi) * sin(theta))^2 + (c * cos(theta))^2)
    
    long = pi - phi
    lat = (pi / 2) - theta
    
    lat = (lat / pi) * 180
    long = (long / pi) * 180
    
    latlongalt = c(lat, long, altitude)
    
    return(latlongalt)
  }
  
  cartesian2spherical = function(cartesian){
    x = cartesian[1]
    y = cartesian[2]
    z = cartesian[3]
    
    r = sqrt(sum(cartesian^2))
    
    theta = acos(z / r)
    
    phi = atan2(y, x)
    if (phi < 0)
      phi = 2*pi + phi
    
    spherical = c(r, phi, theta)
    
    return(spherical)
  }
  
  spherical2cartesian = function(spherical){
    r = spherical[1]
    phi = spherical[2]
    theta = spherical[3]
    
    x = r * cos(phi) * sin(theta)
    y = r * sin(phi) * sin(theta)
    z = r * cos(theta)
    
    cartesian = c(x, y, z)
    
    return(cartesian)
  }
  
  ##### Checking the function arguments required for the computation of the polygon approximating the ellipsoid #####
  
  if (!is.numeric(latlongalt1) || length(latlongalt1) != 3)
    stop('The latitude, longitude and altitude (above the mean sea level) of the first point needs to be given
         as a numeric vector with three components: latitude in degrees (Southern latitude is negative), longitude 
         in degrees (Western longitude is negative) and altitude in meters above the mean sea level.')
  if (!is.numeric(latlongalt2) || length(latlongalt2) != 3)
    stop('The latitude, longitude and altitude (above the mean sea level) of the second point needs to be given
         as a numeric vector with three components: latitude in degrees (Southern latitude is negative), longitude 
         in degrees (Western longitude is negative) and altitude in meters above the mean sea level.')
  
  if (latlongalt1[1] < -90 || latlongalt1[1] > 90)
    stop('The latitude of the first point must be between -90 and 90 degrees.')
  if (latlongalt1[2] < -180 || latlongalt1[2] > 180)
    stop('The longitude of the first point must be between -180 and 180 degrees.')
  if (latlongalt1[3] < 0 || latlongalt1[3] == Inf)
    stop('The altitude (in meters) of the first point must be a positive real number.')
  
  if (latlongalt2[1] < -90 || latlongalt2[1] > 90)
    stop('The latitude of the second point must be between -90 and 90 degrees.')
  if (latlongalt2[2] < -180 || latlongalt2[2] > 180)
    stop('The longitude of the second point must be between -180 and 180 degrees.')
  if (latlongalt2[3] < 0 || latlongalt2[3] == Inf)
    stop('The altitude (in meters) of the second point must be a positive real number.')
  
  if (!is.numeric(segments) || length(segments) != 1)
    stop('The argument `segments` mst be a positive integer.')
  if (segments < 10){
    warning('The value of the argument `segments` must be at least 10. `segments = 10` is set.')
    segments = 10
  }
  
  if (!is.numeric(sections) || length(sections) != 1)
    stop('The argument `sections` mst be a positive integer.')
  if (sections < 10){
    warning('The value of the argument `sections` must be at least 10. `sections = 10` is set.')
    sections = 10
  }
  
  if (fresnel == TRUE){
    if (!is.numeric(fresnel_freq) || length(fresnel_freq) > 1)
      stop('Fresnel frequency must be a real number.')
    if (fresnel_freq <= 0 || fresnel_freq == Inf)
      stop('Fresnel frequency must be a positive real number.')
    if (!is.numeric(fresnel_num) || length(fresnel_num) > 1)
      stop('Fresnel number must be a real number.')
    if (fresnel_num <= 0 || fresnel_num == Inf)
      stop('Fresnel number must be a positive real number.')
  }else{
    if (!is.numeric(max_radius) || length(max_radius) > 1)
      stop('`max_radius` must be a real number.')
    if (max_radius <= 0 || max_radius == Inf)
      stop('`max_radius` must be a positive real number.')
  }
  
  ##### Computation of the vertices of the polygon approximating the ellipsoid #####
  
  p1 = spherical2cartesian(latlong2spherical(latlongalt1))
  p2 = spherical2cartesian(latlong2spherical(latlongalt2))
  
  major_axis = sqrt(sum((p1 - p2)^2))
  
  if (fresnel == TRUE){
    speed_of_light = 299792458
    r = sqrt((fresnel_num * (speed_of_light / fresnel_freq) * (major_axis / 2)^2) / major_axis)
  }else{
    r = max_radius
  }
  
  a = major_axis / 2
  b = r
  c = r
  
  x_vector = seq(from = -a, to = a, length.out = sections)
  
  angles = seq(from = 0, to = 2*pi, length.out = (segments + 1))
  
  radius_vector = r * sqrt(1 - (x_vector^2 / a^2))
  
  y_matrix = matrix(radius_vector, nrow = length(radius_vector), ncol = 1) %*% matrix(cos(angles), nrow = 1, ncol = length(angles))
  
  z_matrix = matrix(radius_vector, nrow = length(radius_vector), ncol = 1) %*% matrix(sin(angles), nrow = 1, ncol = length(angles))
  
  vertices = matrix(nrow = length(radius_vector) * length(angles), ncol = 3)
  index_record_vertices = matrix(nrow = length(radius_vector), ncol = length(angles))
  count = 0;
  for (i in 1:length(radius_vector)){
    for (j in 1:length(angles)){
      count = count + 1;
      vertices[count,] = c(x_vector[i], y_matrix[i,j], z_matrix[i,j])
      
      index_record_vertices[i, j] = count
    }
  }
  
  center = (p1 + p2) / 2
  p2_shifted = p2 - center
  
  roll_angle = 0
  pitch_angle = - asin(p2_shifted[3] / sqrt(sum(p2_shifted^2)))
  yaw_angle = atan2(p2_shifted[2], p2_shifted[1])
  
  Rotation_X = matrix(c(1, 0, 0, 0, cos(roll_angle), -sin(roll_angle), 0, sin(roll_angle), cos(roll_angle)),
                      nrow = 3, ncol = 3, byrow = TRUE)
  Rotation_Y = matrix(c(cos(pitch_angle), 0, sin(pitch_angle), 0, 1, 0, -sin(pitch_angle), 0, cos(pitch_angle)),
                      nrow = 3, ncol = 3, byrow = TRUE)
  Rotation_Z = matrix(c(cos(yaw_angle), -sin(yaw_angle), 0, sin(yaw_angle), cos(yaw_angle), 0, 0, 0, 1),
                      nrow = 3, ncol = 3, byrow = TRUE)
  
  Rotation_matrix = Rotation_Z %*% Rotation_Y %*% Rotation_X
  
  rotated_vertices = vertices %*% t(Rotation_matrix)
  
  final_vertices = rotated_vertices + matrix(center, nrow = nrow(rotated_vertices), ncol = length(center), byrow = TRUE)
  
  spherical_vertices = t(apply(final_vertices, 1, cartesian2spherical))
  
  latlong_vertices = t(apply(spherical_vertices, 1, spherical2latlong))
  
  longlatalt_vertices = cbind(latlong_vertices[,2], latlong_vertices[,1], latlong_vertices[,3])
  
  ##### Creating the kml file #####
  
  if (file.exists(paste(getwd(), paste(filename, 'kml', sep = '.'), sep = '/'))){
    if (overwrite != TRUE){
      stop(paste("'", paste(filename, 'kml', sep = '.'), "'", ' already exists in ', "'", getwd(), "'",
                 '. Check ', "'overwrite' oprion.", sep = ''))
    }else{
      checkremove = file.remove(paste(getwd(), paste(filename, 'kml', sep = '.'), sep = '/'))
      if (checkremove != TRUE)
        stop(paste('Could not remove ', "'", paste(filename, 'kml', sep = '.'), "'", ', please check.', sep = ''))
    }
  }else{
    file.create(paste(getwd(), paste(filename, 'kml', sep = '.'), sep = '/'))
  }
  
  connection = file(paste(getwd(), paste(filename, 'kml', sep = '.'), sep = '/'), open = 'wt')
  
  ##### Reading and fixing the arguments for writing the kml file #####
  
  style_args = list(...)
  
  if (is.null(style_args$indentsymbol) || !(style_args$indentsymbol %in% c('\t', '  '))){
    indentsymbol = '\t'
  }else{
    indentsymbol = style_args$indentsymbol
  }
  
  color_proportion2hex = function(r){
    if (!is.numeric(r))
      stop('`r` must be a number.')
    if (length(r) > 1)
      stop('`r` must be a number, not a vector.')
    if (r < 0 || r > 1)
      stop('`r` must be within 0 and 1.')
    
    r_integer = round(r * 255)
    
    hexdigits = c('0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'a', 'b', 'c', 'd', 'e', 'f')
    hexvalue = paste(hexdigits[floor(r_integer / 16) + 1], hexdigits[r_integer - (16 * floor(r_integer / 16)) + 1], sep = '')
    
    return(hexvalue)
  }
  
  if (is.null(style_args$lineofsight_color) || !is.numeric(style_args$lineofsight_color)){
    lineofsight_color_opacity = 0.7
    lineofsight_color_blue = 1
    lineofsight_color_green = 0
    lineofsight_color_red = 0
    
    lineofsight_color = paste(color_proportion2hex(lineofsight_color_opacity), color_proportion2hex(lineofsight_color_blue),
                              color_proportion2hex(lineofsight_color_green), color_proportion2hex(lineofsight_color_red), sep = '')
  }else{
    if (!is.na(style_args$lineofsight_color[1]) && !(style_args$lineofsight_color[1] >= 0 && style_args$lineofsight_color[1] <= 1)){
      lineofsight_color_opacity = 0.7
    }else{
      lineofsight_color_opacity = style_args$lineofsight_color[1]
    }
    
    if (!is.na(style_args$lineofsight_color[2]) && !(style_args$lineofsight_color[2] >= 0 && style_args$lineofsight_color[2] <= 1)){
      lineofsight_color_blue = 1
    }else{
      lineofsight_color_blue = style_args$lineofsight_color[2]
    }
    
    if (!is.na(style_args$lineofsight_color[3]) && !(style_args$lineofsight_color[3] >= 0 && style_args$lineofsight_color[3] <= 1)){
      lineofsight_color_green = 0
    }else{
      lineofsight_color_green = style_args$lineofsight_color[3]
    }
    
    if (!is.na(style_args$lineofsight_color[4]) && !(style_args$lineofsight_color[4] >= 0 && style_args$lineofsight_color[4] <= 1)){
      lineofsight_color_red = 0
    }else{
      lineofsight_color_red = style_args$lineofsight_color[4]
    }
    
    lineofsight_color = paste(color_proportion2hex(lineofsight_color_opacity), color_proportion2hex(lineofsight_color_blue),
                              color_proportion2hex(lineofsight_color_green), color_proportion2hex(lineofsight_color_red), sep = '')
  }
  
  if (is.null(style_args$lineofsight_width) || !is.numeric(style_args$lineofsight_width)){
    lineofsight_width = 5
  }else{
    if (style_args$lineofsight_width[1] < 0 || style_args$lineofsight_width[1] == Inf){
      lineofsight_width = 5
    }else{
      lineofsight_width = style_args$lineofsight_width[1]
    }
  }
  
  if (is.null(style_args$polygon_color) || !is.numeric(style_args$polygon_color)){
    polygon_color_opacity = 0.3
    polygon_color_blue = 0
    polygon_color_green = 1
    polygon_color_red = 1
    
    polygon_color = paste(color_proportion2hex(polygon_color_opacity), color_proportion2hex(polygon_color_blue),
                          color_proportion2hex(polygon_color_green), color_proportion2hex(polygon_color_red), sep = '')
  }else{
    if (!is.na(style_args$polygon_color[1]) && !(style_args$polygon_color[1] >= 0 && style_args$polygon_color[1] <= 1)){
      polygon_color_opacity = 0.3
    }else{
      polygon_color_opacity = style_args$polygon_color[1]
    }
    
    if (!is.na(style_args$polygon_color[2]) && !(style_args$polygon_color[2] >= 0 && style_args$polygon_color[2] <= 1)){
      polygon_color_blue = 0
    }else{
      polygon_color_blue = style_args$polygon_color[2]
    }
    
    if (!is.na(style_args$polygon_color[3]) && !(style_args$polygon_color[3] >= 0 && style_args$polygon_color[3] <= 1)){
      polygon_color_green = 1
    }else{
      polygon_color_green = style_args$polygon_color[3]
    }
    
    if (!is.na(style_args$polygon_color[4]) && !(style_args$polygon_color[4] >= 0 && style_args$polygon_color[4] <= 1)){
      polygon_color_red = 1
    }else{
      polygon_color_red = style_args$polygon_color[4]
    }
    
    polygon_color = paste(color_proportion2hex(polygon_color_opacity), color_proportion2hex(polygon_color_blue),
                          color_proportion2hex(polygon_color_green), color_proportion2hex(polygon_color_red), sep = '')
  }
  
  if (is.null(style_args$p1_name) || !is.character(style_args$p1_name)){
    p1_name = 'Point 1'
  }else{
    p1_name = style_args$p1_name
  }
  
  if (is.null(style_args$p2_name) || !is.character(style_args$p2_name)){
    p2_name = 'Point 2'
  }else{
    p2_name = style_args$p2_name
  }
  
  ##### Writing text to the kml file #####
  
  indentlevel = 0
  
  write_indented_text = function(text, connection, indentsymbol, indentlevel, indent_increment){
    if (!(indent_increment %in% c(1, 0, -1)))
      stop('`indent_increment` must be either 1 or 0 or -1.')
    if (indentlevel + indent_increment < 0)
      stop('`indentlevel + indent_increment` must be positive.')
    
    indentlevel = indentlevel + indent_increment
    current_indent = paste(rep(indentsymbol, indentlevel), collapse = '')
    writeLines(c(paste(current_indent, text, sep = ''), '\n'), con = connection, sep = '')
    
    return(indentlevel)
  }
  
  writeLines(c('<?xml version="1.0" encoding="UTF-8"?>', '\n', '<kml xmlns="http://www.opengis.net/kml/2.2">', '\n'),
             con = connection, sep = '')
  
  indentlevel = write_indented_text('<Document>', connection, indentsymbol, indentlevel, indent_increment = 0)
  
  indentlevel = write_indented_text('<Style id="EndPoint">', connection, indentsymbol, indentlevel, indent_increment = 1)
  indentlevel = write_indented_text('<IconStyle>', connection, indentsymbol, indentlevel, indent_increment = 1)
  indentlevel = write_indented_text('<Icon></Icon>', connection, indentsymbol, indentlevel, indent_increment = 1)
  indentlevel = write_indented_text('</IconStyle>', connection, indentsymbol, indentlevel, indent_increment = -1)
  indentlevel = write_indented_text('</Style>', connection, indentsymbol, indentlevel, indent_increment = -1)
  
  indentlevel = write_indented_text('<Style id="LineOfSight">', connection, indentsymbol, indentlevel, indent_increment = 0)
  indentlevel = write_indented_text('<LineStyle>', connection, indentsymbol, indentlevel, indent_increment = 1)
  lineofsight_color_line = paste('<color>', lineofsight_color, '</color>', sep = '')
  indentlevel = write_indented_text(lineofsight_color_line, connection, indentsymbol, indentlevel, indent_increment = 1)
  lineofsight_width_line = paste('<width>', lineofsight_width, '</width>', sep = '')
  indentlevel = write_indented_text(lineofsight_width_line, connection, indentsymbol, indentlevel, indent_increment = 0)
  indentlevel = write_indented_text('</LineStyle>', connection, indentsymbol, indentlevel, indent_increment = -1)
  indentlevel = write_indented_text('</Style>', connection, indentsymbol, indentlevel, indent_increment = -1)
  
  indentlevel = write_indented_text('<Style id="PolygonStyle">', connection, indentsymbol, indentlevel, indent_increment = 0)
  indentlevel = write_indented_text('<LineStyle>', connection, indentsymbol, indentlevel, indent_increment = 1)
  polygon_color_line = paste('<color>', polygon_color, '</color>', sep = '')
  indentlevel = write_indented_text(polygon_color_line, connection, indentsymbol, indentlevel, indent_increment = 1)
  indentlevel = write_indented_text('</LineStyle>', connection, indentsymbol, indentlevel, indent_increment = -1)
  indentlevel = write_indented_text('<PolyStyle>', connection, indentsymbol, indentlevel, indent_increment = 0)
  indentlevel = write_indented_text(polygon_color_line, connection, indentsymbol, indentlevel, indent_increment = 1)
  indentlevel = write_indented_text('</PolyStyle>', connection, indentsymbol, indentlevel, indent_increment = -1)
  indentlevel = write_indented_text('</Style>', connection, indentsymbol, indentlevel, indent_increment = -1)
  
  indentlevel = write_indented_text('<Placemark>', connection, indentsymbol, indentlevel, indent_increment = 0)
  p1_name_line = paste('<name>', p1_name, '</name>', sep = '')
  indentlevel = write_indented_text(p1_name_line, connection, indentsymbol, indentlevel, indent_increment = 1)
  endpoint_styleurl_line = paste('<styleUrl>', '#EndPoint', '</styleUrl>', sep = '')
  indentlevel = write_indented_text(endpoint_styleurl_line, connection, indentsymbol, indentlevel, indent_increment = 0)
  indentlevel = write_indented_text('<Point>', connection, indentsymbol, indentlevel, indent_increment = 0)
  endpoint_altitudemode_line = paste('<altitudeMode>', 'absolute', '</altitudeMode>', sep = '')
  indentlevel = write_indented_text(endpoint_altitudemode_line, connection, indentsymbol, indentlevel, indent_increment = 1)
  indentlevel = write_indented_text('<coordinates>', connection, indentsymbol, indentlevel, indent_increment = 0)
  p1_coordinates = paste(c(latlongalt1[2], latlongalt1[1], latlongalt1[3]), collapse = ',')
  indentlevel = write_indented_text(p1_coordinates, connection, indentsymbol, indentlevel, indent_increment = 1)
  indentlevel = write_indented_text('</coordinates>', connection, indentsymbol, indentlevel, indent_increment = -1)
  indentlevel = write_indented_text('</Point>', connection, indentsymbol, indentlevel, indent_increment = -1)
  indentlevel = write_indented_text('</Placemark>', connection, indentsymbol, indentlevel, indent_increment = -1)
  
  indentlevel = write_indented_text('<Placemark>', connection, indentsymbol, indentlevel, indent_increment = 0)
  p2_name_line = paste('<name>', p2_name, '</name>', sep = '')
  indentlevel = write_indented_text(p2_name_line, connection, indentsymbol, indentlevel, indent_increment = 1)
  endpoint_styleurl_line = paste('<styleUrl>', '#EndPoint', '</styleUrl>', sep = '')
  indentlevel = write_indented_text(endpoint_styleurl_line, connection, indentsymbol, indentlevel, indent_increment = 0)
  indentlevel = write_indented_text('<Point>', connection, indentsymbol, indentlevel, indent_increment = 0)
  endpoint_altitudemode_line = paste('<altitudeMode>', 'absolute', '</altitudeMode>', sep = '')
  indentlevel = write_indented_text(endpoint_altitudemode_line, connection, indentsymbol, indentlevel, indent_increment = 1)
  indentlevel = write_indented_text('<coordinates>', connection, indentsymbol, indentlevel, indent_increment = 0)
  p2_coordinates = paste(c(latlongalt2[2], latlongalt2[1], latlongalt2[3]), collapse = ',')
  indentlevel = write_indented_text(p2_coordinates, connection, indentsymbol, indentlevel, indent_increment = 1)
  indentlevel = write_indented_text('</coordinates>', connection, indentsymbol, indentlevel, indent_increment = -1)
  indentlevel = write_indented_text('</Point>', connection, indentsymbol, indentlevel, indent_increment = -1)
  indentlevel = write_indented_text('</Placemark>', connection, indentsymbol, indentlevel, indent_increment = -1)
  
  indentlevel = write_indented_text('<Placemark>', connection, indentsymbol, indentlevel, indent_increment = 0)
  lineofsight_name_line = paste('<name>', 'Line of sight', '</name>', sep = '')
  indentlevel = write_indented_text(lineofsight_name_line, connection, indentsymbol, indentlevel, indent_increment = 1)
  lineofsight_styleurl_line = paste('<styleUrl>', '#LineOfSight', '</styleUrl>', sep = '')
  indentlevel = write_indented_text(lineofsight_styleurl_line, connection, indentsymbol, indentlevel, indent_increment = 0)
  indentlevel = write_indented_text('<LineString>', connection, indentsymbol, indentlevel, indent_increment = 0)
  lineofsight_altitudemode_line = paste('<altitudeMode>', 'absolute', '</altitudeMode>', sep = '')
  indentlevel = write_indented_text(lineofsight_altitudemode_line, connection, indentsymbol, indentlevel, indent_increment = 1)
  indentlevel = write_indented_text('<coordinates>', connection, indentsymbol, indentlevel, indent_increment = 0)
  p1_coordinates = paste(c(latlongalt1[2], latlongalt1[1], latlongalt1[3]), collapse = ',')
  p2_coordinates = paste(c(latlongalt2[2], latlongalt2[1], latlongalt2[3]), collapse = ',')
  indentlevel = write_indented_text(p1_coordinates, connection, indentsymbol, indentlevel, indent_increment = 1)
  indentlevel = write_indented_text(p2_coordinates, connection, indentsymbol, indentlevel, indent_increment = 0)
  indentlevel = write_indented_text('</coordinates>', connection, indentsymbol, indentlevel, indent_increment = -1)
  indentlevel = write_indented_text('</LineString>', connection, indentsymbol, indentlevel, indent_increment = -1)
  indentlevel = write_indented_text('</Placemark>', connection, indentsymbol, indentlevel, indent_increment = -1)
  
  indentlevel = write_indented_text('<Placemark>', connection, indentsymbol, indentlevel, indent_increment = 0)
  zone_name_line = paste('<name>', 'Fresnel zone', '</name>', sep = '')
  indentlevel = write_indented_text(zone_name_line, connection, indentsymbol, indentlevel, indent_increment = 1)
  zone_styleurl_line = paste('<styleUrl>', '#PolygonStyle', '</styleUrl>', sep = '')
  indentlevel = write_indented_text(zone_styleurl_line, connection, indentsymbol, indentlevel, indent_increment = 0)
  
  indentlevel = write_indented_text('<MultiGeometry>', connection, indentsymbol, indentlevel, indent_increment = 0)
  indentlevel = indentlevel + 1
  for (i in 1:(length(radius_vector) - 1)){
    for (j in 1:(length(angles) - 1)){
      indentlevel = write_indented_text('<Polygon>', connection, indentsymbol, indentlevel, indent_increment = 0)
      zone_altitudemode_line = paste('<altitudeMode>', 'absolute', '</altitudeMode>', sep = '')
      indentlevel = write_indented_text(zone_altitudemode_line, connection, indentsymbol, indentlevel, indent_increment = 1)
      indentlevel = write_indented_text('<outerBoundaryIs>', connection, indentsymbol, indentlevel, indent_increment = 0)
      indentlevel = write_indented_text('<LinearRing>', connection, indentsymbol, indentlevel, indent_increment = 1)
      indentlevel = write_indented_text('<coordinates>', connection, indentsymbol, indentlevel, indent_increment = 1)
      
      current_coordinates_1 = paste(longlatalt_vertices[index_record_vertices[i, j],], collapse = ',')
      current_coordinates_2 = paste(longlatalt_vertices[index_record_vertices[i + 1, j],], collapse = ',')
      current_coordinates_3 = paste(longlatalt_vertices[index_record_vertices[i + 1, j + 1],], collapse = ',')
      current_coordinates_4 = paste(longlatalt_vertices[index_record_vertices[i, j + 1],], collapse = ',')
      
      indentlevel = write_indented_text(current_coordinates_1, connection, indentsymbol, indentlevel, indent_increment = 1)
      indentlevel = write_indented_text(current_coordinates_2, connection, indentsymbol, indentlevel, indent_increment = 0)
      indentlevel = write_indented_text(current_coordinates_3, connection, indentsymbol, indentlevel, indent_increment = 0)
      indentlevel = write_indented_text(current_coordinates_4, connection, indentsymbol, indentlevel, indent_increment = 0)
      indentlevel = write_indented_text(current_coordinates_1, connection, indentsymbol, indentlevel, indent_increment = 0)
      
      indentlevel = write_indented_text('</coordinates>', connection, indentsymbol, indentlevel, indent_increment = -1)
      indentlevel = write_indented_text('</LinearRing>', connection, indentsymbol, indentlevel, indent_increment = -1)
      indentlevel = write_indented_text('</outerBoundaryIs>', connection, indentsymbol, indentlevel, indent_increment = -1)
      indentlevel = write_indented_text('</Polygon>', connection, indentsymbol, indentlevel, indent_increment = -1)
    }
  }
  indentlevel = write_indented_text('</MultiGeometry>', connection, indentsymbol, indentlevel, indent_increment = -1)
  indentlevel = write_indented_text('</Placemark>', connection, indentsymbol, indentlevel, indent_increment = -1)
  
  indentlevel = write_indented_text('</Document>', connection, indentsymbol, indentlevel, indent_increment = -1)
  indentlevel = write_indented_text('</kml>', connection, indentsymbol, indentlevel, indent_increment = 0)
  
  close(connection)
  
  return(longlatalt_vertices)
}

Apart from creating the kml file, the function also returns a matrix with three columns, whose rows contain the longitudes, the latitudes and the altitudes of the vertices of the three dimensional polygon approximating the ellipsoid.

8
  • Thanks @Joy! I can't 100% confirm that this solution works since I need to import to kml by hand (I was hoping to automate this process). That said, your explanation is superb, the code runs without error, and the matrix appears to contain the needed information. It seems you are much better at geometry than I am! One note, I did have to alter a tiny section of the customellipsoid function, specifically chanting segments = 20 and sections = 20 to segments = segments and sections = sections so that user inputs for these two arguments are not overwritten.
    – Wu Wei
    Jan 10 at 19:58
  • Thanks for pointing out the mistake! I have corrected it. The two lines segments = 20 and sections = 20 should not be there. I am looking into automatically form the kml file, and I shall be updating this answer with the process to transfer the lat-long and altitude information to a kml format.
    – joy
    Jan 10 at 20:25
  • @WuWei I was curious about one thing. How do you get the maximum radius of the circle? Is there a formula for the Fresnel zone based on some other input?
    – joy
    Jan 10 at 20:26
  • 1
    I was trying to make my OP as simple as possible so left out the Freznal Zone (FZ) formulas. The radius of First FZ at any given point along the long axis is equal to sqrt((1*d1*d2*(C/freq))/(d1+d2)) where d1 is distance (m) from the start of the FZ to a point, d2 is distance from that same point to the end of the FZ, C is the speed of light (299,792,458 m/s), and freq is the frequency of the transmitter (mine is 165,000,000 Hz). The radius is maximized when d1 equals d2 (the midpoint). See this site for more information pagerpower.com/news/fresnel-zone
    – Wu Wei
    Jan 10 at 21:51
  • Thank you @WuWei, that is informative.
    – joy
    Jan 10 at 21:54
0

0) Intro

Ok so first a disclaimer, I don't know R. But the solution I provide is not code related. Moreover, I think your question is (not a bad one but) a "nested" one : answering it completely implies to ask and answer other sub-questions and write a lot of sub-functions to come up with a template, so I will try to give you the major steps I would do.

Second, I do not know why you are trying to achieve this : Is is to manually check for obstacle in Google Earth or do you plan to create an algorithm to automatically find obstacles (which will not be easy as you will need a complet 3D model of the area you are covering). The final goal of the process may influence the way you solve it.

1) General idea

To display a 3D object in Google Earth, you must build a COLLADA file (.dae for the extension). In reality, a collada file is just a particular kind of XML, so you can edit it exactly as KML file using any language with a library handling XML.

The principle is:

  1. Handle the geometry of your capsule (3D ellipsoid) in the .dae file (regardless of projection).
  2. handle the geographic location and altitude in a kml file and link it to the dae

2) Example template

here is a functionning example that put a 20X20m square on top of the Eiffel Tower. By tweeking only a few lines, you can adapt it to your need, I explain this in the next part so that it is not overwhelming

  1. create a file "theKmlFile.kml" containing the following (notice the href to .dae) :
<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://earth.google.com/kml/2.1">
<Placemark>
    <name>name of the project (not important)</name>
    <description>put a description here if you want (not important)</description>
    <LookAt>
        <longitude>2.2945083333333334</longitude>
        <latitude>48.858272222222226</latitude>
        <altitude>300</altitude>
        <range>745</range>
        <tilt>23</tilt>
        <heading>53</heading>
    </LookAt>
    <Model id="anyId_notimportantHere">
        <altitudeMode>relativeToGround</altitudeMode>
        <Location>
            <longitude>2.2945083333333334</longitude>
            <latitude>48.858272222222226</latitude>
            <altitude>300</altitude>
        </Location>
        <Orientation>
            <heading>0</heading>
            <tilt>0</tilt>
            <roll>0</roll>
        </Orientation>
        <Scale>
            <x>1</x>
            <y>1</y>
            <z>1</z>
        </Scale>
        <Link>
            <href>./theGeomFile.dae</href>
        </Link>
    </Model>
</Placemark>
</kml>
  1. create a file "theGeomFile.dae" containing the following:
<?xml version="1.0" encoding="utf-8"?>
<COLLADA xmlns="http://www.collada.org/2005/11/COLLADASchema" version="1.4.1" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <asset>
    <unit name="meter" meter="1"/>
    <up_axis>Z_UP</up_axis>
  </asset>
  <library_geometries>
    <geometry id="Plane-mesh" name="Plane">
      <mesh>
        <source id="Plane-mesh-positions">
          <float_array id="Plane-mesh-positions-array" count="12">-10 -10 0 10 -10 0 -10 10 0 10 10 0</float_array>
          <technique_common>
            <accessor source="#Plane-mesh-positions-array" count="4" stride="3">
              <param name="X" type="float"/>
              <param name="Y" type="float"/>
              <param name="Z" type="float"/>
            </accessor>
          </technique_common>
        </source>
        <vertices id="Plane-mesh-vertices">
          <input semantic="POSITION" source="#Plane-mesh-positions"/>
        </vertices>
        <triangles count="2">
          <input semantic="VERTEX" source="#Plane-mesh-vertices" offset="0"/>
          <p>1 2 0 1 3 2</p>
        </triangles>
      </mesh>
    </geometry>
  </library_geometries>
  <library_visual_scenes>
    <visual_scene id="Scene" name="Scene">
      <node id="Plane" name="Plane" type="NODE">
        <matrix sid="transform">1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1</matrix>
        <instance_geometry url="#Plane-mesh" name="Plane"/>
      </node>
    </visual_scene>
  </library_visual_scenes>
  <scene>
    <instance_visual_scene url="#Scene"/>
  </scene>
</COLLADA>
  1. put them in the same directory, drag and drop the kml into Google Earth, it should work.

4) The structure of your geometry

What you already have is really the complicated part : how to generate the grid around the line. I think your idea to discretize in circles along the sight line with different radius is the way to go.

Then you need to discretize each circle with the same number of points. This will allow you to define a face of the ellipsoïde between two adjacent circles by joining 4 points : 2 on the first circle and 2 on the second. (see the picture below). A face in basic 3D (such as dae file) is a triangle so each of your square faces will be made of 2 triangles.enter image description here
enter image description here
enter image description here

5) how it relates to the DAE file

In the DAE, all you need to do is define each vertex and define the sequence of three vertices that defines a face, repeat the process for each square face of your ellipsoïd and you are done.

The only things you (may) need to modify from the template :

  1. In assets, be aware of the units, everything translates to meter, meter="1" means 1 geometry unit equals a meter. Z_up is the geometry convention.
<unit name="meter" meter="1"/>
    <up_axis>Z_UP</up_axis>
  1. All the vertices of your geometry must be written in one line in <float array> you need to put the number of vertices*nb of coordinates in count and then the coordinates of your vertices.
    They are written in X Y Z order with a space between them (e.g : X1 Y1 Z1 X2 Y2 Z2 X3 Y3 Z3). In the template it is a plane so 4 vertices with 3 coordinates = 12 in the count
<float_array id="Plane-mesh-positions-array" count="12">-10 -10 0 10 -10 0 -10 10 0 10 10 0</float_array>`
  1. You need to list the vertices triples that form a triangle face in the line <p> inside <triangles>, you will also need to update the count. In the template, the plane is made of 2 triangles : triangle 1 = vertices 1,2,0 and triangle 2 = vertices 1,3,2. The numbers refer to the order of the vertices in <float_array>. It is defined at line :
<p>1 2 0 1 3 2</p>
  1. the other lines should remain the same, especially the id and url as they are meant to work together (urls refers to ids).

In the end, this part is just a matter of string formatting

6) how to use the KML file

After you have done this (which is the longest part) You handle the geographic shift in the KML file. Inside the <LookAt>, you will only find parameters for display in GoogleEarth (that is to say: camera placement)
what is important is the location, scale and rotations in the <Model>. it is basically the transform of your cartesian origin and grid. So lat and long are the coordinate of your geometry origin.
parameters seem to be self-explanatory, for more details, see the to sources below

Sources

My solution is based on the descriptions of kml and dae provided here :
https://sites.google.com/site/timlingerscolladasketchup/home/disection-of-the-collada-format https://developers.google.com/kml/documentation/models

1
  • Thanks @ThomasJ and @Spacdman! Apologizes for the complicated question, it seemed simpler when I first posted it. Although I think the answer @ThomasJ provides is thorough, and may represent a potential solution, I lack the skills to apply it to my issue. I can easily create the quad mesh object described in R using turn3d(x = x, y = radius, n=20) but I still do not know how to export this in a meaningful format... I have figured out how to export the sfc_circumference object as a .shp file and was planning to post an update. Sorry I can't apply this answer to this OP!
    – Wu Wei
    Jan 9 at 23:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.