Mercurial > gemma

comparison pkg/imports/sr.go @ 4570:4b3a298b94f8 iso-areas

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Moved area tracing to octree package.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Sun, 06 Oct 2019 11:12:32 +0200
parents 2c49a293f275
children c657dec6b0fa
comparison
equal deleted inserted replaced
4569:efb7851f54dd 4570:4b3a298b94f8
26 "io" 26 "io"
27 "math" 27 "math"
28 "os" 28 "os"
29 "path" 29 "path"
30 "path/filepath" 30 "path/filepath"
31 "runtime"
32 "strconv" 31 "strconv"
33 "strings" 32 "strings"
34 "sync"
35 "time" 33 "time"
36 34
37 "github.com/fogleman/contourmap"
38 shp "github.com/jonas-p/go-shp" 35 shp "github.com/jonas-p/go-shp"
39 36
40 "gemma.intevation.de/gemma/pkg/common" 37 "gemma.intevation.de/gemma/pkg/common"
41 "gemma.intevation.de/gemma/pkg/models" 38 "gemma.intevation.de/gemma/pkg/models"
42 "gemma.intevation.de/gemma/pkg/octree" 39 "gemma.intevation.de/gemma/pkg/octree"
923 defer func() { 920 defer func() {
924 feedback.Info("Generating iso areas took %s.", 921 feedback.Info("Generating iso areas took %s.",
925 time.Since(total)) 922 time.Since(total))
926 }() 923 }()
927 924
928 min, max := tree.Min, tree.Max 925 areas := tree.TraceAreas(heights, isoCellSize)
929
930 width := max.X - min.X
931 height := max.Y - min.Y
932
933 feedback.Info("Width/Height: %.2f / %.2f", width, height)
934
935 xcells := int(math.Ceil(width / isoCellSize))
936 ycells := int(math.Ceil(height / isoCellSize))
937
938 feedback.Info("Raster size: (%d, %d)", xcells, ycells)
939
940 // Add border for closing
941 raster := make([]float64, (xcells+2)*(ycells+2))
942
943 // prefill for no data
944 const nodata = -math.MaxFloat64
945 for i := range raster {
946 raster[i] = nodata
947 }
948
949 // rasterize the height model
950
951 var wg sync.WaitGroup
952
953 rows := make(chan int)
954
955 rasterRow := func() {
956 defer wg.Done()
957 for i := range rows {
958 pos := (i+1)*(xcells+2) + 1
959 row := raster[pos : pos+xcells]
960 py := min.Y + float64(i)*isoCellSize + isoCellSize/2
961 px := min.X + isoCellSize/2
962 for j := range row {
963 v, ok := tree.Value(px, py)
964 if ok {
965 row[j] = v
966 }
967 px += isoCellSize
968 }
969 }
970 }
971
972 start := time.Now()
973
974 for n := runtime.NumCPU(); n >= 1; n-- {
975 wg.Add(1)
976 go rasterRow()
977 }
978
979 for i := 0; i < ycells; i++ {
980 rows <- i
981 }
982 close(rows)
983
984 wg.Wait()
985
986 feedback.Info("Rasterizing took %v.", time.Since(start))
987
988 feedback.Info("Trace iso areas.")
989
990 start = time.Now()
991
992 tracer := contourmap.FromFloat64s(xcells+2, ycells+2, raster)
993
994 areas := make([]wkb.MultiPolygonGeom, len(heights))
995
996 // TODO: Check if this correct!
997 reprojX := common.Linear(1, min.X, float64(xcells+1), max.X)
998 reprojY := common.Linear(1, min.Y, float64(ycells+1), max.Y)
999
1000 cnts := make(chan int)
1001
1002 doContours := func() {
1003 defer wg.Done()
1004 for hIdx := range cnts {
1005 c := tracer.Contours(heights[hIdx])
1006
1007 if len(c) == 0 {
1008 continue
1009 }
1010
1011 // We need to bring it back to the
1012 // none raster coordinate system.
1013 a := make(wkb.MultiPolygonGeom, len(c))
1014
1015 for i, pl := range c {
1016 shell := make(wkb.LinearRingGeom, len(pl))
1017 for j, pt := range pl {
1018 shell[j] = wkb.PointGeom{
1019 X: reprojX(pt.X),
1020 Y: reprojY(pt.Y),
1021 }
1022 }
1023 a[i] = wkb.PolygonGeom{shell}
1024 }
1025
1026 areas[hIdx] = a
1027 }
1028 }
1029
1030 for n := runtime.NumCPU(); n >= 1; n-- {
1031 wg.Add(1)
1032 go doContours()
1033 }
1034
1035 for i := range heights {
1036 cnts <- i
1037 }
1038 close(cnts)
1039
1040 wg.Wait()
1041
1042 feedback.Info("Tracing areas took %v.", time.Since(start))
1043
1044 // Raster and tracer are not needed any more.
1045 raster, tracer = nil, nil
1046 926
1047 return storeAreas( 927 return storeAreas(
1048 ctx, tx, feedback, 928 ctx, tx, feedback,
1049 areas, tree.EPSG, heights, id) 929 areas, tree.EPSG, heights, id)
1050 } 930 }