Plotting temporal trajectories on Google Earth with Python

In this post I’m going to show you how to display temporal trajectories on Google Earth using the python package simplekml. This package enables us to generate KML videos that we can then plot and record on Google Earth. The python script kmlMovie.py described in this post is available on my GitHub account.

This script can be used to make videos like the one below. The video quality is not very high but it gives you an idea of the kind of vizs that you can do with simplekml.

https://www.youtube.com/watch?v=2msQManvbBE

First, we need a temporal network divided into two csv files, Nodes.csv and Links.csv. Nodes.csv is composed of 3 columns providing geographical information about the nodes:

  1. Node ID
  2. Longitude of the node
  3. Latitude of the node

You can use the piece of code below to import the file into the dictionary nodes.

nodes_file = open('Nodes.csv')                               
col = nodes_file.readline().rstrip('\n\r').split(';')
nodes = {}
for line in nodes_file:
        attr = line.rstrip('\n\r').split(';')
        nodes[str(attr[0])] = [float(attr[1]),float(attr[2])]

Links.csv is composed of 5 columns providing information about the “temporal” link between nodes defined in Nodes.csv. The links should be sorted by starting time.

  1. Link ID (numeric)
  2. Node of origin
  3. Node of destination
  4. Departure time of the link from the node of origin (date format ‘2014-01-01T08:00:00Z’)
  5. Arrival time of the link to the node of destination (date format ‘2014-01-01T08:00:00Z’)

You can use the piece of code below to import the file into the dictionary links. If you want to reproduce the example of movie given above, the temporal network is available here.

links_file = open('Links.csv')                               
col = links_file.readline().rstrip('\n\r').split(';')
links = {}
for line in links_file:
        attr = line.rstrip('\n\r').split(';')
        links[attr[0]] = [str(attr[1]),str(attr[2]),str(attr[3]),str(attr[4])]

The script has also 7 parameters used to define the time window (start and end of the movie), the nodes’ and links’ remaining time on screen, the nodes’ and links’ styles and the name of the output file. Their values can be set directly into the Python script kmlMovie.py.

  1. begin: Starting time of the movie (date format ‘2014-01-01T08:00:00Z’)
  2. end: Ending time of the movie (date format ‘2014-01-01T08:00:00Z’)
  3. delay_node: Number of seconds the nodes remain visible after their creation
  4. delay_link: Number of seconds the links remain visible after their creation
  5. hrefnodestyle: Hyperlink leading to the node style image file
  6. rgbcolorlink: Define the links color with an rgb code
  7. output_file: Name of the output file
begin = '2016-01-01T07:59:00Z'  
end = '2016-01-01T08:09:00Z'
delay_node = 50
delay_link = 50
hrefnodestyle = 'http://bit.ly/21awbeB'
rgbcolorlink = [70,147,226]
output_file = 'Movie.kml'

The parameters begin and end are in date format. They form a time interval containing all the departure and arrival times. The parameters delay_node and delay_link can be tuned to adjust the time during which the nodes and links remain visible before disappearing from the screen (50 seconds in the example). Finally, the parameter hrefnodestyle is a hyperlink leading to an image file (Cartman in this example). The color of the links is defined with the parameter rgbcolorline. The size of the nodes and the width of the link can also be tuned (see below). With little effort we could easily modify the code to set the nodes’ and links’ parameters individually in order to obtain different style and remaining times on screen.

Then, we need to define three auxiliary functions to lighten the code. TZtoUNIX and UNIXtoTZ are used to convert date to unix time stamp format and inversely. You will need to import the package datetime to use this two functions.

def TZtoUNIX(date):
  sec = datetime.strptime(date,'%Y-%m-%dT%H:%M:%SZ')
  sec=(sec-datetime(1970,1,1)).total_seconds()      
  return sec
def UNIXtoTZ(sec):
  date = datetime.utcfromtimestamp(sec).strftime('%Y-%m-%dT%H:%M:%SZ') 
  return date

We then define the IntermediatePoint function in order to determine the coordinates of an intermediate point on a great circle path between two points. The location of the intermediate point on the arc can be tuned with the parameter fraction to obtain a point located at a distance fraction*(great circle distance between the nodes of origin and destination) from the origin. We need this function to build progressively the temporal link between the node of origin and the node of destination by increasing the value of the parameter fraction from 0 to 1 by step of 0.05 (see below).

def IntermediatePoint(origin, dest, fraction):
  
  #Earth spherical radius in meters
  r = 6371*1000
  
  #Calculate great circle distance between origin and dest in radians
  lon1 = (math.pi/180.)*origin[0]
  lat1 = (math.pi/180.)*origin[1] 
  lon2 = (math.pi/180.)*dest[0] 
  lat2 = (math.pi/180.)*dest[1]

  dlon = lon2-lon1
  dlat = lat2-lat1

  d = math.pow(math.sin(dlat/2.),2)+math.cos(lat1)*math.cos(lat2)*math.pow(math.sin(dlon/2.),2)    
  d = 2*math.asin(min(1,math.sqrt(d)))
  
  #Calculate longitude and latitude of the intermediate point
  A = math.sin((1-fraction)*d)/math.sin(d)
  B  =math.sin(fraction*d)/math.sin(d)
  
  x = A*math.cos(lat1)*math.cos(lon1) + B*math.cos(lat2)*math.cos(lon2)
  y = A*math.cos(lat1)*math.sin(lon1) + B*math.cos(lat2)*math.sin(lon2)
  z = A*math.sin(lat1)                + B*math.sin(lat2)
  
  lon = (180./math.pi)*math.atan2(y,x)
  lat = (180./math.pi)*math.atan2(z,math.sqrt(math.pow(x,2)+math.pow(y,2)))
  
  #Calculate altitude in meters in a sinusoidal trajectory
  h = d*r*0.3
  h = h*math.sin(math.pi*fraction)
  
  return [lon, lat, h]

We can now define the main function by creating a kml file called, with originality, kml in our case. The simplekml package allow us to add different types of object to this file (points and lines in our case) for which we can define a style and a temporal dimension. Any object in KML can have time data associated with it. This time data has the effect of restricting the visibility of the data set to a given time period or point in time. Although the complete data set is fetched when the KML file is loaded, the time slider in the Google Earth user interface controls which parts of the data are visible (see this webpage for more detail).

kml = simplekml.Kml()

First, we define the nodestyle with the parameter hrefnodestyle, an hyperlink leading to the node style image file (Cartman in our case). We can also adjust the size of the nodes (2 in our example).

nodestyle = simplekml.Style()
nodestyle.iconstyle.icon.href = hrefnodestyle
noestyle.iconstyle.scale = 2

Then we define the time window (starting and ending time of the movie). By default, the time windows of the KML movie depends on the timespan of the KML objects. If we want to set it manually we need to create an hidden node with the dummynodestyle define below.

dummynodestyle = simplekml.Style()
dummynodestyle.iconstyle.icon.href = ''

This hidden node will be “visible” from the beginning of the movie (set with the parameter begin) to the end of the movie (set with the parameter end). Since this node will not appear on screen we can choose its lon/lat coordinates randomly among the nodes (first line of the piece of code below).

dummynode = [nodes[nodes.keys()[0]][0],nodes[nodes.keys()[0]][1]]
dummynode = kml.newpoint(name='', coords=[tuple(dummynode)])
dummynode.timespan.begin = begin
dummynode.timespan.end = end
dummynode.style = dummynodestyle

We can now loop over the links sorted by IDs (i.e. starting time).

for key in sorted(links):

For each link, we extract the departure time, arrival time and the duration of the trip in second using the function TZtoUNIX.

    link = links[key] 
    time_departure = link[2]    
    time_arrival = link[3]
    duration = TZtoUNIX(link[3]) - TZtoUNIX(link[2])

We then create the nodes of origin and destination of this link which will be visible from the departure and arrival time of the link, respectively. Both nodes will remain visible during delay_node seconds.

    node1 =[nodes[link[0]][0],nodes[link[0]][1]]
    node2 =[nodes[link[1]][0],nodes[link[1]][1]]
      
    origin = kml.newpoint(name='', coords=[tuple(node1)])
    origin.timespan.begin = time_departure
    origin.timespan.end = UNIXtoTZ(TZtoUNIX(time_departure) + delay_node)
    origin.style = nodestyle
    
    dest = kml.newpoint(name='', coords=[tuple(node2)])
    dest.timespan.begin = time_arrival
    dest.timespan.end = UNIXtoTZ(TZtoUNIX(time_arrival) + delay_node)
    dest.style = nodestyle  

We finally build the temporal link, progressively, piece by piece. The link is divided into 20 pieces. The first and last points of the link are the node of origin and destination, respectively. The 19 intermediate points are computed with the IntermediatePoint function defined above. We use the parameter rgbcolorlink to define the color of the link. The width of the link can also be adjusted (5 in this example). As it is the case for the nodes, the pieces of link will remain visible on screen during delay_link seconds.

    p0 = (node1[0], node1[1], 0)    #First point 
    for f in np.arange(0,1.05,0.05):
        p = IntermediatePoint(node1, node2, f)
        pi = (p[0],p[1],p[2])      
        line = kml.newlinestring(coords=[p0,pi], tessellate=1, altitudemode='relativeToGround')
        line.style.linestyle.width = 5
        line.style.linestyle.color = simplekml.Color.rgb(rgbcolorlink[0],rgbcolorlink[1],rgbcolorlink[2], a=255) 
        line.timespan.begin = UNIXtoTZ(TZtoUNIX(time_departure) + int(f*duration))
        line.timespan.end = UNIXtoTZ(TZtoUNIX(time_departure) + int(f*duration) + delay_link)     
        p0=pi 

The KML movie is finally saved into the output_file.

kml.save(output_file)
comments powered by Disqus