Another German State Goes Open LiDAR: Saxony

Finally some really good news out of Saxony. ūüė䬆After North Rhine-Westphalia and Thuringia released the first significant amounts of open geospatial data in Germany in a one-two punch in January 2017, we now have a third German state opening their entire tax-payer-funded geospatial data holdings to the tax-paying public via a simple and very easy-to-use online download portal. Welcome to the open data party, Saxony!!!

Currently available via the online portal are the LiDAR-derived raster Digital Terrain Model (DTM) at 1 meter resolution (DGM 1m) for everything flown since 2015 and and at 2 meter resolution (DGM 2m) or 20 meter resolution (DGM 20m) for the entire state. The horizontal coordinates use UTM zone 33 with ETRS89 (aka EPSG code 25833) and the vertical coordinate uses the “Deutsche Haupth√∂hennetz 2016” or “DHHN2016” (aka EPSG code 7837). Also available are orthophotos at 20 cm (!!!) resolution (DOP 20cm).

dgm_1000_rdax_87

Overview of current LiDAR holdings. Areas flown 2015 or later have LAS files and 1 meter rasters. Others have LiDAR as ASCII files and lower resolution rasters.

Offline – by ordering through either this online form or that online form – you can also get the 5 meter DTM and the 10 meter DTM, the raw LiDAR point clouds, LiDAR intensity rasters, hill-shaded DTM rasters, as well as the 1 meter and the 2 meter Digital Surface Model (DSM) for a small administrative fee that ranges between 25 EUR and 500 EUR depending on the effort involved.

Our immediate thought is to get a copy on the entire raw LiDAR points clouds (available as LAS 1.2 files for all  data acquired since 2015 and as ASCII text for earlier acquisitions) and find some portal willing to hosts this data online. We are already in contact with the land survey of Saxony to discuss this option and/or alternate plans.

Let’s have a look at the data. First we download four 2 km by 2 km tiles of the 1 meter DTM raster for an area surrounding the so called “Greifensteine” using the interactive map of the download portal, which are provided as simple XYZ text. Here a look at the contents of one ot these tiles:

more Greifensteine\333525612_dgm1.xyz
352000 5613999 636.26
352001 5613999 636.27
352002 5613999 636.28
352003 5613999 636.27
352004 5613999 636.24
[...]

Note that the elevation are not sampled in the center of every 1 meter by 1 meter cell but exactly on the full meter coordinate pair, which seems especially common  in German-speaking countries. Using txt2las we convert these XYZ rasters to LAZ format and add geo-referencing information for more efficient subsequent processing.

txt2las -i greifensteine\333*_dgm1.xyz ^
        -set_scale 1 1 0.01 ^
        -epsg 25833 ^
        -olaz

Below you see that going from XYZ to LAZ reduces the amount of  data from 366 MB to 10.4 MB, meaning that the data on disk becomes over 35 times smaller. The ability of LASzip to compress elevation rasters was first noted during the search for missing airliner MH370 and resulted in our new LAZ-based compressor for height grid called DEMzip.  The resulting LAZ files now also include geo-referencing information.

96,000,000 333525610_dgm1.xyz
96,000,000 333525612_dgm1.xyz
96,000,000 333545610_dgm1.xyz
96,000,000 333545612_dgm1.xyz
384,000,000 bytes

2,684,820 333525610_dgm1.laz
2,590,516 333525612_dgm1.laz
2,853,851 333545610_dgm1.laz
2,795,430 333545612_dgm1.laz
10,924,617 bytes

Using blast2dem we then create a hill-shaded version of the 1 meter DTM in order to overlay a visual representation of the DTM onto Google Earth.

blast2dem -i greifensteine\333*_dgm1.laz ^
          -merged ^
          -step 1 ^
          -hillshade ^
          -o greifensteine.png

Below the result that nicely shows how the penetrating laser of the LiDAR allows us to strip away the forest to see interesting geological features in the bare-earth terrain.

In a second exercise we use the available RGB orthophoto images to color one of the DTM tiles and explore it using lasview. For this we download the image for the top left of the four tiles that covers the area containing the “Greifensteine” from the interactive download portal¬†for orthophotos. As the resolution of the TIF image is 20 cm and that of the DTM is only 1 meter, we first down-sample the TIF using gdalwarp of¬†GDAL.

gdalwarp -tr 1 1 ^
         -r cubic ^
         greifensteine\dop20c_33352_5612.tif ^
         greifensteine\dop1m_33352_5612.tif

If you are not yet using GDAL today is a good day to start. It nicely complements the point cloud processing functionality of LAStools for raster inputs. Next we use lascolor to give each elevation pixel of the DTM stored in LAZ format its corresponding color from the orthophoto.

lascolor -i greifensteine\333525612_dgm1.laz ^
         -image greifensteine\dop1m_33352_5612.tif ^
         -odix _rgb -olaz

Now we can view the colored DTM in LAZ format interactively with lasview or any other LiDAR viewing software and turn on the RGB colors from the orthophoto as needed to understand the scene.

lasview -i greifensteine\333525612_dgm1_rgb.laz

We thank the¬†“Staatsbetrieb Geobasisinformation und Vermessung Sachsen (GeoSN)” for giving us easy access to the 1 meter DTM and the 20 cm orthophoto¬†that we have used in this article through their new open geodata portal as open data under the user-friendly license “Datenlizenz Deutschland – Namensnennung – Version 2.0.

National Open LiDAR Strategy of Latvia humiliates Germany, Austria, and other European “Closed Data” States

Latvia, officially the Republic of Latvia, is a country in the Baltic region of Northern Europe has around 2 million inhabitants, a territory of 65 thousand square kilometers and Рsince recently Рalso a fabulous open LiDAR policy. Here is a list of 65939 tiles in LAS format available for free download that cover the entire country with airborne LiDAR with a density from 4 to 6 pulses per square meters. The data is classified into ground, building, vegetation, water, low noise, and a few other classifications. It is licensed Creative Commons CC0 1.0 Рmeaning that you can copy, modify, and distribute the data, even for commercial purposes, all without asking permission. And there is a simple and  functional interactive download portal where you can easily download individual tiles.

latvia_open_data_portal_01

Interactive open LiDAR download portal of Latvia.

We downloaded the 5 by 5 block of¬†square kilometer¬†tiles matching “4311-32-XX.las” for checking the quality and creating a 1m DTM and a 1m DSM raster. You can follow along after downloading the latest version of LAStools.

Quality Checking

We first run lasvalidate and lasinfo on the downloaded LAS files and then immediately compress them with laszip because multi-core processing of uncompressed LAS files will quickly overwhelm our file system, make processing I/O bound, and result in overall longer processing times with CPUs waiting idly for data to be loaded from the drives.

lasinfo -i 00_tiles_raw\*.las ^
        -compute_density ^
        -histo z 5 ^
        -histo intensity 256 ^
        -histo user_data 1 ^
        -histo scan_angle 1 ^
        -histo point_source 1 ^
        -histo gps_time 10 ^
        -odir 01_quality -odix _info -otxt ^
        -cores 3
lasvalidate -i 00_tiles_raw\*.las ^
            -no_CRS_fail ^
            -o 01_quality\report.xml

Despite already excluding a missing Coordinate Reference System (CRS) from being a reason to fail (the lasinfo reports show that the downloaded LAS files do not have any geo-referencing information) lasvalidate still reports a few failing files, but scrutinizing the resulting XML file ‘report.xml’ shows only minor issues.

Usually during¬†laszip compression we do not alter the contents of a file, but here we also add the EPSG code 3059 for CRS “LKS92 / Latvia TM” as we turn bulky LAS files into slim LAZ files so we don’t have to specify it in all future processing steps.

laszip -i 00_tiles_raw\*.las ^
       -epsg 3059 ^
       -cores 2

Compression reduces the total size of the 25 tiles from over 4.1 GB to below 0.6 GB.

Next we use lasgrid to visualize the last return density which corresponds to the pulse density of the LiDAR survey. We map each 2 by 2 meter pixel where the last return density is 2 or less to blue and each 2 by 2 meter pixel it is 8 or more to red.

lasgrid -i 00_tiles_raw\*.laz ^
        -keep_last ^
        -step 2 ^
        -density_16bit ^
        -false -set_min_max 2 8 ^
        -odir 01_quality -odix _d_2_8 -opng ^
        -cores 3

This we follow by the mandatory lasoverlap check for flight line overlap and alignment where we map the number of overlapping swaths as well as the worst vertical difference between overlapping swaths to a color that allows for quick visual quality checking.

lasoverlap -i 00_tiles_raw\*.laz ^
           -step 2 ^
           -min_diff 0.1 -max_diff 0.2 ^
           -odir 01_quality -opng ^
           -cores 3

The results of the quality checks with lasgrid and lasoverlap are shown below.

Raster Derivative Generation

Now we use first las2dem to create a Digital Terrain Model (DTM) and a Digital Surface Model (DSM) in RasterLAZ format and then use blast2dem to create merged and hill-shaded versions of both. Because we will use on-the-fly buffering to avoid edge effects along tile boundaries we first spatially index the data using lasindex for more efficient access to the points from neighboring tiles.

lasindex -i 00_tiles_raw\*.laz ^
         -cores 3

las2dem -i 00_tiles_raw\*.laz ^
        -keep_class 2 9 ^
        -buffered 25 ^
        -step 1 ^
        -use_orig_bb ^
        -odir Latvia\02_dtm_1m -olaz ^
        -cores 3

blast2dem -i 02_dtm_1m\*.laz ^
          -merged ^
          -hillshade ^
          -step 1 ^
          -o dtm_1m.png

las2dem -i 00_tiles_raw\*.laz ^
        -drop_class 1 7 ^
        -buffered 10 ^
        -spike_free 1.5 ^
        -step 1 ^
        -use_orig_bb ^
        -odir 03_dsm_1m -olaz ^
        -cores 3

blast2dem -i 03_dsm_1m\*.laz ^
          -merged ^
          -hillshade ^
          -step 1 ^
          -o dsm_1m.png

Because the overlaid imagery does not look as nice in our new Google Earth installation, below are the DTM and DSM at versions down-sampled to 25% of their original size.

Many thanks to SunGIS from Latvia who tweeted us about the Open LiDAR after we chatted about it during the Foss4G 2019 gala dinner. Kudos to the Latvian Geospatial Information Agency (LGIA) for implementing a modern national geospatial policy that created opportunity for maximal return of investment by opening the expensive tax-payer funded LiDAR data for re-purposing and innovation without barriers. Kudos!

Using Open LiDAR to Remove Low Noise from Photogrammetric UAV Point Clouds

We collected drone imagery of the¬†restored “Kance” tavern¬†during the lunch stop of the UAV Tartu summer school field trip (actually the organizer Marko Kohv¬†did that, as I was busy introducing students to SUP boarding). With Agisoft PhotoScan we then processed the images into point clouds below the deck of the historical “J√Ķmmu”¬†barge on the way home (actually Marko did that, because I was busy enjoying the view of the wetlands in the afternoon sun). The resulting data set with 7,855,699 points is shown below and can be downloaded here.

7,855,699 points produced with Agisoft Photoscan

Generating points using photogrammetric techniques in scenes containing water bodies tends to be problematic as dense blotches of noise points above and below the water surface are common as you can see in the picture below. Especially the low points are troublesome as they adversely affect ground classification which results in poor Digital Elevation Models (DTMs).

Clusters of low noise points nearly 2 meters below the actual surface in water areas.

In a previous article we have described a LAStools workflow that can remove excessive low noise. In this article here we use external information about the topography of the area to clean our photogrammetry points. How convenient that the Estonian Land Board has just released their entire LiDAR archives as open data.

Following these instructions here you can download the available open LiDAR for this area, which has the map sheet index 475681. Alternatively you can download the four currently available data sets here flown in spring 2010, in summer 2013, in spring 2014, and in summer 2017. In the following we will use the one flown in spring 2014.

We can view both data sets simultaneously in lasview. By adding ‘-faf’ to the command-line we can switch back and forth between the two data sets by pressing ‘0’ and ‘1’.

lasview -i Kantsi.laz ^
        -i 475681_2014_tava.laz ^
        -points 10000000 ^
        -faf

We find cut the 1 km by 1 km LiDAR tile down to a 250 m by 250 m tile that nicely surrounds our photogrammetric point set using the following las2las command-line:

las2las -i 475681_2014_tava.laz ^
        -inside_tile 681485 6475375 250 ^
        -o LiDAR_Kantsi.laz

lasview -i Kantsi.laz ^
        -i LiDAR_Kantsi.laz ^
        -points 10000000 ^
        -faf

Scrutinizing the two data sets we quickly find that there is a miss-alignment between the dense imagery-derived and the comparatively sparse LiDAR point clouds. With lasview we investigate the difference between the two point clouds by hovering over a point from one point cloud and pressing <i> and then hovering over a somewhat corresponding point from the other point cloud and pressing <SHIFT>+<i>. We measure displacements of around 2 meters vertically and of around 3 to 3.5 meter in total.

Before we can use the LiDAR points to remove the low noise from the photogrammetric points we must align them properly. For simple translation errors this can be done with a new feature that was recently added to lasview. Make sure to download the latest version (190404 or newer) of LAStools to follow the steps shown in the image sequence below.

las2las -i Kantsi.laz ^
        -translate_xyz 0.89 -1.90 2.51 ^
        -o Kantsi_shifted.laz

lasview -i Kantsi_shifted.laz ^
        -i LiDAR_Kantsi.laz ^
        -points 10000000 ^
        -faf

The result looks good in the sense that both sides of the photogrammetric roof are reasonably well aligned with the LiDAR. But there is still a shift along the roof so we repeat the same thing once more as shown in the next image sequence:

We use a suitable displacement vector and apply it to the photogrammetry points, shifting them again:

las2las -i Kantsi_shifted.laz ^
        -translate_xyz -1.98 -0.95 0.01 ^
        -o Kantsi_shifted_again.laz

lasview -i Kantsi_shifted_again.laz ^
        -i LiDAR_Kantsi.laz ^
        -points 10000000 ^
        -faf

The result is still not perfect as there is also some rotational error and you may find another software such as Cloud Compare more suited to align the two point clouds, but for this exercise the alignment shall suffice. Below you see the match between the photogrammetry points and the LiDAR TIN before and after shifting the photogrammetry points with the two (interactively determined) displacement vectors.

The final steps of this exercise use las2dem and the already ground-classified LiDAR compute a 1 meter DTM, which we then use as input to lasheight. We classify the photogrammetry points using their height above this set of ground points with 1 meter spacing: points that are 40 centimeter or more below the LiDAR DTM are classified as noise (7), points that are between 40 below to 1 meter above the LiDAR DTM are classified to a temporary class (here we choose 8) that has those points that could potentially be ground points. This will help, for example, with subsequent ground classification as large parts of the photogrammetry points – namely those on top of buildings and in higher vegetation – can be ignored from the start by a ground classification algorithm such as lasground.

las2dem -i LiDAR_Kantsi.laz ^
        -keep_class 2 ^
        -kill 1000 ^
        -o LiDAR_Kantsi_dtm_1m.bil
lasheight -i Kantsi_shifted_again.laz ^
          -ground_points LiDAR_Kantsi_dtm_1m.bil ^
          -classify_below -0.4 7 ^
          -classify_between -0.4 1.0 8 ^
          -o Kantsi_cleaned.laz

Below the results we have achieved after “roughly” aligning the two point clouds with some new lasview tricks and then using the LiDAR elevations to classify the photogrammetry points into “low noise”, “potential ground”, and “all else”.

We thank the Estonian Land Board for providing¬†open data with a permissive¬†license. Special thanks also go to the organizers of the UAV Summer School in Tartu, Estonia and the European Regional Development Fund for funding this event. Especially fun was the fabulous excursion to the Emaj√Ķe-Suursoo Nature Reserve and through to Lake Peipus aboard, overboard and aboveboard the historical barge “J√Ķmmu”. If you look carefully you can also find the barge in the photogrammetry point cloud. The photogrammetry data used here was acquired during our lunch stop.

Fun aboard and overboard the historical barge “J√Ķmmu”.

In Sweden, all they Wanted for Christmas was National LiDAR as Open Data

Let’s heat up some sweet, warm and spicy Gl√∂gg in celebration! They must have been good boys and girls up there in Sweden. Because “Jultomten” or simply ‚ÄĚTomten‚ÄĚ – how Sweden’s Santa Clause is called – is assuring a “God Jul” for all the Swedish LiDAR lovers this Christmas season.

Only a few weeks ago this tweet of ours had (mistakenly) included Sweden in a list of European countries that had released their national LiDAR archives as open data for public reuse over the past six years.

Turns out we were correct after all. Sweden has just opened their LiDAR data for free and unencumbered download. To get the data simply create a user account and browse to the ftp site for download as shown in the image sequence below.

The released LiDAR data was collected with a density of 1 to 2 pulses per square meter and is distributed in LASzip compressed LAZ tiles of 2500 by 2500 meters. The returns are classified into four classes: ground (2), water (9), low noise (7) and high noise (18). All items that can not be classified as any of the first four classes coded as left unclassified (1). The LAZ files do not contain CRS information, but this can easily be added with horizontal coordinates in SWERED99 TM (EPSG code 3006) and elevations in RH2000 height (EPSG code 5613).

Below a look with lasview at a 5 km by 5 km area that composed of the four tiles ‘18P001_67100_5800_25.laz‘, ‘18P001_67100_5825_25.laz‘, ‘18P001_67125_5800_25.laz‘ and ‘18P001_67125_5825_25.laz‘ with several of the different color modes available.

 

Some more details: The data was acquired at flying altitude of around 3000 meter with a¬†maximum scan angle of ¬Ī 20¬ļ and a minimum side overlap of 10% between the flightlines. The laser footprint on ground is below 75 centimeters with slight variation based on the flying altitude. The laser scanning survey was performed with LiDAR instruments that can provide at least three returns from the same pulse. All LiDAR returns are preserved throughout the entire production chain.

The LiDAR data comes with the incredibly Creative Commons РCC0 license, which means that you can use, disseminate, modify and build on the data Рeven for commercial purposes Рwithout any restrictions. You are free to acknowledge the source when you distribute the data further, but it is not required.

The LiDAR data will eventually cover approximately 75% of Sweden and new point clouds will continuously be added as additional scanning is performed according to the schedule shown below. The survey will be returning to scan every spot again after about 7 years.

2018-2022 LiDAR acquisition plan for Sweden

Below a lasinfo report for tile ‘18P001_67125_5825_25.laz‘. One noticeable oddity is the distribution of intensities. The histogram across all intensities with bins of size 256 shows two clearly distinct sets of intensities each with their own peak and a void of values between 3000 and 10000.

lasinfo -i 18P001_67125_5825_25.laz -cd -histo intensity 256
reporting all LAS header entries:
  file signature:             'LASF'
  file source ID:             0
  global_encoding:            1
  project ID GUID data 1-4:   00000000-0000-0000-0000-000000000000
  version major.minor:        1.2
  system identifier:          ''
  generating software:        'TerraScan'
  file creation day/year:     303/2018
  header size:                227
  offset to point data:       227
  number var. length records: 0
  point data format:          1
  point data record length:   28
  number of point records:    20670652
  number of points by return: 13947228 4610837 1712043 358397 42147
  scale factor x y z:         0.01 0.01 0.01
  offset x y z:               0 0 0
  min x y z:                  582500.00 6712500.00 64.56
  max x y z:                  584999.99 6714999.99 136.59
LASzip compression (version 3.2r2 c2 50000): POINT10 2 GPSTIME11 2
reporting minimum and maximum for all LAS point record entries ...
  X            58250000   58499999
  Y           671250000  671499999
  Z                6456      13659
  intensity          32      61406
  return_number       1          5
  number_of_returns   1          5
  edge_of_flight_line 0          1
  scan_direction_flag 0          1
  classification      1         18
  scan_angle_rank   -19         19
  user_data           0          1
  point_source_ID  1802       1804
  gps_time 222241082.251248 222676871.876191
number of first returns:        13947228
number of intermediate returns: 2110980
number of last returns:         13952166
number of single returns:       9339722
covered area in square units/kilounits: 5923232/5.92
point density: all returns 3.49 last only 2.36 (per square units)
      spacing: all returns 0.54 last only 0.65 (in units)
overview over number of returns of given pulse: 9339722 5797676 4058773 1263967 210514 0 0
histogram of classification of points:
        10888520  unclassified (1)
         9620725  ground (2)
           22695  noise (7)
          138147  water (9)
             565  Reserved for ASPRS Definition (18)
intensity histogram with bin size 256.000000
  bin [0,256) has 1753205
  bin [256,512) has 3009640
  bin [512,768) has 2240861
  bin [768,1024) has 1970696
  bin [1024,1280) has 1610647
  bin [1280,1536) has 1285858
  bin [1536,1792) has 974475
  bin [1792,2048) has 790480
  bin [2048,2304) has 996926
  bin [2304,2560) has 892755
  bin [2560,2816) has 164142
  bin [2816,3072) has 57367
  bin [3072,3328) has 18
         [void]
  bin [10752,11008) has 589317
  bin [11008,11264) has 3760
  bin [11264,11520) has 99653
  bin [11520,11776) has 778739
  bin [11776,12032) has 1393569
  bin [12032,12288) has 1356850
  bin [12288,12544) has 533202
  bin [12544,12800) has 140223
  bin [12800,13056) has 16195
  bin [13056,13312) has 2319
  bin [13312,13568) has 977
  bin [13568,13824) has 765
  bin [13824,14080) has 648
  bin [14080,14336) has 289
  bin [14336,14592) has 513
  bin [14592,14848) has 383
  bin [14848,15104) has 178
  bin [15104,15360) has 526
  bin [15360,15616) has 108
  bin [15616,15872) has 263
  bin [15872,16128) has 289
  bin [16128,16384) has 69
  bin [16384,16640) has 390
  bin [16640,16896) has 51
  bin [16896,17152) has 186
  bin [17152,17408) has 239
  bin [17408,17664) has 169
  bin [17664,17920) has 58
  bin [17920,18176) has 227
  bin [18176,18432) has 169
  bin [18432,18688) has 40
  bin [18688,18944) has 401
  bin [18944,19200) has 30
  bin [19200,19456) has 411
  bin [19456,19712) has 34
  bin [19712,19968) has 34
  bin [19968,20224) has 398
  bin [20224,20480) has 24
  bin [20480,20736) has 108
  bin [20736,20992) has 267
  bin [20992,21248) has 29
  bin [21248,21504) has 318
  bin [21504,21760) has 26
  bin [21760,22016) has 59
  bin [22016,22272) has 184
  bin [22272,22528) has 52
  bin [22528,22784) has 18
  bin [22784,23040) has 116
  bin [23040,23296) has 55
  bin [23296,23552) has 89
  bin [23552,23808) has 250
  bin [23808,24064) has 24
  bin [24064,24320) has 52
  bin [24320,24576) has 14
  bin [24576,24832) has 29
  bin [24832,25088) has 71
  bin [25088,25344) has 74
  bin [25344,25600) has 2
  bin [25600,25856) has 17
  bin [25856,26112) has 2
  bin [26368,26624) has 9
  bin [26624,26880) has 1
  bin [26880,27136) has 1
  bin [27136,27392) has 1
  bin [27392,27648) has 1
  bin [27648,27904) has 3
  bin [28416,28672) has 2
  bin [29184,29440) has 4
  bin [30720,30976) has 1
  bin [30976,31232) has 2
  bin [31232,31488) has 1
  bin [32512,32768) has 1
  bin [36864,37120) has 1
  bin [58368,58624) has 1
  bin [61184,61440) has 1
  average intensity 3625.2240208968733 for 20670652 element(s)

CyArk partners with Google, takes over “Don’t be Evil” Mantra, opens LiDAR Archive

One of our most popular (and controversial) blog articles was “Can You Copyright LiDAR“. It was written after we saw the then chief executive director at CyArk¬†commenting ‚ÄúSweeeet use of CyArk data‚ÄĚ on an article describing the creation of a sugary fudge replica of Guatemala’s¬†Tikal temple promoting a series of sugars by multinational agribusiness Tate & Lyle. Yet just a few months earlier our CEO’s university was instructed to take down his Web pages that – using the same data set – were demonstrating how to realize efficient 3D content delivery across the Web. CyArk¬†told university administrators in an email that he was ‚Äú[‚Ķ]¬†hosting unauthorized content from CyArk¬†[‚Ķ]‚ÄĚ. The full story is here.

Back then, the digital preservation strategy of¬†CyArk¬†was to keep their archaeological scans safe through their partnership with Iron Mountain. In the comment section of “Can You Copyright LiDAR” you can find several entries that are critical of this approach. But that was five years ago.¬†Earlier this year and just after Google removed the¬†“Don’t be Evil” mantra from their code of conduct,¬†CyArk¬†stepped up to¬†take it over and¬†completely changed their tight data control policies. Through their “Open Heritage initiative”¬†CyArk¬†released for the first time their raw¬†LiDAR and imagery with an¬†open license. Here in their own words:

In 2018, CyArk launched the Open Heritage initiative, a
collaboration with Google Arts and Culture to make available
our archive to a broader audience. This was the first time
CyArk has made available primary data sets, including lidar
scans, photogrammetric imagery and corresponding metadata
in a standardized format on a self-serve platform. We are
committed to opening up our archive further as we collect
new data and publishing existing projects where permissions
allow. The data is made available for education, research
and other non-commercial uses via a a Creative Commons
Attribution-NonCommercial 4.0 International License.

This is a HUGE change from the situation in 2013 that resulted in the deletion of our CEO’s Web pages. So we went to download Guatemala’s¬†Tikal temple¬†– the one that got him into trouble back then. It is provided as a single E57 file called ‘Tikal.e57’ with a size of 1074 MB¬†that contains 35,551,759 points in 118 individual scan positions. Using the e572las.exe tool that is part of LAStools¬†we converted this into a single LAZ file ‘Tikal.laz’ with a size of 164 MB.

C:\LAStools\bin>e572las -i c:\data\Tikal\Tikal.e57 ^
                        -o c:\data\Tikal\Tikal.laz

We were not able to find information about the Coordinate Reference System (CRS), but after looking at the coordinate bounding box (see lasinfo report at the end of the article) and the set of projections covering Guatemala, one can make an educated guess that it might be UTM 16 north. Generating a false-colored highest-return 0.5 meter raster with lasgrid and loading it into Google Earth quickly confirms that this is correct.

lasgrid -i c:\data\Tikal\Tikal.laz ^
        -step 0.5 ^
        -highest ^
        -false ^
        -utm 16north ^
        -odix _elev -opng

Now we can laspublish the file with the command line below to create an interactive 3D Web portal using Potree. Unlike five years ago we should now be permitted to create an online portal without the headaches of last time. The CC BY-NC 4.0 license allows to copy and redistribute the material in any medium or format.

laspublish -i c:\data\Tikal\Tikal.laz ^
           -rgb ^
           -utm 16north ^
           -o tikal.html ^
           -title "CyArk's LiDAR Scan of Tikal" ^
           -description "35,551,759 points from 118 individual scans (licensed CC BY-NC 4.0)" ^
           -odir C:\data\Tikal\Tikal -olaz ^
           -overwrite

Below are two screenshots of the online portal that we have just created including some quick distance measurements. This is amazing data. Wow!

Looking at “Templo del Gran Jaguar” from “La Gran Plaza” after taking two measurements.

Overlooking “La Gran Plaza” out of the upper opening of “Templo del Gran Jaguar” with “Templo del las Mascaras” in the back.

We congratulate CyArk to their new Open Heritage initiative and thank them for providing easy access to the Tikal temple LiDAR scans as open data with a useful Creative Commons Attribution-NonCommercial 4.0 International license. Thank you, CyArk, for your contribution to open data and open science. Kudos!

C:\LAStools\bin>lasinfo -i c:\data\Tikal\Tikal.laz
lasinfo (181119) report for 'c:\data\Tikal\Tikal.laz'
reporting all LAS header entries:
  file signature:             'LASF'
  file source ID:             0
  global_encoding:            0
  project ID GUID data 1-4:   00000000-0000-0000-0000-000000000000
  version major.minor:        1.2
  system identifier:          'LAStools (c) by Martin Isenburg'
  generating software:        'e572las.exe (version 180919)'
  file creation day/year:     0/0
  header size:                227
  offset to point data:       227
  number var. length records: 0
  point data format:          2
  point data record length:   26
  number of point records:    35551759
  number of points by return: 35551759 0 0 0 0
  scale factor x y z:         0.001 0.001 0.001
  offset x y z:               220000 1900000 0
  min x y z:                  220854.951 1905881.781 291.967
  max x y z:                  221115.921 1906154.829 341.540
LASzip compression (version 3.2r4 c2 50000): POINT10 2 RGB12 2
reporting minimum and maximum for all LAS point record entries ...
  X              854951    1115921
  Y             5881781    6154829
  Z              291967     341540
  intensity       24832      44800
  return_number       1          1
  number_of_returns   1          1
  edge_of_flight_line 0          0
  scan_direction_flag 0          0
  classification      0          0
  scan_angle_rank     0          0
  user_data           0          0
  point_source_ID     1        118
  Color R 0 65280
        G 0 65280
        B 0 65280
number of first returns:        35551759
number of intermediate returns: 0
number of last returns:         35551759
number of single returns:       35551759
overview over number of returns of given pulse: 35551759 0 0 0 0 0 0
histogram of classification of points:
        35551759  never classified (0)

City of Guadalajara creates first Open LiDAR Portal of Latin America

Small to medium sized LiDAR data sets can easily be published online for exploration and download with laspublish of LAStools, which is an easy-to-use wrapper around the powerful¬†Potree¬†open source software for which rapidlasso GmbH has been a major sponsor. During a workshop on LiDAR processing at CICESE in Ensenada, Mexico we learned that Guadalajara – the city with five “a” in its name – has recently published its LiDAR holdings online for download using an interactive 3D portal based on¬†Potree.

There is a lot more data available in Mexico but only Guadalajara seems to have an interactive download portal at the moment with open LiDAR. Have a look at the map below to get an idea of the LiDAR holdings that are held in the archives of the Instituto Nacional de Estadística y Geografía (INEGI). You can request this data either by filling out this form or by sending an email to atencion.usuarios@inegi.org.mx. You will need to explain the use of the information, but apparently INEGI has a fast response time. I was given the KML files you see below and told that each letter in scale 1: 50,000 is divided into 6 regions (a-f) and each region subdivided into 4 parts. Contact me if you want the KML files or if you can provide further clarification on this indexing scheme and/or the data license.

LiDAR available at the Instituto Nacional de Estadística y Geografía (INEGI)

But back to Guadalajara’s open LiDAR. The tile names become visible when you zoom in closer on the map with the tiling overlay as seen below. An¬†individual tile can easily be downloaded by first clicking so that it becomes highlighted and then pressing the “D” button in the lower left corner. We download the two tiles called ‘F08C04.laz’ and ‘F08C05.laz’ and use lasinfo to determine that their average density is 9.0 and 8.9 last returns per per square meter. This means on average 9 laser pulses were fired at each square meter in those two tiles.

lasinfo -i F08C04.laz -cd
lasinfo -i F08C05.laz -cd

Selecting a tile on the map and pressing the “D” button will download the highlighted tile.

The minimal quality check that we recommend doing for any newly obtained LiDAR data is to verify proper alignment of the flightlines using lasoverlap. For tiles with properly populated ‘point source ID’ fields this can be done using the command line shown below.

lasoverlap -i F08C04.laz F08C05.laz ^
           -min_diff 0.1 -max_diff 0.3 ^
           -odir quality -opng ^
           -cores 2

We notice some slight miss-alignments in the difference image (see other tutorials such as this one for how to interpret the resulting color images). We suggest you follow the steps done there to take a closer look at some of the larger strip-like areas that exhibit some systematic disscolorization (compared to other areas) into overly blueish or reddish tones of with lasview. Overlaying one of the resulting *_diff.png files in the GUI of LAStools makes it easy to pick a suspicious area.

We use the “pick” functionality to view only the building of interest.

Unusual are also the large red and blue areas where some of the taller buildings are. Usually those are just one pixel wide which has to do with the laser of one flightline not being able to see the lower area seen by the laser of the other flightline because the line-of-sight is blocked by the structure. We have a closer look at one of these unusual building colorization by picking the building shown above and viewing it with the different visualization options that are shown in the images below.

No. Those are not the “James Bond movie” kind of lasers that burn holes into the building to get ground returns through several floors. The building facade is covered with glass so that the lasers do not scatter photons when they hit the side of the building. Instead they reflect by the usual rule “incidence angle equals reflection angle” of perfectly specular surfaces and eventually hit the ground next to the building. Some of the photons travel back the same way to the receiver on the plane where they get registered as returns. The LiDAR system has no way to know that the photons did not travel the usual straight path. It only measures the time until the photons return and generates a return at the¬†range¬†corresponding to this time along the direction vector that this laser shot was fired at. If the specular reflection of the photons¬†hits a truck or a tree situated next to to building, then we should find that truck or that tree – mirrored by the glossy surface of the building – on the inside of the building. If you look careful at the “slice” through the building below you may find an example … (-:

Some objects located outside the building are mirrored into the building due to its glossy facade.

Kudos to the City of Guadalajara for becoming – to my knowledge – the first city in Latin America to both open its entire LiDAR holdings and also making it available for download in form of a nice and functional interactive 3D portal.

Scrutinizing LiDAR Data from Leica’s Single Photon Scanner SPL100 (aka SPL99)

We show how simple reordering and clever remapping of single photon LiDAR data can reduce file size by a whopping 50%. We also show that there is at least one¬†Leica’s SPL100 sensor¬†out there that should be called¬†SPL99 because one of its 100 beamlets (the one with beamlet ID 53) does not seem to produce any data … (-:

Closeup on the returns of two beamlet shots colored by beamlet ID from 1 (blue) to 100 (red). Beamlet ID 53 is missing.

Following up on a recent discussion in the LAStools user forum we take a closer look at some of the single photon LiDAR collected with Leica’s SPL100 sensor¬†made available as open data by the USGS¬†in form of LASzip-compressed tiles in LAS 1.4 format of point type 6. This investigation was sparked by the curiosity of what value was stored to the “scanner channel” field that was added to the new point types 6 to 10 in the LAS 1.4 specification.

lasview -i USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP120310_LAS_2018.laz ^
        -copy_scanner_channel_into_point_source ^
        -color_by_flightline

Visualizing this 2 bit number whose value can range from 0 to 3 for the first tile we downloaded¬†resulted in this non-conclusive “magic eye” visualization. What do you see? A sailboat?

Visualizing the “scanner channel” field by mapping its four different values to different colors.

Jason Stoker from the USGS suggested that this is the truncated “beamlet” ID. Leica’s SPL100 sensor¬†uses 100 beamlets rather than one or two laser beams to collect data. Storing the beamlet IDs between 1 and 100 to this 2 bit field that can only hold numbers between 0 and 3 is kind of pointless and should be avoided.¬†LASzip switches prediction contexts based on this field resulting in slower compression speed and lower compression rates. The beamlet ID is also stored in the 8 bit “user data” field, so that we can simply zero the “scanner channel” field. To investigate this further we downloaded these nine tiles from this FTP site of the USGS:

Whenever we download LAZ files we first run laszip with the ‘-check’ option which performs a sanity check to make sure that the files are not truncated or otherwise corrupted. In our case we get nine solid reports of SUCCESS.

laszip -i USGS_LPC_SD_MORiver_Woolpert_B1_*_2018.laz -check

A visual inspection with lasview tells us that there are a number of flightlines in the data.

lasview -i USGS_LPC_SD_MORiver_Woolpert_B1_*_2018.laz ^
        -points 15000000 ^
        -color_by_flightline

We use las2las to extract flightline 2003 and lasinfo to produce a histogram of GPS times which we use in turn to decide on which quarter second of GPS time worth of data we want to extract again with las2las.

las2las -i USGS_LPC_SD_MORiver_Woolpert_B1_2016_*_LAS_2018.laz ^
        -merged ^
        -keep_point_source 2003 ^
        -o USGS_LPC_SD_MORiver_Woolpert_B1_ps_2002.laz

lasinfo -i USGS_LPC_SD_MORiver_Woolpert_B1_ps_2002.laz ^
        -cd ^
        -histo gps_time 1 ^
        -odix _info -otxt

las2las -i USGS_LPC_SD_MORiver_Woolpert_B1_ps_2002.laz ^
        -keep_gps_time 176475495 176475495.25 ^
        -o USGS_LPC_SD_MORiver_Woolpert_B1_gps176475495_quarter.laz

It always helps to give your LAZ files meaningful names in case you find them again two years later or so. We can nicely see the circular scanning pattern¬†Leica’s SPL100 sensor. With lasview we measure that this single flightline has an extent of about 2000 meters on the ground. The lasinfo report shows a pulse density of around 19 last returns per square meter. We then sort the points by GPS time using lassort. This groups together all the returns that are the result of one “shot” of the laser with 100 beamlets as we can nicely see in the las2txt output below. Each group of returns has slightly below 100 points and there is one group every 0.00002 seconds. This means the¬†SPL100 is firing once every 20 microseconds.

lassort -i USGS_LPC_SD_MORiver_Woolpert_B1_gps176475495_quarter.laz ^
        -gps_time ^
        -odix _sorted -olaz

las2txt -i USGS_LPC_SD_MORiver_Woolpert_B1_gps176475495_quarter_sorted.laz ^
        -parse tuxyz ^
        -stdout | more
176475495.000008 4 514408.78 4830989.78 487.79
176475495.000008 9 514410.38 4830987.49 487.70
176475495.000008 47 514411.49 4830987.71 487.70
        [ ... 86 lines deleted ... ]
176475495.000008 39 514408.53 4830991.81 487.80
176475495.000008 50 514407.97 4830991.69 487.80
176475495.000008 16 514409.24 4830991.46 487.85
176475495.000028 55 514413.51 4830985.79 487.61
176475495.000028 97 514411.10 4830990.03 487.74
176475495.000028 72 514411.30 4830989.53 487.74
        [ ... 82 lines deleted ... ]
176475495.000028 45 514410.30 4830986.19 487.70
176475495.000028 3 514409.15 4830987.52 487.73
176475495.000028 96 514411.81 4830985.46 487.67
176475495.000048 66 514411.35 4830985.15 487.67
176475495.000048 83 514411.59 4830984.65 487.61
176475495.000048 64 514413.09 4830983.93 487.61
        [ ... 78 lines deleted ... ]
176475495.000048 4 514407.30 4830984.82 487.70
176475495.000048 34 514408.65 4830983.01 487.70
176475495.000048 21 514408.11 4830982.90 487.70
176475495.000068 13 514408.25 4830981.13 487.66
176475495.000068 92 514410.53 4830984.23 487.68
176475495.000068 44 514407.17 4830980.88 487.67
        [ ... 80 lines deleted ... ]
176475495.000068 76 514408.67 4830984.37 487.71
176475495.000068 47 514409.23 4830980.27 487.67
176475495.000068 87 514412.11 4830981.93 487.61
176475495.000088 97 514408.80 4830982.62 487.70
176475495.000088 33 514407.24 4830980.68 487.64
176475495.000088 30 514407.36 4830981.77 487.68
[ ... ]

Now we can “play back” the returns in acquisition order. We map returns from one group to the same color in¬†lasview¬†with the new ‘-bin_gps_time_into_point_source 0.00002’ option (that will be available in the next LAStools release).¬†For a slower playback we add ‘-steps 5000’. Press the ‘c’ key to switch through the coloring options. Press the ‘s’ key repeatedly to incrementally show the points. To take a step back press <SHIFT>+’s’.

lasview -i USGS_LPC_SD_MORiver_Woolpert_B1_gps176475495_quarter_sorted.laz ^
        -bin_gps_time_into_point_source 0.00002 ^
        -scale_user_data 2.5 ^
        -steps 5000 ^
        -win 1024 384

This slideshow requires JavaScript.

The last image colors the points by the values in the user data field (multiplied by 2.5), which essentially maps the beamlet IDs between 1 and 100 to a rainbow color ramp from blue to red. This tells us how the numbering of the beamlets from 1 to 100 corresponds to their layout in space. The next sequence of images takes a closer look at that.

This slideshow requires JavaScript.

From a compression point of view it makes sense to (1) zero the meaningless scanner channel, (2) order the points by GPS time stamps to groups beamlet returns together, and (3) order the points with the same time stamp by the user data field. The compression gain is enormous with the 9 tiles going from over 3 GB to under 2 GB:

ORIGINAL:
337,156,981 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP120300_LAS_2018.laz
331,801,150 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP120310_LAS_2018.laz
358,928,274 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP120320_LAS_2018.laz
328,597,628 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP130300_LAS_2018.laz
355,997,013 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP130310_LAS_2018.laz
360,403,079 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP130320_LAS_2018.laz
355,399,781 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP140300_LAS_2018.laz
354,523,659 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP140310_LAS_2018.laz
357,248,968 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP140320_LAS_2018.laz
  3,140,056,533 Bytes

IMPROVED:
197,641,087 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP120300_LAS_2018_sorted.laz
194,750,096 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP120310_LAS_2018_sorted.laz
210,013,408 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP120320_LAS_2018_sorted.laz
190,687,275 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP130300_LAS_2018_sorted.laz
206,447,730 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP130310_LAS_2018_sorted.laz
209,580,551 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP130320_LAS_2018_sorted.laz
205,827,197 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP140300_LAS_2018_sorted.laz
203,808,113 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP140310_LAS_2018_sorted.laz
206,789,959 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP140320_LAS_2018_sorted.laz
  1,825,545,416 Bytes

Enumerating the 100 beamlets with a geometrically more coherent order would improve compression even more. Can anyone convince Leica to do this? The simple mapping of beamlet IDs shown below that arranges the beamlets into a zigzag order another huge compression gain of 15 percent. Altogether reordering and remapping lower the compressed file size by a whopping 50 percent.

Beamlet ID mapping table to improve spatial coherence.

168,876,666 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP120300_LAS_2018_mapped_sorted.laz
165,241,508 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP120310_LAS_2018_mapped_sorted.laz
176,524,959 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP120320_LAS_2018_mapped_sorted.laz
163,679,216 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP130300_LAS_2018_mapped_sorted.laz
176,086,559 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP130310_LAS_2018_mapped_sorted.laz
178,909,108 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP130320_LAS_2018_mapped_sorted.laz
174,735,634 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP140300_LAS_2018_mapped_sorted.laz
171,679,105 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP140310_LAS_2018_mapped_sorted.laz
174,997,090 USGS_LPC_SD_MORiver_Woolpert_B1_2016_14TNP140320_LAS_2018_mapped_sorted.laz
  1,550,729,845 Bytes

Once this is done a final space-filling sort into a Hilbert-curve or a Morton-order with lassort or lasoptimize would improve spatial coherence for efficient spatial indexing with lasindex.

Oh yes … the SPL100 was not firing on all cylinders. The beamlet ID 53 that would have mapped to 61 in our table was not present in any of the 9 tiles with 355,047,478 points that we had downloaded as the lasinfo histogram below shows.

lasinfo -i USGS_LPC_SD_MORiver_Woolpert_B1_2016_*_2018.laz -merged -histo user_data 1
lasinfo (180911) report for 9 merged files
reporting all LAS header entries:
  file signature:             'LASF'
  file source ID:             0
  global_encoding:            17
  project ID GUID data 1-4:   194774FA-35FE-4591-D484-010AFA13F6D9
  version major.minor:        1.4
  system identifier:          'Woolpert LAS'
  generating software:        'GeoCue LAS Updater'
  file creation day/year:     332/2017
  header size:                375
  offset to point data:       1376
  number var. length records: 1
  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:                  512000.00 4830000.00 286.43
  max x y z:                  514999.99 4832999.99 866.81
  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: 355047478
  extended number of points by return: 298476060 52480771 3929583 157365 3699 0 0 0 0 0 0 0 0 0 0
variable length header record 1 of 1:
  reserved             0
  user ID              'LASF_Projection'
  record ID            2112
  length after header  943
  description          'OGC WKT Coordinate System'
    WKT OGC COORDINATE SYSTEM:
    COMPD_CS["NAD83(2011) / UTM zone 14N + NAVD88 height - Geoid12B (metre)",PROJCS["NAD83(2011) / UTM zone 14N",GEOGCS["NAD83(2011)",DATUM["NAD83_National_Spat
ial_Reference_System_2011",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","1116"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","
8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","6318"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0]
,PARAMETER["central_meridian",-99],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EP
SG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","6343"]],VERT_CS["NAVD88 height - Geoid12B (metre)",VERT_DATUM["North American Vertica
l Datum 1988",2005,AUTHORITY["EPSG","5103"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","5703"]]]
the header is followed by 4 user-defined bytes
LASzip compression (version 3.1r0 c3 50000): POINT14 3
reporting minimum and maximum for all LAS point record entries ...
  X            51200000   51499999
  Y           483000000  483299999
  Z               28643      86681
  intensity        3139      12341
  return_number       1          5
  number_of_returns   1          5
  edge_of_flight_line 0          0
  scan_direction_flag 0          1
  classification      1         10
  scan_angle_rank  -127        127
  user_data           1        100
  point_source_ID  1061       2005
  gps_time 176475467.194000 176496233.636563
  extended_return_number          1      5
  extended_number_of_returns      1      5
  extended_classification         1     10
  extended_scan_angle        -21167  21167
  extended_scanner_channel        0      3
number of first returns:        298476060
number of intermediate returns: 6282
number of last returns:         355000765
number of single returns:       298435629
overview over extended number of returns of given pulse: 298435629 52515017 3935373 157750 3709 0 0 0 0 0 0 0 0 0 0
histogram of classification of points:
       138382030  unclassified (1)
       207116732  ground (2)
         9233160  noise (7)
          310324  water (9)
            5232  rail (10)
 +-> flagged as withheld:  9233160
 +-> flagged as extended overlap: 226520346
user data histogram with bin size 1.000000
  bin 1 has 3448849
  bin 2 has 3468566
  bin 3 has 3721848
  bin 4 has 3376990
  bin 5 has 3757996
  bin 6 has 3479546
  bin 7 has 3799930
  bin 8 has 3766887
  bin 9 has 3448383
  bin 10 has 3966036
  bin 11 has 3232086
  bin 12 has 3686789
  bin 13 has 3763869
  bin 14 has 3847765
  bin 15 has 3659059
  bin 16 has 3666918
  bin 17 has 3427468
  bin 18 has 3375320
  bin 19 has 3222116
  bin 20 has 3598643
  bin 21 has 3108323
  bin 22 has 3553625
  bin 23 has 3782185
  bin 24 has 3577792
  bin 25 has 3063871
  bin 26 has 3451800
  bin 27 has 3518763
  bin 28 has 3845852
  bin 29 has 3366980
  bin 30 has 3797986
  bin 31 has 3623477
  bin 32 has 3606798
  bin 33 has 3762737
  bin 34 has 3861023
  bin 35 has 3821228
  bin 36 has 3738173
  bin 37 has 3902190
  bin 38 has 3726752
  bin 39 has 3910989
  bin 40 has 3771132
  bin 41 has 3718437
  bin 42 has 3609113
  bin 43 has 3339941
  bin 44 has 3003191
  bin 45 has 3697140
  bin 46 has 2329171
  bin 47 has 3398836
  bin 48 has 3511882
  bin 49 has 3719592
  bin 50 has 2995275
  bin 51 has 3673925
  bin 52 has 3535992
  bin 54 has 3799430
  bin 55 has 3613345
  bin 56 has 3761436
  bin 57 has 3296831
  bin 58 has 3810146
  bin 59 has 3768464
  bin 60 has 3520871
  bin 61 has 3833149
  bin 62 has 3639778
  bin 63 has 3623008
  bin 64 has 3581480
  bin 65 has 3663180
  bin 66 has 3661434
  bin 67 has 3684374
  bin 68 has 3723125
  bin 69 has 3552397
  bin 70 has 3554207
  bin 71 has 3535494
  bin 72 has 3621334
  bin 73 has 3633928
  bin 74 has 3631845
  bin 75 has 3526502
  bin 76 has 3605631
  bin 77 has 3452006
  bin 78 has 3796382
  bin 79 has 3731841
  bin 80 has 3683314
  bin 81 has 3806024
  bin 82 has 3749709
  bin 83 has 3808218
  bin 84 has 3634032
  bin 85 has 3631015
  bin 86 has 3712206
  bin 87 has 3627775
  bin 88 has 3674966
  bin 89 has 3231151
  bin 90 has 3780037
  bin 91 has 3621958
  bin 92 has 3623264
  bin 93 has 3853536
  bin 94 has 3623380
  bin 95 has 3418309
  bin 96 has 3374827
  bin 97 has 3464734
  bin 98 has 3562560
  bin 99 has 3078686
  bin 100 has 3426924