[R] Generating time progressing line for Google Earth

[fbielejec]

Dear,

I am trying to visualise a time-progressing line (it's supposed to
represent spread patterns) using brew package and Google Earth. 

The idea is to have a function which takes start and end point
geographic coordinates, as well as number of intervals to chop the path
up, and returns the collection of points segmenting this line.

Unfortunately my calculations fail for large distances, as the
generated lines are nowhere near being straigth (run the code to see the
example). 

My R code so far:  

############
#---LIBS---#
############
library(brew)

###############################
#---GREAT CIRCULAR DISTANCE---#
###############################
#degrees to radians
Radians <- function (degree) {
    radian = degree * (pi/180.0)
    return(radian)
}

#radians to degrees
Degrees <- function (radian) {
    degree = radian * (180.0/pi)
    return(degree)
}

# Calculates the distance between two points using the
# Spherical Law of Cosines
gcd.slc <- function(lon1, lat1, lon2, lat2) {

	  R = 6371.0 # Earth mean radius [km]

          lon1 = Radians(lon1)  
          lat1 = Radians(lat1)
          lon2 = Radians(lon2)
          lat2 = Radians(lat2)

	  d = acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) *
	  cos(lon2-lon1)) * R 
	  return(d) # Distance in km
	}

########################
#---PROGRESSING LINE---#
########################
GenerateLineSegments <- function(start_point.Long, start_point.Lat,
end_point.Long, end_point.Lat, numberOfIntervals) {

coords <- matrix(NA, numberOfIntervals, 2)

full_distance = gcd.slc(lat1 = start_point.Lat, lon1 =
start_point.Long, lat2 = end_point.Lat, lon2 = end_point.Long)
distance_slice = full_distance / numberOfIntervals

for(i in 1 : numberOfIntervals) {

distance = i * distance_slice

ang_dist = distance / 6371.0

lon_1 = Radians(start_point.Long);
lat_1 = Radians(start_point.Lat);
lat_2 = Radians(end_point.Lat);
lon_diff = Radians(end_point.Long - start_point.Long);

# First calculate the bearing 
bearing = atan2( sin(lon_diff) * cos(lat_2), (cos(lat_1) * sin(lat_2))
- (sin(lat_1) * cos(lat_2) * cos(lon_diff)) );

# Then use the bearing and the start point to find the destination
new_lat_rad = asin(sin(lat_1) * cos(ang_dist) + cos(lat_1) *
sin(ang_dist) * cos(bearing)); new_lon_rad = lon_1 +
atan2( sin(bearing) * sin(ang_dist) * cos(lat_1), cos(ang_dist) -
sin(lat_1) * sin(lat_2) );

# Convert from radians to degrees
new_lat = Degrees(new_lat_rad);
new_lon = Degrees(new_lon_rad);

coords[i, 2] = new_lat
coords[i, 1] = new_lon
  }
  
  return(coords)
}


###################
#---BREWING KML---#
###################
#A	39.5	-4.5
#E	44.75	-107.5

startLong = -4.5
startLat = 39.5

endLong = -107.5
endLat = 44.75

numberOfIntervals = 100

coords <- GenerateLineSegments(startLong, startLat, endLong, endLat,
numberOfIntervals) coords <- as.data.frame(coords)
names(coords) <- c("Long", "Lat")

seg <- data.frame(matrix(NA, nrow(coords) - 1, 5))
names(seg) <- c("x", "y", "xend", "yend", "segment")

for (i in 1 : (nrow(coords) - 1 ) ) {
seg[i, ]$x = coords[i, ]$Long 
seg[i, ]$y = coords[i, ]$Lat
seg[i, ]$xend = coords[i + 1, ]$Long
seg[i, ]$yend = coords[i + 1, ]$Lat 
seg[i, ]$segment = paste(i)
}

#seg

brew(file = "LinesTemplate.xml", output = "RKMLoutput.kml" )


...and the kml template file to be used with brewer:


<?xml version="1.0" encoding="utf-8" ?>
<kml xmlns="http://www.opengis.net/kml/2.2">
<Document>
<Folder><name>lines</name>
<% for(i in 1:nrow(seg)){ %>
  <Placemark>
  <TimeSpan>
	<begin><%= paste(1980+i, "01", "01", sep="-")%></begin>
    </TimeSpan>
          <name><%=paste(seg$x[i], ",", seg$y[i]," : ", seg$xend[i],
",", seg$yend[i], sep="")%></name>
<Style><LineStyle><color>7f9bffff</color><width>3.5</width></LineStyle></Style>
<LineString> <tessellate>1</tessellate>
                  <coordinates><%=seg[i,]$x%>, <%=seg[i,]$y%>, 0.0 
                               <%=seg[i,]$xend%>, <%=seg[i,]$yend%>, 0.0
      </coordinates></LineString>
  </Placemark>
<% } %> 
     </Folder> 
</Document></kml>


What is going wrong here?


-- 
while(!succeed) { try(); }