Kirk Waters

by Monday, February 3, 2014 @ .

I’m a physical scientist at the NOAA Coastal Services Center. In my spare time, when I’m not torturing co-workers, I try to fit in some technical work on lidar processing and distribution. I also try to figure out ways to improve the Digital Coast’s data offerings in general. Somewhere in the back of my head there are still a few brain cells that remember satellite ocean color, oceanographic field work, and something about the ozone hole.

Read More >

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://csc.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).

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

#!/usr/bin/python

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=float(values[2]);
	gpstime += 1e9; # unadjusted GPS time
	nleaps=countleaps(gpstime);
	unixtime = gpstime + offset - nleaps;
	datetimestr = datetime.datetime.fromtimestamp(unixtime).strftime('%Y-%m-%d %H:%M:%S');
	outstr = values[0] + ' ' + values[1] + ' ' + datetimestr;
	print outstr;

3 Responses to “Mapping Lidar Collection Dates from LAS/LAZ Points”

  1. Kirk Waters

    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. Kirk Waters

    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.