Downloading and processing LiDAR data and aerial photograph

This notebook gives a toturial which treats the subjects of getting LiDAR data, and give it a color using aerial photgraphs.


In [1]:
# Import packages

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import owslib.wfs
import owslib.wms
import json
import shapely.geometry
from matplotlib import collections
import timeit
import liblas

LiDAR

The lidar data from the dutch ahn2 dataset is split up into tiles of 5 by 6.25 kilometers. To find out which tile we need to download in order to get the data we want, we can use the WFS. This is extensively discribed in the ahn2 notebook available in the github/notebooks repositry.


In [2]:
# Enter WFS properties
wfsurl = 'http://geodata.nationaalgeoregister.nl/ahn2/wfs'
wfs = owslib.wfs.WebFeatureService(wfsurl, version="2.0.0")
wfslayer = wfs.contents['ahn2:ahn2_bladindex']

# Get boxes from WFS
f = wfs.getfeature(typename=[wfslayer.id], outputFormat="json")
data = json.load(f)

shapes = []
tiles = []

for feature in data['features']:
    shapes.append(shapely.geometry.asShape(feature['geometry'])[0])
    tiles.append(shapes[-1].exterior.coords[:])
print len(tiles)
# Plot boxes using a collection
col = collections.LineCollection(tiles, colors = 'k')
fig, ax  = plt.subplots()
fig.set_size_inches(10,10)
ax.add_collection(col, autolim = True)
ax.set_aspect('equal')
ax.autoscale_view()


1372

The coordinates are in the Dutch RijksDriehoekstelsel.


In [3]:
# Center of Delft in RD-coordinates:
x = 84468 
y = 447557

# Create a shapely point
p = shapely.geometry.Point(x, y)

# Check in which box the point lies
for i in range(len(shapes)):
    if p.within(shapes[i]):
        bladnr = data['features'][i]['properties']['bladnr']
        delftshape = shapes[i]
        print 'The point is in box {}.'.format(bladnr)
        break


The point is in box 37en1.

Now we need to download the lidar tile 37en1. We can do this by surfing to: http://geodata.nationaalgeoregister.nl/ahn2/atom/ahn2_uitgefilterd.xml and http://geodata.nationaalgeoregister.nl/ahn2/atom/ahn2_gefilterd.xml and search for our box. There are two links since there are two datasets. One, 'gefilterd', contains the point cloud with all the objects (houses, trees, etcetera) removed. The other 'uitgefilterd' contains these objects. We will need both of them. Let's download 37en1 for both. This might take a wile because the total is almost 1GB.

When we have downloaded and unpacked the files, we get two laz files. The laz files contain the x, y and z coordinates of the point cloud. Let's see what's inside the file


In [4]:
fileloc = r'D:\rongen\Documents\IPython Notebooks\\'
files = [fileloc+'g37en1.laz', fileloc+'u37en1.laz']

count = 0
lasfiles = []
for f in files:
    lasfiles.append(liblas.file.File(f, mode = 'r'))
    count += lasfiles[-1].get_header().get_count()
    
print 'The number of points in the cloud is {}.'.format(count)
print 'The number of points per square meter is {}.'.format(count/delftshape.area)


The number of points in the cloud is 606665728.
The number of points per square meter is 19.413303296.

That are over 600 million points! Since this number is to high for the visualization we aim te make, we will reduce our area to an area of 500 x 500 m in the center of Delft. This will give us 5 million points which should be processable.

The problem with getting our five million points is that we cannot know the location of the point on beforehand. Also the points have to be extracted one by one, so building in an if-statement to check the location will greatly increase the computation time. Extracting all point will probably take more than an hour and some memory problems so we wont do that.

We do however know that the point are somewhat grouped in subtiles within the las files (Trust me on that one). So if we take a sample of 10000 points, and check if they are in our range plus some buffer, we roughly know which part of the las file is in our range. This way our computation time will be under 10 minutes, and we will cover most points. Let's select the bounds and extract a sample. This will probably take about two minutes.


In [5]:
bounds      = (x-250, y-250, x+250, y+250)
outerbounds = (x-500, y-500, x+500, y+500)

In [6]:
# Create a function to get the index and coordinates from the laz file
def las2rows(lasfile, start, stop, step):
    
    for i in range(start, stop, step):
        p = lasfile.read(i)
        yield i
        yield p.x
        yield p.y
        yield p.z
        
starttime = timeit.default_timer()

xyz = []

# Get sample from las file
for las in lasfiles:
    npoints = las.get_header().get_count()

    start = 10000
    stop = npoints+1
    step = 20000
    
    xyz.append(np.fromiter(las2rows(las, start, stop, step), np.float32, npoints/step*4).reshape(npoints/step, 4))

print 'Extracting time: {} s.'.format(timeit.default_timer()-starttime)


Extracting time: 138.624283207 s.

Now filter out the point within our bounds and make a scatterplot to see if the selection is succesfull


In [7]:
for i in range(len(lasfiles)):
    index = ((xyz[i][:,1] > outerbounds[0]) *
             (xyz[i][:,1] < outerbounds[2]) *
             (xyz[i][:,2] > outerbounds[1]) *
             (xyz[i][:,2] < outerbounds[3]))

    xyz[i] = xyz[i][index,:]

    plt.scatter(xyz[i][:,1], xyz[i][:,2], marker='.', linewidths = 0, vmax=5, vmin=-5)
plt.plot([bounds[0], bounds[0], bounds[2], bounds[2], bounds[0]], [bounds[1], bounds[3], bounds[3], bounds[1], bounds[1]], 'b-')


Out[7]:
[<matplotlib.lines.Line2D at 0xbfb9a90>]

Around the selected points we want the 10000 points before and 10000 points after.


In [8]:
subranges = []
for i in range(len(lasfiles)):
    subrange = []
    for j in xyz[i]:
        subrange.append([int(j[0]-step/2.), int(j[0]+step/2.)])
    subranges.append(subrange)

And get the points from the subrange, I suggest you take a break:


In [9]:
# Create another function to retrieve points. Note that this one is different since it does not return the index.
def las2rows(lasfile, start, stop):
    
    for i in xrange(start, stop, 1):
        p = lasfile.read(i)
        yield p.x
        yield p.y
        yield p.z

starttime = timeit.default_timer()
xyz = []
for i in range(len(lasfiles)):
    for subr in subranges[i]:
        npoints = subr[1] - subr[0]
        xyz.append(np.fromiter(las2rows(lasfiles[i], subr[0], subr[1]), np.float32, npoints*3).reshape(npoints, 3))
    
print 'Extracting time: {} s.'.format(timeit.default_timer()-starttime)

xyz = np.vstack(xyz).astype(np.float32)
print 'There are {} points in our extracted set.'.format(len(xyz))


Extracting time: 381.780448828 s.
There are 18480000 points in our extracted set.

And finally filter out the points we need again:


In [10]:
index = ((xyz[:,0] > bounds[0]) *
         (xyz[:,0] < bounds[2]) *
         (xyz[:,1] > bounds[1]) *
         (xyz[:,1] < bounds[3]))

xyz = xyz[index,:]
print 'There are {} points in our bounded set.'.format(len(xyz))


There are 4266958 points in our bounded set.

Aerial photographs

We will also need the color values to make a nice visualisation. For this we will use luchtfotos which we will get via WMS. The maximum size for our request is about 16 MP. The accuracy of the aerial photgraphs is 0.5 meter. We choose 1250 x 1250 pizels for our 500 x 500 area, which is the maximum reolution we can request. The reason for this has probably something to do with the accuracy.


In [11]:
# Get WMS layer
wmsurl   = 'http://geodata1.nationaalgeoregister.nl/luchtfoto/wms'
wms      = owslib.wms.WebMapService(wmsurl)
wmslayer = ['luchtfoto_png']            

# Download aerial photograps
for i in range(1):
    f = wms.getmap(layers = wmslayer,
                   srs    = 'EPSG:28992',
                   bbox   = bounds,
                   size   = (1250, 1250),
                   format = 'image/png')
    rgb = (plt.imread(f)*255).astype(np.uint8)

We increase the saturation a bit to make the colors a little less gray and plot the image to see the result.


In [12]:
# Increase saturation
from skimage import color
rgb = color.rgb2hsv(rgb)
rgb[:,:,1] *= 1.5
rgb = color.hsv2rgb(rgb)

# Plot figure
fig, ax  = plt.subplots()
fig.set_size_inches(10,10)
ax.imshow(rgb)


Out[12]:
<matplotlib.image.AxesImage at 0x3fca4a8>