Fixing Intensity Differences between Flightlines (“quick and dirty”)

Visiting our users on-site, such as last week at Mariano Marcos State University in Ilocos Norte in the Philippines, we sometimes come across situations as pictured below where the intensity values of the returns of one flightline are drastically different from that of other flightlines.

The intensity of returns in the left most flightline is different from that of other flightlines.

The intensity of returns in the left most flightline is different from that of other flightlines.

Using intensity rasters with such dark strips as an additional input for land cover classification may likely make things worse. Radiometrically “correct” intensity calibration is a complex topic and may not always be possible to do using only the LAZ files without meta information such as the internals of the scanning system and the aircraft trajectory. However, we now describe a “quick and dirty” fix to the situation shown above so that the intensity grids (that were computed as averages of first return intensities) at least “look” as sensible as for the one square tile (shown below) that was corrected by a simple multiplication with 5 for all intensities of the dark strip.

Simply multiplying all intensities of the dark flightline with 5 seems to "fix" the issue.

Simply multiplying all intensities of the dark flightline with 5 seems to “fix” the issue for our test tile.

The number 5 was determined by a quick glance at the intensity histograms that we can generate with lasinfo. We decide to only look at single returns as we expect them to have a higher correlation: Their locations are more likely to be “seen similarly” from and their energy is more likely “reflected similarly” to different flightlines compared to that of multiple returns.

lasinfo -i strip1.laz strip2.laz strip3.laz ^
        -keep_single ^
        -histo intensity 1 ^
        -nmm -nh -nv ^
        -odix _histo_int -otxt

The resulting histograms for the dark ‘strip1.laz’ is quite different from that of the much brighter ‘strip2.laz’ and ‘strip3.laz’. The average single return intensity for the dark ‘strip1.laz’ is a meager 5.13 whereas the  brighter ‘strip2.laz’ and ‘strip3.laz’ have similar averages of 24.15 and 24.50 respectively.

Draw your own conclusion about which scale factor to use. We have the choice to match either the peak of the histograms or their averages. Scaling the peak of 3 for ‘strip1.laz’ to match the 25 of the other two strips is too much of an upscaling. But the average 24.15 divided by 5.13 gives a potential scale of 4.71 and the average 24.50 divided by 5.13 gives a potential scale of 4.77 and we already saw a multiplication by 5 giving reasonable results. So this is how we can fix the intensity:

las2las -i strip1.laz ^
        -scale_intensity 4.75 ^
        -odix _corr_int -olaz

But what if your data is already in tiles? How can you adjust only the intensity of those returns that are from the flightline 1? Assuming that your flightline information is properly stored in the point source ID field of every point this is easily done with the new ‘-filtered_transform’ in LAStools using las2las on as many cores as you have as follows:

las2las -i tiles/*.laz ^
        -keep_point_source 1 ^
        -filtered_transform ^
        -scale_intensity 4.75 ^
        -odir tiles_corr -olaz ^
        -cores 8

This is not currently exposed in the GUI of las2las but you can simply add it by typing it into the ‘RUN’ pop-up window as shown below.

Scaling only the intensities of flightline 1 by 4.9 using the new '-filtered_transform'.

Scaling only the intensities of flightline 1 by 4.9 using the new ‘-filtered_transform’.

After this “quick and dirty” intensity correction we again ran lasgrid as follows:

lasgrid -i tiles_corr/*.laz ^
        -gray -set_min_max 0 60 ^
        -odir tiles_int_rasters -opng ^
        -cores 8

And the result is shown below. The obvious flightline-induced discontinuity in the intensities has pretty much disappeared. Do you have similar flightline-related intensity issues? We like to hear from you whether this technique works or if we need to implement something more clever in the future …

lasgrid_intensity_differences_3

Create proper LAS 1.4 files with LAStools (for free)

So your LiDAR processing workflow produces beautiful LAS or LAZ output. You think the files are nothing short of perfect. But they are in LAS 1.2 format and the tender document explicitly requests delivery in the latest LAS 1.4 format. The free and open source las2las module of LAStools – also known as the Swiss Army knife of LiDAR processing – can easily up-convert them.

There are two possible scenarios: (1) The client wants the data as LAS 1.4 files but does not care about which point type is used. (2) The client wants LAS 1.4 and specifically asks for one of the new point types such as point type 6.

Let us assume you have LAS 1.2 content with point type 1 like these LAZ files from Martin county in Minnesota. Use these four tiles 5126-05-42.laz5126-05-43.laz5126-06-42.laz, and 5126-06-43.laz to repeat the LAS 1.4 up-conversion we are doing in the following. You will need to use version 160110 of LAStools (or newer) that you can download here.

A report generated with lasinfo shows geo-referencing information stored as GeoTIFF keys in the first VLR. Then there are 10 useless numbers stored in the second VLR that seem left-over from an earlier version of geo-referencing without EPSG codes. The two vendor-specific ‘NIIRS10’ VLRs should probably be removed unless they are meaningful to you. The LAS header is followed by 1564 user-defined bytes (sometimes used as “padding”) that we can also delete before delivery.

D:\LAStools\bin>lasinfo -i martin\5126-05-42.laz
reporting all LAS header entries:
 file signature: 'LASF'
 file source ID: 0
 global_encoding: 0
 project ID GUID data 1-4: E603549E-B74A-4696-3BAD-EA49F643C221
 version major.minor: 1.2
 system identifier: 'LAStools (c) by Martin Isenburg'
 generating software: 'lassort (110815) unlicensed'
 file creation day/year: 119/2011
 header size: 227
 offset to point data: 2163
 number var. length records: 4
 point data format: 1
 point data record length: 28
 number of point records: 13772409
 number of points by return: 13440460 240598 78743 12608 0
 scale factor x y z: 0.01 0.01 0.01
 offset x y z: 0 0 0
 min x y z: 361818.33 4855898.31 349.89
 max x y z: 364400.98 4859420.79 404.22
variable length header record 1 of 4:
 reserved 43707
 user ID 'LASF_Projection'
 record ID 34735
 length after header 40
 description 'by LAStools of Martin Isenburg'
 GeoKeyDirectoryTag version 1.1.0 number of keys 4
 key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected
 key 3072 tiff_tag_location 0 count 1 value_offset 26915 - ProjectedCSTypeGeoKey: NAD83 / UTM 15N
 key 3076 tiff_tag_location 0 count 1 value_offset 9001 - ProjLinearUnitsGeoKey: Linear_Meter
 key 4099 tiff_tag_location 0 count 1 value_offset 9001 - VerticalUnitsGeoKey: Linear_Meter
variable length header record 2 of 4:
 reserved 43707
 user ID 'LASF_Projection'
 record ID 34736
 length after header 80
 description 'GeoTiff double parameters'
 GeoDoubleParamsTag (number of doubles 10)
 500000 0 -93 0.9996 0 1 6.37814e+006 298.257 0 0.0174533
variable length header record 3 of 4:
 reserved 43707
 user ID 'NIIRS10'
 record ID 1
 length after header 26
 description 'NIIRS10 Tile Index'
variable length header record 4 of 4:
 reserved 43707
 user ID 'NIIRS10'
 record ID 4
 length after header 10
 description 'NIIRS10 Timestamp'
the header is followed by 1564 user-defined bytes
LASzip compression (version 2.1r0 c2 50000): POINT10 2 GPSTIME11 2
reporting minimum and maximum for all LAS point record entries ...
 X 36181833 36440098
 Y 485589831 485942079
 Z 34989 40422
 intensity 1 4780
 return_number 1 4
 number_of_returns 1 4
 edge_of_flight_line 0 0
 scan_direction_flag 0 1
 classification 1 10
 scan_angle_rank -24 24
 user_data 0 32
 point_source_ID 5401 8015
 gps_time 243852.663596 403514.323142
number of first returns: 13440460
number of intermediate returns: 91408
number of last returns: 13440349
number of single returns: 13199808
overview over number of returns of given pulse: 13199808 323647 198449 50505 0 0 0
histogram of classification of points:
 6560507 unclassified (1)
 6744481 ground (2)
 401464 medium vegetation (4)
 55433 building (6)
 100 noise (7)
 10157 water (9)
 267 rail (10)

(1) Create LAS 1.4 files with point type 1

The las2las command shown below turns a folder of LAS 1.2 files into a folder of  LAS 1.4 files with option ‘-set_version 1.4’ without changing the point type. It removes the last three VLRs with option ‘-remove_vlrs_from_to 1 3’ and the user-defined bytes in the LAS header with option ‘-remove_padding’ (*). The result is as LAZ with option ‘-olaz’ and the task is run on up to 4 CPUs in parallel with option ‘-cores 4’.

(*) option ‘-remove_padding’ is called  ‘-remove_extra’ in older versions of LAStools.

las2las -i martin/*.laz ^
        -remove_vlrs_from_to 1 3 ^
        -remove_padding ^
        -set_version 1.4 ^
        -odir martin_LAS14_pt1 -olaz ^
        -cores 4

(2) Create LAS 1.4 files with point type 6

The las2las command shown below turns a folder of LAS 1.2 files into a folder of  LAS 1.4 files with option ‘-set_version 1.4’ and changes the point type from 1 to 6 with option ‘-set_point_type 6’. As before some VLRs and the header padding is removed. It also adds a second description of the georeferencing information by rewriting the existing information stored as GeoTIFF tags into a properly formatted OGC WKT string with option ‘-set_ogc_wkt’. Thanks to ESRI’s assault on the LAZ format the output still needs to be in LAS. It can be compressed with laszip in a subsequent step using the ‘LAS 1.4 compatibility mode‘ sponsored by NOAA.

las2las -i martin/*.laz ^
        -remove_vlrs_from_to 1 3 ^
        -remove_padding ^
        -set_version 1.4 ^
        -set_point_type 6 ^
        -set_ogc_wkt ^
        -odir martin_LAS14_pt6 -olas ^
        -cores 4

laszip -i martin_LAS14_pt6/*.las -cores 4

Below the output of lasinfo for the “upgraded” LAS 1.4 file with the OGC WKT string. Some of the key changes have been marked in red.

D:\LAStools\bin>lasinfo -i martin_LAS14_pt6\5126-05-42.las
reporting all LAS header entries:
 file signature: 'LASF'
 file source ID: 0
 global_encoding: 16
 project ID GUID data 1-4: E603549E-B74A-4696-3BAD-EA49F643C221
 version major.minor: 1.4
 system identifier: 'LAStools (c) by rapidlasso GmbH'
 generating software: 'las2las (version 160110)'
 file creation day/year: 119/2011
 header size: 375
 offset to point data: 1135
 number var. length records: 2
 point data format: 6
 point data record length: 30
 number of point records: 0
 number of points by return: 0 0 0 0 0
 scale factor x y z: 0.01 0.01 0.01
 offset x y z: 0 0 0
 min x y z: 361818.33 4855898.31 349.89
 max x y z: 364400.98 4859420.79 404.22
 start of waveform data packet record: 0
 start of first extended variable length record: 0
 number of extended_variable length records: 0
 extended number of point records: 13772409
 extended number of points by return: 13440460 240598 78743 12608 0 0 0 0 0 0 0 0 0 0 0
variable length header record 1 of 2:
 reserved 43707
 user ID 'LASF_Projection'
 record ID 34735
 length after header 40
 description 'by LAStools of Martin Isenburg'
 GeoKeyDirectoryTag version 1.1.0 number of keys 4
 key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected
 key 3072 tiff_tag_location 0 count 1 value_offset 26915 - ProjectedCSTypeGeoKey: NAD83 / UTM 15N
 key 3076 tiff_tag_location 0 count 1 value_offset 9001 - ProjLinearUnitsGeoKey: Linear_Meter
 key 4099 tiff_tag_location 0 count 1 value_offset 9001 - VerticalUnitsGeoKey: Linear_Meter
variable length header record 2 of 2:
 reserved 43707
 user ID 'LASF_Projection'
 record ID 2112
 length after header 612
 description 'by LAStools of rapidlasso GmbH'
 WKT OGC COORDINATE SYSTEM:
 PROJCS["NAD83 / UTM 15N",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-93],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","26915"]]
reporting minimum and maximum for all LAS point record entries ...
 X 36181833 36440098
 Y 485589831 485942079
 Z 34989 40422
 intensity 1 4780
 return_number 1 4
 number_of_returns 1 4
 edge_of_flight_line 0 0
 scan_direction_flag 0 1
 classification 1 10
 scan_angle_rank -24 24
 user_data 0 32
 point_source_ID 5401 8015
 gps_time 243852.663596 403514.323142
 extended_return_number 1 4
 extended_number_of_returns 1 4
 extended_classification 1 10
 extended_scan_angle -4000 4000
 extended_scanner_channel 0 0
number of first returns: 13440460
number of intermediate returns: 91408
number of last returns: 13440349
number of single returns: 13199808
overview over extended number of returns of given pulse: 13199808 323647 198449 50505 0 0 0 0 0 0 0 0 0 0 0
histogram of classification of points:
 6560507 unclassified (1)
 6744481 ground (2)
 401464 medium vegetation (4)
 55433 building (6)
 100 noise (7)
 10157 water (9)
 267 rail (10)