Mapping Lidar Collection Dates from LAS/LAZ Points

In an earlier post, we talked about finding the time a lidar data set or point had been collected. However, suppose you’d like to see something like a map showing collection times instead of just trying to query points. Of course, if you have the data showing the flightlines and dates, you don’t need to do this, but many end users won’t have the flightline data. Using some free open source tools and a little python scripting, we can make a table of the points with location and time and then use a GIS to map it.

Step 0: Get Your Lidar Data

This is step zero because you’ve probably already done it. If not, there are plenty of places to find lidar data. Some elevation specific repositories are given in this earlier post and you can also try data.gov. Make sure your LAS or LAZ format lidar has time recorded in adjusted GPS time. You can check that with lasinfo from LAStools. Just run lasinfo on the lidar file and look for the report on gpstime. It will look something like:

gps_time -14462805.050203 -11263004.749858

Those two numbers will be the beginning and end GPS times in seconds after subtracting a billion. If instead, your file is in GPS seconds of the week, you’ll see numbers between 0 and 604800. There is something in the file that is supposed to tell you how the time was stored, but it is often wrong. The lasinfo program will mention right after giving the GPS times if the times are outside the expected bounds. Don’t worry if you see that, it’s actually a good sign for what we need. It means you have adjusted GPS time.

Step 1: Extract

All we need for our task is the point locations and the time, so we’ll extract that from the lidar file using las2txt (again from LAStools). We really don’t need or want all the points, so we’ll just sample the data set as we extract. That can be done on the command line with

las2txt -parse xyt -i lidar_file -keep_every_nth 50 -o point_file.txt

Step 2: Convert the Times

At this point you’ve got a text file with x, y, and some weird number of seconds that somehow represents the date. It’s time to pull out some coding and do a conversion. I’m going to use python for a couple of reasons. First, if you’re running ArcGIS, you probably already have it on your machine. Second, I haven’t done much python coding, so it’s an opportunity for you to ridicule my code. The code is at the end of this post and available at ftp://ftp.coast.noaa.gov/pub/crs/kwaters/xysec2xytime.py.

You can run it on the command line with

python xysec2xytime.py point_file > xytime.txt

Or, if you want to combine steps 1 and 2 you can do

las2txt -parse xyt -i lidar_file -keep_every_nth 50 -stdout | python xysec2xytime.py > xytime.txt

Now we have a file that might have lines that look like (if x and y are geographic):
-118.8437907 34.0238348 2009-09-24 16:16:06
-118.8434954 34.0246657 2009-09-24 16:16:06
-118.8437089 34.0244767 2009-09-24 16:16:07
-118.8439095 34.0235908 2009-09-24 16:16:07

By the way, you may notice that while the python code is short, it is considerably longer than it really needs to be as it accounts for leap seconds. Those aren’t really an issue when we’re just looking for days, but I included it anyway.

Step 3: Map It

Now we’re going to bring that into ArcGIS to map it. There are other packages you could use, but since ESRI is the most common, we’ll use that. The rest of the steps may be old hat if you work in ArcGIS regularly.

Step 3a: Import the Table

Use ArcCatalog and right-click to get a context menu. Select import and then the table (either single or multiple).

This brings up the table to table tool. The file you made in step 2 is your input rows. Note that the field map did recognize the date field.

Step 3b: Display the Points

Once the data is imported and added to your ArcMap table of contents, right click for another context menu so you can display the XY points.

This will give you another menu to fill out. The X Field should already be right, but the Y Field needs to be set to Field 2. You’ll probably also want to set the coordinate system of your data so you can add a basemap later.

Hitting the OK button will display the points, but they will all be one color. All you have to do now is open the symbology, change from a single value to Categories->Unique Values, and choose Field 3 as the Value Field. Use either “Add All Values” to color all the dates or just “Add Values …” if you only want to color a few dates differently. In my example, I added all the values.

After adding a basemap to give a bit of spatial reference, my points look like the figure below.

Step 4: Done

You can see the points from this one file were flown on six different days over a three month span, even though the metadata indicated a single month for the whole project. Note that if we were trying to figure out the tide stage for the points, we would also look at the time of day, but we have that in field 4 of the table. Remember the dates and times are in UTC (Greenwich).

#!/usr/bin/env python
"""
  Author: Kirk Waters with cleanup suggestions from Kurt Schwehr.

  Converts time in adjusted GPS seconds (seconds since Jan 6, 1980 minus 1e9)
  to a date and time of day in Universal Time Coordinates.
  Inputs are expected to be a space separted stream of 
  x y seconds
  ...

  Output is a stream of 
  x y date time
  ...

  Note: This was designed to work with lidar data and the input may have
  come from a command like 'las2txt -parse xyt ...'. If the lidar file
  had time stored in seconds of the GPS week, your output indicate a week of
  September, 2011.
"""

import datetime
import fileinput

#number of seconds between the start of unix time (Jan 1, 1970) and gps time (Jan 6, 1980)
offset = 315964800


def countleaps(gpsTime):
  """Count number of leap seconds that have passed."""

  # a tuple of the gps times where leap seconds were added
  leaps = (
    46828800, 78364801, 109900802, 173059203, 252028804,
    315187205, 346723206, 393984007, 425520008, 457056009,
    504489610, 551750411, 599184012, 820108813, 914803214, 
    1025136015
    )

  nleaps = 0
  for leap in leaps:  
    if gpsTime >= leap:
       nleaps += 1
  
  return nleaps


if __name__ == '__main__':
  for line in fileinput.input():
    xvalue, yvalue, timevalue=line.split(' ')[:3]
    gpstime=float(timevalue)
    gpstime += 1e9; # unadjusted GPS time
    unixtime = gpstime + offset - countleaps(gpstime);
    datetimestr = datetime.datetime.fromtimestamp(unixtime).strftime('%Y-%m-%d %H:%M:%S')
    
    print xvalue, yvalue, datetimestr