--- a/pkg/imports/sr.go	Sun Oct 06 10:29:43 2019 +0200
+++ b/pkg/imports/sr.go	Sun Oct 06 11:12:32 2019 +0200
@@ -28,13 +28,10 @@
 	"os"
 	"path"
 	"path/filepath"
-	"runtime"
 	"strconv"
 	"strings"
-	"sync"
 	"time"
 
-	"github.com/fogleman/contourmap"
 	shp "github.com/jonas-p/go-shp"
 
 	"gemma.intevation.de/gemma/pkg/common"
@@ -925,124 +922,7 @@
 			time.Since(total))
 	}()
 
-	min, max := tree.Min, tree.Max
-
-	width := max.X - min.X
-	height := max.Y - min.Y
-
-	feedback.Info("Width/Height: %.2f / %.2f", width, height)
-
-	xcells := int(math.Ceil(width / isoCellSize))
-	ycells := int(math.Ceil(height / isoCellSize))
-
-	feedback.Info("Raster size: (%d, %d)", xcells, ycells)
-
-	// Add border for closing
-	raster := make([]float64, (xcells+2)*(ycells+2))
-
-	// prefill for no data
-	const nodata = -math.MaxFloat64
-	for i := range raster {
-		raster[i] = nodata
-	}
-
-	// rasterize the height model
-
-	var wg sync.WaitGroup
-
-	rows := make(chan int)
-
-	rasterRow := func() {
-		defer wg.Done()
-		for i := range rows {
-			pos := (i+1)*(xcells+2) + 1
-			row := raster[pos : pos+xcells]
-			py := min.Y + float64(i)*isoCellSize + isoCellSize/2
-			px := min.X + isoCellSize/2
-			for j := range row {
-				v, ok := tree.Value(px, py)
-				if ok {
-					row[j] = v
-				}
-				px += isoCellSize
-			}
-		}
-	}
-
-	start := time.Now()
-
-	for n := runtime.NumCPU(); n >= 1; n-- {
-		wg.Add(1)
-		go rasterRow()
-	}
-
-	for i := 0; i < ycells; i++ {
-		rows <- i
-	}
-	close(rows)
-
-	wg.Wait()
-
-	feedback.Info("Rasterizing took %v.", time.Since(start))
-
-	feedback.Info("Trace iso areas.")
-
-	start = time.Now()
-
-	tracer := contourmap.FromFloat64s(xcells+2, ycells+2, raster)
-
-	areas := make([]wkb.MultiPolygonGeom, len(heights))
-
-	// TODO: Check if this correct!
-	reprojX := common.Linear(1, min.X, float64(xcells+1), max.X)
-	reprojY := common.Linear(1, min.Y, float64(ycells+1), max.Y)
-
-	cnts := make(chan int)
-
-	doContours := func() {
-		defer wg.Done()
-		for hIdx := range cnts {
-			c := tracer.Contours(heights[hIdx])
-
-			if len(c) == 0 {
-				continue
-			}
-
-			// We need to bring it back to the
-			// none raster coordinate system.
-			a := make(wkb.MultiPolygonGeom, len(c))
-
-			for i, pl := range c {
-				shell := make(wkb.LinearRingGeom, len(pl))
-				for j, pt := range pl {
-					shell[j] = wkb.PointGeom{
-						X: reprojX(pt.X),
-						Y: reprojY(pt.Y),
-					}
-				}
-				a[i] = wkb.PolygonGeom{shell}
-			}
-
-			areas[hIdx] = a
-		}
-	}
-
-	for n := runtime.NumCPU(); n >= 1; n-- {
-		wg.Add(1)
-		go doContours()
-	}
-
-	for i := range heights {
-		cnts <- i
-	}
-	close(cnts)
-
-	wg.Wait()
-
-	feedback.Info("Tracing areas took %v.", time.Since(start))
-
-	// Raster and tracer are not needed any more.
-	raster, tracer = nil, nil
+	areas := tree.TraceAreas(heights, isoCellSize)
 
 	return storeAreas(
 		ctx, tx, feedback,

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pkg/octree/areas.go	Sun Oct 06 11:12:32 2019 +0200
@@ -0,0 +1,146 @@
+// This is Free Software under GNU Affero General Public License v >= 3.0
+// without warranty, see README.md and license for details.
+//
+// SPDX-License-Identifier: AGPL-3.0-or-later
+// License-Filename: LICENSES/AGPL-3.0.txt
+//
+// Copyright (C) 2018 by via donau
+//   – Österreichische Wasserstraßen-Gesellschaft mbH
+// Software engineering by Intevation GmbH
+//
+// Author(s):
+//  * Sascha L. Teichmann <sascha.teichmann@intevation.de>
+
+package octree
+
+import (
+	"log"
+	"math"
+	"runtime"
+	"sync"
+	"time"
+
+	"github.com/fogleman/contourmap"
+
+	"gemma.intevation.de/gemma/pkg/common"
+	"gemma.intevation.de/gemma/pkg/wkb"
+)
+
+func (tree *Tree) TraceAreas(
+	heights []float64,
+	cellSize float64,
+) []wkb.MultiPolygonGeom {
+	min, max := tree.Min, tree.Max
+
+	width := max.X - min.X
+	height := max.Y - min.Y
+
+	log.Printf("info: Width/Height: %.2f / %.2f\n", width, height)
+
+	xcells := int(math.Ceil(width / cellSize))
+	ycells := int(math.Ceil(height / cellSize))
+
+	log.Printf("info: Raster size: (%d, %d)\n", xcells, ycells)
+
+	start := time.Now()
+
+	// Add border for closing
+	raster := make([]float64, (xcells+2)*(ycells+2))
+
+	// prefill for no data
+	const nodata = -math.MaxFloat64
+	for i := range raster {
+		raster[i] = nodata
+	}
+
+	// rasterize the height model
+
+	var wg sync.WaitGroup
+
+	rows := make(chan int)
+
+	rasterRow := func() {
+		defer wg.Done()
+		for i := range rows {
+			pos := (i+1)*(xcells+2) + 1
+			row := raster[pos : pos+xcells]
+			py := min.Y + float64(i)*cellSize + cellSize/2
+			px := min.X + cellSize/2
+			for j := range row {
+				v, ok := tree.Value(px, py)
+				if ok {
+					row[j] = v
+				}
+				px += cellSize
+			}
+		}
+	}
+
+	for n := runtime.NumCPU(); n >= 1; n-- {
+		wg.Add(1)
+		go rasterRow()
+	}
+
+	for i := 0; i < ycells; i++ {
+		rows <- i
+	}
+	close(rows)
+
+	wg.Wait()
+	log.Printf("info: Rastering took %v\n", time.Since(start))
+
+	start = time.Now()
+
+	tracer := contourmap.FromFloat64s(xcells+2, ycells+2, raster)
+
+	areas := make([]wkb.MultiPolygonGeom, len(heights))
+
+	// TODO: Check if this correct!
+	reprojX := common.Linear(1, min.X, float64(xcells+1), max.X)
+	reprojY := common.Linear(1, min.Y, float64(ycells+1), max.Y)
+
+	cnts := make(chan int)
+
+	doContours := func() {
+		defer wg.Done()
+		for hIdx := range cnts {
+			c := tracer.Contours(heights[hIdx])
+
+			if len(c) == 0 {
+				continue
+			}
+
+			// We need to bring it back to the
+			// none raster coordinate system.
+			a := make(wkb.MultiPolygonGeom, len(c))
+
+			for i, pl := range c {
+				shell := make(wkb.LinearRingGeom, len(pl))
+				for j, pt := range pl {
+					shell[j] = wkb.PointGeom{
+						X: reprojX(pt.X),
+						Y: reprojY(pt.Y),
+					}
+				}
+				a[i] = wkb.PolygonGeom{shell}
+			}
+
+			areas[hIdx] = a
+		}
+	}
+
+	for n := runtime.NumCPU(); n >= 1; n-- {
+		wg.Add(1)
+		go doContours()
+	}
+
+	for i := range heights {
+		cnts <- i
+	}
+	close(cnts)
+
+	wg.Wait()
+	log.Printf("info: Tracing areas took %v\n", time.Since(start))
+
+	return areas
+}