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 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

You can run it on the command line with

python 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 > 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).

Screengrab of Arccatalog context menu.


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.

Screengrab of table to table import.


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.

Screengrab of context menu displaying 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.

Screengrab of context menu.


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.

Screengrab of symbology set.


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

Screengrab of mapped points with dates.


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).


import datetime
import fileinput

# Count number of leap seconds that have passed
def countleaps(gpsTime):
 	# 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;  # number of leap seconds prior to gpsTime
  for leap in leaps :  
    if (gpsTime >= leap):
       nleaps += 1
  return nleaps;

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

for line in fileinput.input():
	values=line.split(' ');
	gpstime += 1e9; # unadjusted GPS time
	unixtime = gpstime + offset - nleaps;
	datetimestr = datetime.datetime.fromtimestamp(unixtime).strftime('%Y-%m-%d %H:%M:%S');
	outstr = values[0] + ' ' + values[1] + ' ' + datetimestr;
	print outstr;


  1. This method works well. But there is a problem that different images turns out the same date: 2011-09-13. I used different study areas for testing ( which means different flight date for lidar), however, those different study areas got the same flight date, which is 2011-09-13. I can’t figure it out where goes wrong. Could you give me some suggestions? Thanks!


  2. If you happen to get dates in September 2011, especially if you know that’s not the right date as in Yan’s case, the problem is likely that you don’t have adjusted GPS time, you have GPS seconds of the week. If the data was recorded in that way, as much of the older data was, you need to figure out which week the data was in and use something like the LAStools free program las2las to fix it. Prior to LAS 1.2, there was no official way in LAS to record anything but seconds of the week with no way to indicate which week.


  3. Can we get the posts saved as referenced pdf for professional papers sourcing too. Every time using IE print preview crashes the system. Thanks


    • Earl,
      We can look into that. Our process to publish pdfs is definitely more involved than posting to the blog. I don’t think we’d want to go through that in general. However, if there are particular posts that people find would be useful to be redone as a pdf (and potentially expanded), we’d like to hear it.



Leave a Reply. Comments are moderated.

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s