--- a/pkg/models/cross.go	Fri Sep 14 10:18:06 2018 +0200
+++ b/pkg/models/cross.go	Fri Sep 14 11:44:27 2018 +0200
@@ -8,6 +8,7 @@
 	"fmt"
 	"log"
 	"math"
+	"sort"
 	"time"
 )
 
@@ -234,6 +235,107 @@
 	}
 }
 
+func deg2rad(d float64) float64 { return d * math.Pi / 180.0 }
+
+func (cz GeoJSONCoordinateZ) Distance(other GeoJSONCoordinateZ) float64 {
+
+	const EarthRadius = 6378137.0
+
+	dLat := deg2rad(other.Lat - cz.Lat)
+	dLng := math.Abs(deg2rad(other.Lon - cz.Lon))
+
+	if dLng > math.Pi {
+		dLng = 2*math.Pi - dLng
+	}
+
+	x := dLng * math.Cos(deg2rad((cz.Lat+other.Lat)/2.0))
+	return math.Sqrt(dLat*dLat+x*x) * EarthRadius
+}
+
+func (cz GeoJSONCoordinateZ) quant() GeoJSONCoordinateZ {
+	const (
+		xyScale = 1e6
+		zScale  = 1e5
+	)
+	return GeoJSONCoordinateZ{
+		Lon: math.Round(cz.Lon*xyScale) / xyScale,
+		Lat: math.Round(cz.Lat*xyScale) / xyScale,
+		Z:   math.Round(cz.Z*zScale) / zScale,
+	}
+}
+
+func (lcz GeoJSONLineCoordinatesZ) join(other GeoJSONLineCoordinatesZ) GeoJSONLineCoordinatesZ {
+	l1, l2 := len(lcz), len(other)
+	n := make(GeoJSONLineCoordinatesZ, l1+l2-1)
+	copy(n, lcz[:l1-1])
+	n[l1-1] = lcz[l1-1].mid(other[0])
+	copy(n[l1:], other[1:])
+	return n
+}
+
+func (ml GeoJSONMultiLineCoordinatesZ) Join() GeoJSONMultiLineCoordinatesZ {
+
+	if len(ml) < 2 {
+		return ml
+	}
+
+	start := time.Now()
+
+	type value struct {
+		coords GeoJSONLineCoordinatesZ
+		order  int
+	}
+
+	heads := make(map[GeoJSONCoordinateZ]*value)
+
+	for i, coords := range ml {
+		key := coords[len(coords)-1].quant()
+		if v := heads[key]; v != nil {
+			delete(heads, key)
+			v.coords = coords.join(v.coords)
+			heads[v.coords[0].quant()] = v
+		} else {
+			heads[coords[0].quant()] = &value{coords, i}
+		}
+	}
+
+again:
+	for {
+		for k, v1 := range heads {
+			key := v1.coords[len(v1.coords)-1].quant()
+			if k == key { // cycle detected
+				continue
+			}
+			if v2 := heads[key]; v2 != nil {
+				delete(heads, key)
+				v1.coords = v1.coords.join(v2.coords)
+				continue again
+			}
+		}
+		break
+	}
+
+	// To make it deterministic sort in input order.
+	tmp := make([]*value, len(heads))
+	var i int
+	for _, v := range heads {
+		tmp[i] = v
+		i++
+	}
+	sort.Slice(tmp, func(i, j int) bool { return tmp[i].order < tmp[j].order })
+
+	out := make(GeoJSONMultiLineCoordinatesZ, len(tmp))
+	for i, v := range tmp {
+		out[i] = v.coords
+	}
+
+	log.Printf("joining took: %s\n", time.Since(start))
+	log.Printf("segments before/after: %d %d\n", len(ml), len(out))
+
+	return out
+}
+
+/*
 func (ml GeoJSONMultiLineCoordinatesZ) Join() GeoJSONMultiLineCoordinatesZ {
 
 	if len(ml) < 2 {
@@ -278,3 +380,4 @@
 
 	return lists
 }
+*/