lidar

Visualize Lidar Scattered Point Elevation Data

This notebook uses datashader to visualize Lidar elevation data from the Puget Sound Lidar consortium , a source of Lidar data for the Puget Sound region of Washington, U.S.

Setup

Run the download_sample_data.py script to download Lidar data from the S3 datashader examples budget. The script downloads data as a .zip and automatically unzips it to 25 three-column text files with the extension .gnd . and the zip Puget Sound LiDAR consortium and other example data sets.

From your local clone of the datashader repository:

cd examples
conda env create
source activate ds 
python download_sample_data.py

Note on Windows, replace source activate ds with activate ds .

Lidar Elevation Data

Example X,Y,Z scattered point elevation data from the unpacked 7zip files (unpacked as .gnd files)

! head ../data/q47122d2101.gnd
X,Y,Z
1291149.60,181033.64,467.95
1291113.29,181032.53,460.24
1291065.38,181035.74,451.41
1291113.16,181037.32,455.51
1291116.68,181037.42,456.20
1291162.42,181038.90,467.81
1291111.90,181038.15,454.89
1291066.62,181036.73,451.41
1291019.10,181035.20,451.64

The Seattle area example below loads 25 .gnd elevation files like the one above.

In [1]:
import os

from bokeh.models import WMTSTileSource
from dask.distributed import Client
from holoviews.operation.datashader import datashade
from pyproj import Proj, transform
import dask
import dask.dataframe as dd
import geoviews as gv
import glob
import holoviews as hv
import pandas as pd
client = Client()
In [2]:
if not os.path.exists('../data'):
    raise ValueError('Run python download_sample_data.py from the examples directory first')
LIDAR_XYZ_FILES = glob.glob(os.path.join('../data', '*.gnd'))
if not LIDAR_XYZ_FILES:
    raise ValueError('Run python download_sample_data.py from the examples directory first')
LIDAR_XYZ_FILES[:2]
Out[2]:
['../data/q47122d2113.gnd', '../data/q47122d2107.gnd']

Coordinate System Metadata (for this example)

Grid_Coordinate_System_Name : State Plane Coordinate System

State_Plane_Coordinate_System : SPCS_Zone_Identifier Washington North, FIPS 4601

Lambert_Conformal_Conic :

  • Standard_Parallel: 47.500000
  • Standard_Parallel: 48.733333
  • Longitude_of_Central_Meridian: -120.833333
  • Latitude_of_Projection_Origin: 47.000000
  • False_Easting: 1640416.666667
  • False_Northing: 0.000000

http://www.spatialreference.org/ref/esri/102348/

Washington State Plane North - FIPS 4601

In [3]:
washington_state_plane = Proj(init='epsg:2855')   # Washington State Plane North (see metadata above)
web_mercator = Proj(init='epsg:3857')             # Mercator projection EPSG code
In [4]:
FT_2_M = 0.3048    
def convert_coords(ddf):
    lon, lat = transform(washington_state_plane, web_mercator, ddf.X.values * FT_2_M, ddf.Y.values * FT_2_M)
    ddf['meterswest'], ddf['metersnorth'] = lon, lat 
    ddf2 = ddf[['meterswest', 'metersnorth', 'Z']].copy()
    del ddf
    return ddf2

@dask.delayed
def read_gnd(fname):
    return convert_coords(pd.read_csv(fname))

Use web_mercator (from above) to hard-code the bounding box

In [5]:
left, bottom = web_mercator(-122.32, 47.42)
right, top = web_mercator(-122.22, 47.52) 
x_range, y_range = ((left, right), (bottom, top))
In [6]:
df = dd.from_delayed([read_gnd(f) for f in LIDAR_XYZ_FILES])
kdims=['meterswest', 'metersnorth',]
dataset = gv.Dataset(df, kdims=kdims, vdims=['Z'])
shade_defaults = dict(x_range=x_range, y_range=y_range, x_sampling=1, y_sampling=1, width=800, height=455)
tri = hv.Points(dataset, kdims=kdims, vdims=['Z'])
shaded = datashade(tri, **shade_defaults)
df.head()
Out[6]:
meterswest metersnorth Z
0 -1.360602e+07 6.019953e+06 31.46
1 -1.360602e+07 6.019953e+06 31.76
2 -1.360602e+07 6.019953e+06 32.32
3 -1.360602e+07 6.019953e+06 32.35
4 -1.360602e+07 6.019953e+06 32.32

Alternatively we could have done the following dask compute operations to get the bounds of the region:

minn, maxx = df.min().compute(), df.max().compute()
left, bottom = map(float, (minn.meterswest, minn.metersnorth))
right, top = map(float, (maxx.meterswest, maxx.metersnorth))
In [7]:
hv.notebook_extension('bokeh', width=95)

%opts Overlay [width=800 height=800 xaxis=None yaxis=None show_grid=False] 
%opts Shape (fill_color=None line_width=1.5) [apply_ranges=False] 
%opts Points [apply_ranges=False] WMTS (alpha=0.5) NdOverlay [tools=['tap']]