--- a/pkg/imports/sr.go	Wed Oct 02 16:43:34 2019 +0200
+++ b/pkg/imports/sr.go	Wed Oct 02 18:35:09 2019 +0200
@@ -28,15 +28,19 @@
 	"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"
 	"gemma.intevation.de/gemma/pkg/models"
 	"gemma.intevation.de/gemma/pkg/octree"
+	"gemma.intevation.de/gemma/pkg/wkb"
 )
 
 // SoundingResult is a Job to import sounding reults
@@ -70,9 +74,17 @@
 )
 
 const (
+	// pointsPerSquareMeter is the average number of points
+	// when generating a artifical height model for single beam scans.
 	pointsPerSquareMeter = 2
 )
 
+const (
+	// isoCellSize is the side length of a raster cell when tracing
+	// iso areas.
+	isoCellSize = 0.5
+)
+
 // SRJobKind is the unique name of this import job type.
 const SRJobKind JobKind = "sr"
 
@@ -177,6 +189,31 @@
 FROM waterway.sounding_results sr
 WHERE id = $1
 `
+	insertIsoAreasSQL = `
+INSERT INTO waterway.sounding_results_iso_areas (
+  sounding_result_id,
+  height,
+  areas
+)
+SELECT
+  $1,
+  $2,
+  ST_Transform(
+    ST_Multi(
+      ST_CollectionExtract(
+        ST_SimplifyPreserveTopology(
+          ST_Multi(ST_Collectionextract(
+            ST_MakeValid(ST_GeomFromWKB($4, $3::integer)), 2)),
+          $5
+        ),
+        2
+      )
+    ),
+    4326
+  )
+FROM waterway.sounding_results sr
+WHERE id = $1
+`
 
 	selectGaugeLDCSQL = `
 SELECT
@@ -619,15 +656,10 @@
 		return nil, err
 	}
 	feedback.Info("Storing octree index took %s.", time.Since(start))
-	feedback.Info("Generate contour lines")
-
-	start = time.Now()
-	err = generateContours(ctx, tx, feedback, builder.Tree(), id)
+	err = generateIsos(ctx, tx, feedback, builder.Tree(), id)
 	if err != nil {
 		return nil, err
 	}
-	feedback.Info("Generating contour lines took %s.",
-		time.Since(start))
 
 	// Store for potential later removal.
 	if err = track(ctx, tx, importID, "waterway.sounding_results", id); err != nil {
@@ -827,23 +859,19 @@
 	return shapeToPolygon(s)
 }
 
-func generateContours(
+func generateIsos(
 	ctx context.Context,
 	tx *sql.Tx,
 	feedback Feedback,
 	tree *octree.Tree,
 	id int64,
 ) error {
-	stmt, err := tx.PrepareContext(ctx, insertContourSQL)
-	if err != nil {
-		return err
-	}
-	defer stmt.Close()
 
 	heights, err := octree.LoadClassBreaks(
 		ctx, tx,
 		"morphology_classbreaks",
 	)
+
 	if err != nil {
 		feedback.Warn("Loading class breaks failed: %v", err)
 		feedback.Info("Using default class breaks")
@@ -871,6 +899,189 @@
 
 	heights = common.DedupFloat64s(heights)
 
+	if err := generateContours(ctx, tx, feedback, tree, heights, id); err != nil {
+		return err
+	}
+
+	return generateIsoAreas(ctx, tx, feedback, tree, heights, id)
+}
+
+func generateIsoAreas(
+	ctx context.Context,
+	tx *sql.Tx,
+	feedback Feedback,
+	tree *octree.Tree,
+	heights []float64,
+	id int64,
+) error {
+	feedback.Info("Generate iso areas")
+	total := time.Now()
+	defer func() {
+		feedback.Info("Generating iso areas took %s.",
+			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))
+
+	// TODO: Implement me!
+	feedback.Info("Trace iso areas")
+
+	start = time.Now()
+
+	tracer := contourmap.FromFloat64s(xcells+2, ycells+2, raster)
+
+	areas := make([]wkb.MultiPolygonGeom, len(heights))
+
+	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])
+			a := make(wkb.MultiPolygonGeom, len(c))
+
+			for i, pl := range c {
+
+				shell := make(wkb.LinearRingGeom, len(pl))
+				for j, pt := range pl {
+					x := reprojX(pt.X)
+					y := reprojY(pt.Y)
+					shell[j] = wkb.PointGeom{X: x, Y: 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
+
+	return storeAreas(ctx, tx, feedback, areas, heights)
+}
+
+func storeAreas(
+	ctx context.Context,
+	tx *sql.Tx,
+	feedback Feedback,
+	areas []wkb.MultiPolygonGeom,
+	heights []float64,
+) error {
+	feedback.Info("Storing iso areas")
+	total := time.Now()
+	defer func() {
+		feedback.Info("Storing iso areas took %v",
+			time.Since(total))
+	}()
+
+	stmt, err := tx.PrepareContext(ctx, insertIsoAreasSQL)
+	if err != nil {
+		return err
+	}
+	defer stmt.Close()
+
+	// TODO: Implement me!
+
+	return nil
+}
+
+func generateContours(
+	ctx context.Context,
+	tx *sql.Tx,
+	feedback Feedback,
+	tree *octree.Tree,
+	heights []float64,
+	id int64,
+) error {
+	feedback.Info("Generate contour lines")
+	start := time.Now()
+	defer func() {
+		feedback.Info("Generating contour lines took %v",
+			time.Since(start))
+	}()
+
+	stmt, err := tx.PrepareContext(ctx, insertContourSQL)
+	if err != nil {
+		return err
+	}
+	defer stmt.Close()
+
 	octree.DoContours(tree, heights, func(res *octree.ContourResult) {
 		if err == nil && len(res.Lines) > 0 {
 			_, err = stmt.ExecContext(