USGS LiDAR point cloud of New York

retrieve, process and view the LiDAR data of New York city

December 25, 2018

In Winter 2012, Storm Sandy hit New York, claimed 147 lives, and caused $60 billion of damage. Subsequent to the storm, the US Geological Survey conducted a LiDAR mission to map the region in 3D for the purpose of post-disaster assessment. The mission completed in 2014. The data portion covering New York City is publically available on the NOAA's website as "2014 USGS CMGP LiDAR: Post Sandy - New York City". This post will provide some insight into the dataset and guides you through a simple procedure to retrieve, process, and view the data. No commercial software will be needed if you want to follow along.

Storm Sandy (1) (2) USGS 2013-2014 LiDAR

Data

The 2014 USGS CMGP is about 25.55 GB, containing 4,583,325,042 points, and covering approximately 787 km2. The project aimed for a point density of 2 points/m2. The actual density varies quite a lot and is as high as 9 points/m2 in some areas. Compared to our 2015 LiDAR data collections in Dublin, this dataset is more than an order of magnitude sparser as a result of the higher flying altitude, higher flying speed, slower pulse frequency, and lower level of overlapping as outlined in the table below.

Compare NY LiDAR and Dublin LiDAR

2014 New York 2015 Dublin
Flight altitude (m) 2,286 300
Flying speed (knots) 173 50
Pulse repetition frequency (kHz) 239 400
Side lap 25% 2 x 67%
Number of points 4,583,325,042 1,420,984,613
Size in LAZ format (GB) 25.55 7.18
Coverage area (km2) 787 4

Retrieve, process and view


Retrieve

25.55 gigabytes of data does not seem to be a large amount. However, not many software can handle that much of points. The number of points matters more than the disk consumption as far as point clouds are concerned. With point clouds at that size, we often need to subset the data before doing anything with the data.

The 2014 USGS CMGP dataset is organized as a set of disjoint tiles as showed in the map below. Each tile is stored as an LAZ file. LAZ is a compressed version of the LAS format, the standard file format for airborne LiDAR data developed by the American Society of Photogrammetry and Remote Sensing. The dataset on the NOAA's website comes with a tile index in the ERSI's shapefile format. For convenience, I have transformed the tile index to an interactive map as showed below.

LAS format LAZ

To use the map to retrieve data:

  1. Click the button at the top right corner of the map to maximize it.
  2. Use the search tool (the magnifying glass icon at the top left corner of the map) to search for an address or a point of interest of your choice.
  3. Click on the tiles intersecting your area of interest. The URLs to download the corresponding files will pop up.
  4. Download the files to somewhere you can remember.



Process

The files downloaded from the USGS website are in a relatively atypical format. The elevation values of the data are in meters above the Mean Sea Level. More precisely, it is in an orthometric vertical datum NAVD88 using GEOID12B. The horizontal datum of the dataset is UTM, Zone 18, North American Datum of 1983 (2011) [EPSG:4269]. In that projection, the x- and y-values are in degrees. Due to the unit mismatch, many point cloud visualization software fail to display the data.

Geodetic datum

Here is how the coordinates look like when converted to a text format.

    
        x(deg),      y(deg),     z(meter)
        -73.7260988, 40.7654425, 57.277
        -73.7260977, 40.7654422, 57.356
        -73.7260972, 40.7654422, 57.286
        -73.7260962, 40.7654421, 57.347
        -73.7260955, 40.7654421, 57.257
        -73.7258715, 40.7654410, 57.228
        -73.7258774, 40.7654415, 57.108
        -73.7258829, 40.7654419, 57.117
        -73.7258884, 40.7654421, 57.108
        -73.7258937, 40.7654426, 57.167
    

Here is what you would see when plotting the x (deg), y (deg), z (m) coordinates directly.

To properly plot the point data, you need to reproject it so that the x-, y-, and z-values are in the same unit. If you are in the metric system like me. Use the following libLAS command. libLAS is a free and open-source C++ library for LAS data processing. In addition to reprojecting the data to EPSG:32118, the command also modify the offset and scale parameters to match the new projection.

libLAS


las2las  -v [INPUT_FILE] [OUTPUT_FILE] --offset 277000,33000,0 --scale .01 .01 .01 --a_srs EPSG:4269 --t_srs EPSG:32118

If you prefer the US-foot system, you can use the projection EPSG:4269, and the libLAS command is:


las2las -v [INPUT_FILE] [OUTPUT_FILE] --offset 970000,180000,0 --scale .01 .01 .01 --a_srs EPSG:4269 --t_srs EPSG:2263

If you prefer using FME instead of libLAS, the equivalent workspaces can be downloaded at the link below.

Equivalent FME workspaces

After the re-projection, you can use any point cloud data viewer (e.g. CloudCompare, Fugro Viewer) to display the data. The point cloud should no longer look 'weird'. You can do whatever you could imagine with the 3D point cloud. Here are some suggestions:

  • Measure tree heights, building heights
  • The point cloud is discrete, can you convert that into something continuous [e.g. a DSM, or a TIN]
  • Can you pour some water onto the model and simulate flood inundation?


A not-so-fun fact

The post was written during the US government shutdown in Dec 2018. Attempt to access the NOAA website brought me to the pages below. Luckily, I had a copy of the dataset stored locally.

hutdown-1

hutdown-1