Discussion: Correcting LIDAR data for use in Kartapullautin

Trying to convert some CT lidar (NAD83 CT State Plane) which is in feet, to metric for processing in Kartapullautin. Using las2las.exe but it keeps resulting in tiles with far too many contours. Not sure if I'm using the wrong command line arguments, or perhaps in the wrong order to be interpreted correctly.

Anyone here process NAD83 data and care to share your las2las syntax so I can deduce where I'm going wrong? Thanks..

"Far too many contours" ... do you mean after KP processing, and are you looking at the "contours03.dxf" (raw file in ... 0.25m (maybe?) interval) or "out2.dxf", which is the final contours with (by default) 2.5m interval contours? Use "out2.dxf", if you're not.

That said: to answer your question, I have successfully used:

las2las -i INFILE.laz -o OUTFILE.laz -target_utm 18N -feet -elevation_feet

18N is something different where you live, which you probably already know. :)

Does that help? If not, get us the output of:

lasinfo FILE.laz


This should work for a quick fix: The NOAA data viewer has the most recent CT lidar data, and that site lets you specify what options you want before downloading, e.g., meters or feet, state plane or UTM, las or laz, etc. The url is:

https://www.coast.noaa.gov/dataviewer/#/imagery/se...

I'm guessing you're getting the CT lidar from the UCONN Magic site (at least it was the Magic site last time I used it). The NOAA site has the exact same lidar files, just with better user options.

Best answer, RickD! Assuming it’s public data.

It is public--your tax dollars at work....

Thanks Rick, I'll check the NOAA site. I was pulling 2016 data from UConn (now CTECO). Just trying to come up with a quick and dirty map for some explorations I need to do on Saturday.

{edit} seemed pretty painless, and easier to box a region than downloading individual tiles... We'll see what it looks like when it shows up...

moving laz file to meters:
las2txt -i tile_1300000_195000.laz -o full.laz.xyz -scale_x 0.3048 -scale_y 0.3048 -scale_z 0.3048 -parse xyzcnri
output is xyz text file, but pullauta.exe can understand it.

What is xyzcnri?
xyz -- coordinate/altitude
cnri?
c - classification -- this is one of the LIDAR classifications from the laz file, e.g. - "Lidar Return types to import"

This will work for you specifically.

las2las -i "FILEPATH\FILENAME.las" -target_utm 18N -target_precision 0.00025 -olaz -sp83 CT -survey_feet -vertical_navd88_geoid12b

I recommend if you're doing multiple tiles that you merge the tiles prior to transforming then split if desired or you may have some discontinuities along the tile edges.

At least the one file I looked at is messed up - the separate metadata is correct but the file is wrong which is why it was causing you issues.

Getting the data from NOAA is much better - the national sources do a much better job providing data.

@RickD: hah! Yeah, I meant that the data he had was public, not a private source.

@hughmac4 - I hadn't even looked at the dxf contour files. I was just trying to get a quick PNG to have in hand for a quick exploratory run later this week.

@RickD - the pre-corrected data from NOAA did the trick. It was large enough that they sent it as numerous tiles anyway, but I merged them into one large las file and processed overnight - I'll lay a few of the obvious features on top in OOM and should be good to go.

Thanks all for the input..

But of course you could have just used OCAD 12 or OCAD 2019 and the DEM Import Wizard function would have taken the .laz files as is and produced your map in minutes, not hours, I'm guessing.

I have been recently filling some holes in mapant.fi - just for fun - by processing missing tiles with my home laptop. Original 2016 coverage at left, current at right. About 15000 km2 new map now. If I got figures right Karttapullautin seems to process about 7 tiles (3x3 km tiles) per hour, that makes about 1 min/km2. I've been processing 16 tiles parallel with octacore laptop. Slow, but with thin enough data (~1 pt/m2), parallel processing and enough cores makes it not that bad.

But of course you could have just used OCAD 12 or OCAD 2019 and the DEM Import Wizard function would have taken the .laz files as is and produced your map in minutes, not hours, I'm guessing.

I couldn't figure out how to get OCAD to handle planar feet, just elevation feet, so I thought it was necessary to reproject first, unless you are never going to get open street stuff at 3.xx times scale to help with trails etc. If there's a way to skip that, please let me know!

As I understand it the original LiDAR was in the Grid coordinate system US State plane projection system (ft US). Right?
OCAD wants it in the equivalent UTM in meters so it will do that conversion for you.
There is an article by David Fisher, SLOC that explains the process. He is using OCAD 12 but the process is pretty well identical for OCAD 2019.

Gord - I am not running Windows. OCAD 2019 does run under Wine, but there are some issues so its not my first choice. In particular it doesn't play nice with HiDPI monitors so the menu text is barely legible on my system.

Too bad. I'm sorry for your difficulty which I definitely do not understand.

Chris Gkikas and I were messaging the other day about a problematic set of downloaded lidar that was causing KP to lock up in the vegetation generation step.

I noticed it had negative easting coordinates, and I guessed that might be the issue, and also it just seemed wrong.

Chris had already done what I normally recommend, which is to run lasinfo:

lasinfo -i input_file.laz -otxt

That creates a report. Sometimes and in this case, there is metadata embedded in the laz file. This was one of the report headings:

description 'WKT Projection'

We googled WKT, and there is a wiki page. It means Well-Known Text, and it's a standard way of including projection and metadata information. So I looked for projection, saw EPSG 8901, and thought I was home free. But then I looked for it again and found a different EPSG number and then several. I hope this pastes (I've added spaces and underlined the instances of EPSG):

OGC COORDINATE SYSTEM WKT:
COMPD_CS["NAD83 / Conus Albers + NAVD88 height - Geoid12B (meters)",PROJCS["NAD83 / Conus Albers",GEOGCS["NAD83",DATUM["North_American_Datum_ 1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]] ,TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM[ "Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree", 0.0174532925199433,AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4269"]],PROJECTION[ "Albers_Conic_Equal_Area"],PARAMETER ["standard_parallel_1",29.5],PARAMETER["standard_parallel_2", 45.5],PARAMETER["latitude_of_center",23],PARAMETER ["longitude_of_center",-96],PARAMETER["false_easting",0], PARAMETER["false_northing",0],UNIT["meter",1,AUTHORITY ["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY ["EPSG","5070"]],VERT_CS["NAVD88 height - Geoid12B (meters)",VERT_DATUM["North American Vertical Datum 1988",2005,AUTHORITY["EPSG","5103"]],UNIT["meter",1, AUTHORITY["EPSG","9001"]],AXIS["Up",UP],AUTHORITY ["EPSG","5703"]]]

I guessed (without evidence) that this tile was converted from a downloaded bigger file, and that something was messed up, or it was tiled from several different datasets.

If you know the "input file" projection, you can write a command line similar to this to reproject it:

las2las -i *.las -sp83 OR_N -feet -elevation_feet -target_utm 17N -target_meter -target_elevation_meter -olaz

(If you don't know, this is saying

-i *.las The input file is all the las files

-sp83 OR_N The input file's projection is Oregon North, as defined in SP83

-feet -elevation_feet The input file's coordinates are feet and the elevations are in feet (international feet, not survey feet!)

-target_utm 17N I want the output file to be in UTM Zone 17 Northern Hemisphere

-target_meter -target_elevation_meter (I want the output units to be meters---this isn't actually needed, because the default output for UTM is in meters)

-olaz Make the output file be in .laz format (saves *lots* of space on your harddrive)

Anyway, Chris tried this command line, which *doesn't specify* the "input file's projection", only the target---and it worked. So LAStools was somehow able to know what projection the data was in because of the embedded WKT metadata.

las2las -i *.laz -olaz -target_utm 12N

For those who don't know, this is saying:

-i *.laz the input file is all the "laz" files
-olaz the output should be in laz format
-target_utm 12N "the target projection to convert the input file into" is UTM zone 12 Northern Hemisphere (with meters the default units).

Two comments:

1. In the second comment in this thread, hughmac⁴ uses this "don't specify the input projection", but he does specify the units of the input are in feet. If you want to specify the units of the output file, especially if they're different from the default units, you need "target" commands:

-target_meter -target_elevation_meter

or something similar (feet or surveyfeet) (see the las2las readme file)

2. ledusledus suggested this command line to create the txt file needed by KP:

las2txt -i tile_1300000_195000.laz -o full.laz.xyz -scale_x 0.3048 -scale_y 0.3048 -scale_z 0.3048 -parse xyzcnri

I'd suggest this is bad for two reasons.

a. It's better to have the input in laz format for space reasons, and to just paste a copy of las2txt.exe into the KP directory. KP will use las2txt.exe to create whatever txt input file it wants from the laz file. (I often delete the xyz files after I run KP just to free up space.)

b. It's not a good practice to just use the 0.3048 conversion factor for feet to meters. The problem is it doesn't take into account the false easting of 500,000m for meters and 1,500,000 feet for feet. 500000m does not equal 1500000ft. If you want the output files to automatically align with aerials and other data, you need to convert to a known projection. I use UTM because it's so well known and its default units are meters. Also, the 0.3048 is only exact for international feet (1" = 2.54cm). For survey feet it's 0.3048006096. Sure, it's close enough for nearly everything, but if you can easily reproject the data and have everything line up, it's worth taking that time.

And I will admit to using 0.3048 conversions if I'm planning to do it right later. I've even just created a grid of 3m or 6m, but given the program data in feet, which means the grid is actually 3ft and 6ft. I did that back when I had no idea how to reproject, and I still do it if I want to see some rocket fast output (I use 3m/6m/10m grids even on data in both feet and meters just to quickly see outputs and diagnose problems.)