--- a/pkg/imports/sr.go	Tue Jun 04 18:11:23 2019 +0200
+++ b/pkg/imports/sr.go	Tue Jun 04 18:37:15 2019 +0200
@@ -28,11 +28,13 @@
 	"os"
 	"path"
 	"path/filepath"
+	"sort"
 	"strconv"
 	"strings"
 	"time"
 
 	shp "github.com/jonas-p/go-shp"
+	"gonum.org/v1/gonum/stat"
 
 	"gemma.intevation.de/gemma/pkg/common"
 	"gemma.intevation.de/gemma/pkg/misc"
@@ -176,11 +178,12 @@
 `
 
 	reprojectPointsSingleBeamSQL = `
-SELECT ST_AsBinary(
-  ST_Transform(
-    ST_GeomFromWKB($1, $2::integer),
-	best_utm(CAST(ST_GeomFromWKB($1, $2::integer) AS geography))),
-	best_utm(CAST(ST_GeomFromWKB($1, $2::integer) AS geography))
+SELECT
+  ST_AsBinary(
+    ST_Transform(
+	  ST_GeomFromWKB($1, $2::integer),
+	  best_utm(CAST(ST_Transform(ST_GeomFromWKB($1, $2::integer), 4326) AS geography)))),
+  best_utm(CAST(ST_Transform(ST_GeomFromWKB($1, $2::integer), 4326) AS geography))
 `
 )
 
@@ -205,6 +208,7 @@
 	defer z.Close()
 
 	feedback.Info("Looking for 'meta.json'")
+
 	mf := common.FindInZIP(z, "meta.json")
 	if mf == nil && !sr.completeOverride() {
 		return nil, errors.New("Cannot find 'meta.json'")
@@ -342,8 +346,57 @@
 		return nil, err
 	}
 
+	if err := xyz.FromWKB(reproj); err != nil {
+		return nil, err
+	}
+
 	feedback.Info("Reprojecting points to EPSG %d took %v.",
 		epsg, time.Since(start))
+	feedback.Info("Number of reprojected points: %d", len(xyz))
+
+	start = time.Now()
+	tri, err := octree.Triangulate(xyz)
+	if err != nil {
+		return nil, err
+	}
+	feedback.Info("Triangulation took %v.", time.Since(start))
+	feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
+
+	if boundary == nil {
+		feedback.Info("No boundary given.")
+		feedback.Info("Eliminate triangles with long edges.")
+
+		lens := make([]float64, len(tri.Triangles)/3)
+		for i := 0; i < len(tri.Triangles); i += 3 {
+			vs := tri.Triangles[i : i+3]
+			var max float64
+			for j := 0; j < 3; j++ {
+				k := (j + 1) % 3
+				if l := tri.Points[vs[j]].SquaredDistance2D(tri.Points[vs[k]]); l > max {
+					max = l
+				}
+			}
+			lens[i/3] = max
+		}
+		qs := make([]float64, len(lens))
+		copy(qs, lens)
+		sort.Float64s(qs)
+		q := stat.Quantile(0.98, stat.LinInterp, qs, nil)
+		feedback.Info("98%% quantile of the edge lengths: %.2f\n", math.Sqrt(q))
+
+		ts := tri.Triangles
+
+		for i := len(lens) - 1; i >= 0; i-- {
+			if lens[i] >= q {
+				if i != len(lens)-1 {
+					copy(ts[i*3:], ts[(i+1)*3:])
+				}
+				ts = ts[:len(ts)-3]
+			}
+		}
+		feedback.Info("Remaining triangles: %d (%.2f%%)\n",
+			len(ts)/3, float64(len(ts))/float64(len(tri.Triangles))*100.0)
+	}
 
 	// TODO: Implement me!
 	return nil, errors.New("Not implemented, yet!")