```
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

```
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()

```
```

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

```
```

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)

```
```

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)

```
```

```
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]:
```

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))

```
```

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))

```
```

```
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)

```
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]:
```