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.
Best Answer
The question concerns what kinds of curves deserve an implicitly exact representation rather than a discretized approximation. The crux of the matter is this: to be successful, the class of curves you support in this manner must be closed under the class of curve- and polygon-creation operations supported in the GIS.
These operations include:
Buffering. In this process, you need to construct curves that are parallel to features. ("Parallel" means in the sense of maintaining a fixed distance.) This includes circles and portions thereof (for buffering points), oblique parallels (which are curves equidistant to geodesics on the spheroid, and can reduce to isolated points in special cases), and concentric circles. On the sphere (but not, generally, on an ellipsoid) the oblique parallels are themselves circles.
Polygons of influence (Thiessen polygons; Voronoi polygons; Dirichlet cells). To construct the Thiessen polygons for a collection of point features we need to find bisecting lines, which are geodesics (they are straight); but for a collection of other kinds of features, such as points and segments, the boundaries of the Thiessen polygons include portions of parabolas (in the plane). Maybe you don't want to support this...
Set-theoretic overlays (intersection, union, difference, complement). These operations do not create any new kinds of curves.
Parallel translation and rotation. These are usually not possible to perform exactly on an ellipsoid (because it is not a homogeneous space), but are straightforward on the sphere. On the sphere, these operations do not create new kinds of curves.
The really problematic class of curves you propose consists of the general rhumb lines (loxodromes). Lines of latitude are rhumb lines but (on the sphere at least) they are also circles, so they present no additional problem. But general rhumb lines are complicated beasts: if they are not meridians or parallels, they spiral into one pole or the other. Buffers and parallel translations of rhumb lines will be genuinely new types of curves. You would have to represent these results as broken segments of lines and circles, which would defeat your purpose (and be fairly difficult to compute). Therefore I suggest not trying to support rhumb lines exactly.
In sum, it looks like you can be successful in your program if (a) you work on a spherical model of the earth rather than the more general ellipsoidal ("spheroidal") model and (b) you limit certain constructions such as Thiessen polygons (and medial axes, which are closely related) to collections of points.