Mercurial > gemma
view pkg/imports/sr.go @ 5741:09d27eff745a default tip
Bumped version for development.
author | Sascha Wilde <wilde@sha-bang.de> |
---|---|
date | Mon, 09 Sep 2024 18:47:11 +0200 |
parents | 6270951dda28 |
children |
line wrap: on
line source
// 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> // * Tom Gottfried <tom@intevation.de> package imports import ( "archive/zip" "bufio" "context" "crypto/sha1" "database/sql" "encoding/hex" "errors" "fmt" "io" "math" "os" "path" "path/filepath" "regexp" "strconv" "strings" "sync" "time" shp "github.com/jonas-p/go-shp" "gemma.intevation.de/gemma/pkg/common" "gemma.intevation.de/gemma/pkg/log" "gemma.intevation.de/gemma/pkg/mesh" "gemma.intevation.de/gemma/pkg/models" "gemma.intevation.de/gemma/pkg/wkb" ) // SoundingResult is a Job to import sounding reults // from a ZIP file into the database. type SoundingResult struct { // Dir is the folder in the file system where the // 'sr.zip' is expect to be. Dir string `json:"dir"` // Override data // Date if given overrides the date value from the meta.json. Date *models.Date `json:"date,omitempty"` // Date if given overrides the name of the bottleneck from the meta.json. Bottleneck *string `json:"bottleneck,omitempty"` // EPSG if given overrides the EPSG code from the meta.json. // Defaults to WGS84. EPSG *uint `json:"epsg,omitempty"` // DepthReference if given overides the DepthReference value // from the meta.json. DepthReference *string `json:"depth-reference,omitempty"` // SingleBeam is kept in for compat. SingleBeam *bool `json:"single-beam,omitempty"` // SurveyType indicates that the sounding is a single beam scan. SurveyType *models.SurveyType `json:"survey-type,omitempty"` // NegateZ indicated that the Z values of thy XYZ input should be // multiplied by -1. NegateZ *bool `json:"negate-z,omitempty"` } const ( contourStepWidth = 0.1 contourTolerance = 0.1 ) const ( // multiBeamThreshold is the number of points m² which // is assumed that greater values are indicators for // an already interpolated point cloud. multiBeamThreshold = 1.0 / 5.0 ) const ( // pointsPerSquareMeter is the average number of points // when generating a artifical height model for single beam scans. pointsPerSquareMeter = 1.0 ) 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" type srJobCreator struct{} func init() { RegisterJobCreator(SRJobKind, srJobCreator{}) } func (srJobCreator) Description() string { return "sounding results" } func (srJobCreator) AutoAccept() bool { return false } func (srJobCreator) Create() Job { return new(SoundingResult) } // LoadingLogs ensures that log lines are loaded when import is exported. func (srJobCreator) LoadingLogs() bool { return true } func (srJobCreator) Depends() [2][]string { return [2][]string{ {"sounding_results", "sounding_results_iso_areas", "sounding_results_marking_points"}, {"bottlenecks"}, } } // StageDone moves the imported sounding result out of the staging area. func (srJobCreator) StageDone( ctx context.Context, tx *sql.Tx, id int64, _ Feedback, ) error { _, err := tx.ExecContext(ctx, srStageDoneSQL, id) return err } const ( srStageDoneSQL = ` UPDATE waterway.sounding_results SET staging_done = true WHERE id = ( SELECT key from import.track_imports WHERE import_id = $1 AND relation = 'waterway.sounding_results'::regclass)` insertHullSQL = ` INSERT INTO waterway.sounding_results ( bottleneck_id, date_info, depth_reference, area, surtyp, zpg_exception ) VALUES ( $1, $2::date, $3, (SELECT CASE WHEN $5::bytea IS NULL THEN ST_Transform(ST_ConcaveHull(ST_Force2D(ST_GeomFromWKB($4, $6::integer)), 0.7), 4326)::geography ELSE ST_MakeValid(ST_Transform(ST_GeomFromWKB($5, $6::integer), 4326))::geography END), $7, $8 ) RETURNING id, ST_X(ST_Centroid(area::geometry)), ST_Y(ST_Centroid(area::geometry)), best_utm(area), ST_AsBinary(ST_Transform(area::geometry, best_utm(area))) ` reprojectPointsBufferedSQL = ` SELECT ST_AsBinary(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer)), ST_AsBinary(ST_Buffer(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer), 0.1)), ST_Area(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer))` insertMeshSQL = ` UPDATE waterway.sounding_results SET mesh_checksum = $2, mesh_index = $3, mesh_index_version = $4 WHERE id = $1` repairBoundarySQL = ` SELECT ST_AsBinary(ST_Buffer(ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 0.0)), ST_AsBinary(ST_Buffer(ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 0.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_GeomFromWKB($4, $3::integer), $5), 3 ) ), 4326 ) FROM waterway.sounding_results sr WHERE id = $1 ` insertMarkingPointsSQL = ` INSERT INTO waterway.sounding_results_marking_points ( sounding_result_id, height, points ) SELECT $1, $2, ST_Transform(ST_GeomFromWKB($4, $3::integer), 4326) FROM waterway.sounding_results sr WHERE id = $1 ` selectGaugeLDCSQL = ` SELECT grwl.value, grwl.depth_reference FROM waterway.gauges_reference_water_levels grwl JOIN waterway.bottlenecks bns ON grwl.location = bns.gauge_location AND grwl.validity @> CAST($2 AS timestamptz) WHERE bns.bottleneck_id = $1 AND bns.validity @> CAST($2 AS timestamptz) AND grwl.depth_reference like 'LDC%' ` selectZPGExceptionAllowedSQL = ` SELECT rolname = 'sys_admin' OR country IN ('BG', 'RO') FROM users.list_users WHERE username = current_user` reprojectPointsSingleBeamSQL = ` 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)) ` ) var ( bottleneckRe, surveyTypeRe, dateRe, negateZRe *regexp.Regexp recoverRegsOnce sync.Once ) func compileRecoverRegs() { bottleneckRe = regexp.MustCompile(`Bottleneck:\s*(.+)\s*$`) surveyTypeRe = regexp.MustCompile(`Processing as\s+([^\s]+)\s+beam scan.`) dateRe = regexp.MustCompile(`Survey date:\s*(\d{4}-\d{2}-\d{2})`) negateZRe = regexp.MustCompile(`Z values will be negated\.`) } // Description gives a short info about relevant facts of this import. func (sr *SoundingResult) Description(msgs []string) (string, error) { recoverRegsOnce.Do(compileRecoverRegs) //log.Debugln(strings.Join(msgs, "\n") + "\n\n") var ( bottleneck, st, date string negZ bool ) for _, msg := range msgs { if m := bottleneckRe.FindStringSubmatch(msg); m != nil { bottleneck = m[1] continue } if m := surveyTypeRe.FindStringSubmatch(msg); m != nil { st = m[1] continue } if m := dateRe.FindStringSubmatch(msg); m != nil { date = m[1] continue } if negateZRe.MatchString(msg) { negZ = true } } var descs []string if sr.Bottleneck != nil { descs = append(descs, *sr.Bottleneck) } else if bottleneck != "" { log.Debugf("bottleneck recovered: %s\n", bottleneck) descs = append(descs, bottleneck) } if sr.Date != nil { descs = append(descs, (*sr).Date.Format(common.DateFormat)) } else if date != "" { log.Debugf("date recovered: %s\n", date) descs = append(descs, date) } if sr.NegateZ != nil && *sr.NegateZ { descs = append(descs, "negateZ") } else if negZ { log.Debugln("negateZ recovered") descs = append(descs, "negateZ") } if sr.SurveyType != nil { descs = append(descs, string(*sr.SurveyType)) } else if st != "" { log.Debugf("survey type recovered: %s\n", st) descs = append(descs, st) } return strings.Join(descs, "|"), nil } func (sr *SoundingResult) negateZ() bool { return sr.NegateZ != nil && *sr.NegateZ } // Do executes the actual sounding result import. func (sr *SoundingResult) Do( ctx context.Context, importID int64, conn *sql.Conn, feedback Feedback, ) (any, error) { start := time.Now() zpath := filepath.Join(sr.Dir, "sr.zip") z, err := zip.OpenReader(zpath) if err != nil { feedback.Info("%v", err) feedback.Info("Falling back to TXT file mode.") z = nil } if z != nil { 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'") } m, err := sr.loadMeta(mf) if err != nil { return nil, err } feedback.Info("Bottleneck: %s", m.Bottleneck) feedback.Info("Survey date: %s", m.Date.Format(common.DateFormat)) var xform vertexTransform if sr.negateZ() { feedback.Info("Z values will be negated.") xform = negateZTransform } else { xform = identityTransform } var zpgException bool if m.DepthReference == "ZPG" { feedback.Info("Found ZPG as reference system -> translating Z values to LDC") var ldc float64 var depthReference string err := conn.QueryRowContext(ctx, selectGaugeLDCSQL, m.Bottleneck, m.Date.Time, ).Scan( &ldc, &depthReference, ) switch { case err == sql.ErrNoRows: if err := conn.QueryRowContext( ctx, selectZPGExceptionAllowedSQL).Scan(&zpgException); err != nil { return nil, err } if !zpgException { return nil, errors.New("cannot load LDC value") } feedback.Warn("Not LCD found, but ZPG exception granted") case err != nil: return nil, err } if !zpgException { // LDC is cm. The data is in m. ldc /= 100 xform = chainTransforms( xform, func(v mesh.Vertex) mesh.Vertex { return mesh.Vertex{X: v.X, Y: v.Y, Z: v.Z + ldc} }) m.DepthReference = depthReference } } if err := m.Validate(ctx, conn); err != nil { return nil, common.ToError(err) } var xyz mesh.MultiPointZ if z != nil { // Scanning ZIP file for *.xyz file. var xyzf *zip.File for _, ext := range []string{".xyz", ".txt"} { feedback.Info("Looking for '*%s'", ext) if xyzf = common.FindInZIP(z, ext); xyzf != nil { break } } if xyzf == nil { return nil, errors.New("cannot find any *.xyz or *.txt file") } xyz, err = loadXYZ(xyzf, feedback, xform) } else { // TXT file mode xyz, err = loadXYZFile(zpath, feedback, xform) } if err != nil { return nil, err } if len(xyz) == 0 { return nil, errors.New("XYZ does not contain any vertices") } // Is there a boundary shapefile in the ZIP archive? boundary, err := loadBoundary(z) if err != nil { return nil, err } tx, err := conn.BeginTx(ctx, nil) if err != nil { return nil, err } defer tx.Rollback() summary, err := sr.processScan( ctx, tx, feedback, importID, m, xyz, boundary, zpgException, ) if err != nil { return nil, err } if err = tx.Commit(); err != nil { feedback.Error( "Storing sounding result failed after %v.", time.Since(start)) return nil, err } feedback.Info( "Storing sounding result was successful after %v.", time.Since(start)) return summary, nil } func (sr *SoundingResult) processScan( ctx context.Context, tx *sql.Tx, feedback Feedback, importID int64, m *models.SoundingResultMeta, xyz mesh.MultiPointZ, boundary polygonSlice, zpgException bool, ) (any, error) { feedback.Info("Processing as %s beam scan.", m.SurveyType) feedback.Info("Reproject XYZ data.") start := time.Now() var reproj []byte var epsg uint32 if err := tx.QueryRowContext( ctx, reprojectPointsSingleBeamSQL, xyz.AsWKB(), m.EPSG, ).Scan(&reproj, &epsg); err != nil { 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)) start = time.Now() xyz = xyz.QuantizeXY(mesh.QuantScale) feedback.Info("Quantizing points in X/Y to 1/%d took %v.", mesh.QuantScale, time.Since(start)) feedback.Info("Number of reprojected points: %d", len(xyz)) feedback.Info("Triangulate XYZ data.") start = time.Now() tri, err := mesh.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) var ( clippingPolygon mesh.Polygon clippingPolygonBuffered mesh.Polygon removed map[int32]struct{} polygonArea float64 clippingPolygonWKB []byte ) if boundary == nil { feedback.Info("No boundary given. Calulate from XYZ data.") tooLongEdge := tri.EstimateTooLong() feedback.Info("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge) var polygon mesh.LineStringZ start = time.Now() polygon, removed = tri.ConcaveHull(tooLongEdge) polygonArea = polygon.Area() if polygonArea < 0.0 { // counter clockwise polygonArea = -polygonArea polygon.Reverse() } clippingPolygon.FromLineStringZ(polygon) var buffered []byte if err := tx.QueryRowContext( ctx, repairBoundarySQL, clippingPolygon.AsWKB(), epsg, ).Scan(&clippingPolygonWKB, &buffered); err != nil { return nil, err } if err := clippingPolygon.FromWKB(clippingPolygonWKB); err != nil { return nil, err } if err := clippingPolygonBuffered.FromWKB(buffered); err != nil { return nil, err } feedback.Info("Calculating took %v.", time.Since(start)) } else { // has Boundary feedback.Info("Using uploaded boundary polygon.") var buffered []byte if err = tx.QueryRowContext( ctx, reprojectPointsBufferedSQL, boundary.asWKB(), m.EPSG, epsg, ).Scan(&clippingPolygonWKB, &buffered, &polygonArea); err != nil { return nil, err } if err := clippingPolygon.FromWKB(clippingPolygonWKB); err != nil { return nil, err } if err := clippingPolygonBuffered.FromWKB(buffered); err != nil { return nil, err } tin := tri.Tin() tin.EPSG = epsg var str mesh.STRTree str.Build(tin) removed = str.Clip(&clippingPolygon) } if m.SurveyType == models.SurveyTypeSingleBeam { origDensity := float64(len(xyz)) / polygonArea feedback.Info("Boundary area: %.2fm²", polygonArea) feedback.Info("Original point density: %.2f points/m²", origDensity) if origDensity > multiBeamThreshold { feedback.Warn("The density is greater than %.2f points/m².", multiBeamThreshold) feedback.Warn("It is assumed that the data is already interpolated.") goto multibeam } // Build the first mesh to generate random points on. feedback.Info("Build virtual DEM based on original XYZ data.") start = time.Now() tin := tri.Tin() var virtual mesh.STRTree virtual.BuildWithout(tin, removed) feedback.Info("Building took %v", time.Since(start)) numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter)) feedback.Info("Generate %d random points for an average density of ~%.2f points/m².", numPoints, pointsPerSquareMeter) start = time.Now() generated := make(mesh.LineStringZ, 0, numPoints+clippingPolygon.NumVertices(0)) mesh.GenerateRandomVertices( numPoints, tin.Min, tin.Max, virtual.Value, func(vertices []mesh.Vertex) { generated = append(generated, vertices...) }) feedback.Info("Generating %d points took %v.", len(generated), time.Since(start)) // Add the boundary to new point cloud. dupes := map[[2]float64]struct{}{} clippingPolygon.Vertices(0, func(x, y float64) { key := [2]float64{x, y} if _, found := dupes[key]; found { return } dupes[key] = struct{}{} if z, ok := virtual.Value(x, y); ok { generated = append(generated, mesh.Vertex{X: x, Y: y, Z: z}) } }) feedback.Info("Triangulate new point cloud.") xyz = mesh.MultiPointZ(generated) start = time.Now() tri, err = mesh.Triangulate(xyz) if err != nil { return nil, err } feedback.Info("Second triangulation took %v.", time.Since(start)) feedback.Info("Number triangles: %d.", len(tri.Triangles)/3) feedback.Info("Clipping triangles from new mesh.") } multibeam: final := mesh.STRTree{Entries: 16} if m.SurveyType != models.SurveyTypeMarking { start = time.Now() tin := tri.Tin() tin.EPSG = epsg var str mesh.STRTree str.Build(tin) feedback.Info("Building clipping index took %v", time.Since(start)) start = time.Now() clippingPolygonBuffered.Indexify() removed = str.Clip(&clippingPolygonBuffered) feedback.Info("Clipping took %v.", time.Since(start)) feedback.Info("Number of triangles to clip %d.", len(removed)) start = time.Now() final.BuildWithout(tin, removed) feedback.Info("Building final mesh took %v.", time.Since(start)) feedback.Info("Store mesh.") } else { start = time.Now() clippingPolygonBuffered.Indexify() before := len(xyz) xyz = xyz.Filter(func(v mesh.Vertex) bool { return clippingPolygonBuffered.IntersectionBox2D(v.Box2D()) != mesh.IntersectionOutSide }) feedback.Info("Clipping took %v.", time.Since(start)) feedback.Info("Number of points to clip: %d.", before-len(xyz)) } start = time.Now() var ( id int64 dummy uint32 lat, lon float64 hull []byte ) switch err := tx.QueryRowContext( ctx, insertHullSQL, m.Bottleneck, m.Date.Time, m.DepthReference, nil, clippingPolygonWKB, epsg, m.SurveyType, zpgException, ).Scan( &id, &lat, &lon, &dummy, &hull, ); { case err == sql.ErrNoRows: return nil, fmt.Errorf( "no matching bottleneck of given name or time available: %v", err) case err != nil: return nil, err } if m.SurveyType != models.SurveyTypeMarking { final.OptimizeForSerialization(-1) var index []byte var version int index, version, err = final.Bytes(-1) if err != nil { return nil, err } h := sha1.New() h.Write(index) checksum := hex.EncodeToString(h.Sum(nil)) _, err = tx.ExecContext(ctx, insertMeshSQL, id, checksum, index, version) if err != nil { return nil, err } feedback.Info("Storing mesh index took %s.", time.Since(start)) err = generateIsoAreas( ctx, tx, feedback, &final, loadClassBreaks(ctx, tx, feedback, final.Min().Z, final.Max().Z), id) } else { // SurveyTypeMarking err = generateMarkingPoints( ctx, tx, feedback, xyz, epsg, id) } if err != nil { return nil, err } // Store for potential later removal. if err = track(ctx, tx, importID, "waterway.sounding_results", id); err != nil { return nil, err } summary := struct { Bottleneck string `json:"bottleneck"` Date models.Date `json:"date"` Lat float64 `json:"lat"` Lon float64 `json:"lon"` }{ Bottleneck: m.Bottleneck, Date: m.Date, Lat: lat, Lon: lon, } return &summary, nil } // CleanUp removes the folder containing the ZIP file with the // the sounding result import. func (sr *SoundingResult) CleanUp() error { return os.RemoveAll(sr.Dir) } func (sr *SoundingResult) completeOverride() bool { // sr.EPSG == nil -> WGS84 return sr.Bottleneck != nil && sr.Date != nil && sr.DepthReference != nil && (sr.SingleBeam != nil || sr.SurveyType != nil) && sr.NegateZ != nil } func (sr *SoundingResult) loadMeta(f *zip.File) (*models.SoundingResultMeta, error) { if f == nil { var epsg uint if sr.EPSG != nil { epsg = *sr.EPSG } else { epsg = models.WGS84 } st := models.SurveyTypeMultiBeam if sr.SingleBeam != nil && *sr.SingleBeam { st = models.SurveyTypeSingleBeam } if sr.SurveyType != nil { st = *sr.SurveyType } return &models.SoundingResultMeta{ Date: *sr.Date, Bottleneck: *sr.Bottleneck, EPSG: epsg, DepthReference: *sr.DepthReference, SurveyType: st, NegateZ: sr.negateZ(), }, nil } r, err := f.Open() if err != nil { return nil, err } defer r.Close() var m models.SoundingResultMeta if err := m.Decode(r); err != nil { return nil, err } // Apply overrides if sr.Date != nil { m.Date = *sr.Date } if sr.Bottleneck != nil { m.Bottleneck = *sr.Bottleneck } if sr.EPSG != nil { m.EPSG = *sr.EPSG } if sr.DepthReference != nil { m.DepthReference = *sr.DepthReference } // Kept in for compat if sr.SingleBeam != nil { if *sr.SingleBeam { m.SurveyType = models.SurveyTypeSingleBeam } else { m.SurveyType = models.SurveyTypeMultiBeam } } if sr.SurveyType != nil { m.SurveyType = *sr.SurveyType } if sr.NegateZ != nil { m.NegateZ = *sr.NegateZ } return &m, nil } type vertexTransform func(mesh.Vertex) mesh.Vertex func identityTransform(v mesh.Vertex) mesh.Vertex { return v } func negateZTransform(v mesh.Vertex) mesh.Vertex { return mesh.Vertex{X: v.X, Y: v.Y, Z: -v.Z} } func chainTransforms(a, b vertexTransform) vertexTransform { return func(v mesh.Vertex) mesh.Vertex { return b(a(v)) } } func loadXYZReader(r io.Reader, feedback Feedback, xform vertexTransform) (mesh.MultiPointZ, error) { mpz := make(mesh.MultiPointZ, 0, 250000) s := bufio.NewScanner(r) warnLimiter := common.WarningLimiter{Log: feedback.Warn, MaxWarnings: 100} warn := warnLimiter.Warn defer warnLimiter.Close() for line := 1; s.Scan(); line++ { text := s.Text() var p mesh.Vertex // fmt.Sscanf(text, "%f,%f,%f") is 4 times slower. idx := strings.IndexByte(text, ',') if idx == -1 { warn("format error in line %d", line) continue } var err error if p.X, err = strconv.ParseFloat(text[:idx], 64); err != nil { warn("format error in line %d: %v", line, err) continue } text = text[idx+1:] if idx = strings.IndexByte(text, ','); idx == -1 { feedback.Warn("format error in line %d", line) continue } if p.Y, err = strconv.ParseFloat(text[:idx], 64); err != nil { warn("format error in line %d: %v", line, err) continue } text = text[idx+1:] if p.Z, err = strconv.ParseFloat(text, 64); err != nil { warn("format error in line %d: %v", line, err) continue } p = xform(p) mpz = append(mpz, p) } if err := s.Err(); err != nil { return nil, err } return mpz, nil } func loadXYZ(f *zip.File, feedback Feedback, xform vertexTransform) (mesh.MultiPointZ, error) { r, err := f.Open() if err != nil { return nil, err } defer r.Close() return loadXYZReader(r, feedback, xform) } func loadXYZFile(f string, feedback Feedback, xform vertexTransform) (mesh.MultiPointZ, error) { r, err := os.Open(f) if err != nil { return nil, err } defer r.Close() return loadXYZReader(r, feedback, xform) } func loadBoundary(z *zip.ReadCloser) (polygonSlice, error) { shpF := common.FindInZIP(z, ".shp") if shpF == nil { return nil, nil } prefix := strings.TrimSuffix(shpF.Name, path.Ext(shpF.Name)) dbfF := common.FindInZIP(z, prefix+".dbf") if dbfF == nil { return nil, fmt.Errorf("no DBF file found for %s", shpF.Name) } shpR, err := shpF.Open() if err != nil { return nil, err } dbfR, err := dbfF.Open() if err != nil { shpR.Close() return nil, err } sr := shp.SequentialReaderFromExt(shpR, dbfR) defer sr.Close() if !sr.Next() { return nil, sr.Err() } _, s := sr.Shape() if s == nil { return nil, sr.Err() } return shapeToPolygon(s) } func defaultClassBreaks(min, max float64) mesh.ClassBreaks { var heights mesh.ClassBreaks h := contourStepWidth * math.Ceil(min/contourStepWidth) for ; h <= max; h += contourStepWidth { heights = append(heights, h) } return heights } func generateMarkingPoints( ctx context.Context, tx *sql.Tx, feedback Feedback, xyz mesh.MultiPointZ, epsg uint32, id int64, ) error { heights, err := mesh.LoadClassBreaks( ctx, tx, "morphology_classbreaks") if err != nil { feedback.Warn("Loading class breaks failed: %v", err) feedback.Info("Using default class breaks") min, max := xyz.MinMax() heights = defaultClassBreaks(min.Z, max.Z) } heights = heights.Dedup() classes := heights.Classify(xyz) stmt, err := tx.PrepareContext(ctx, insertMarkingPointsSQL) if err != nil { return err } defer stmt.Close() for i, class := range classes { // Ignore empty classes if len(class) == 0 { continue } if _, err := stmt.ExecContext( ctx, id, heights[i], epsg, class.AsWKB(), ); err != nil { return err } } return nil } func loadClassBreaks( ctx context.Context, tx *sql.Tx, feedback Feedback, minZ, maxZ float64, ) mesh.ClassBreaks { heights, err := mesh.LoadClassBreaks( ctx, tx, "morphology_classbreaks") if err != nil { feedback.Warn("Loading class breaks failed: %v", err) feedback.Info("Using default class breaks") heights = defaultClassBreaks(minZ, maxZ) } else { heights = heights.ExtrapolateClassBreaks(minZ, maxZ) } return heights.Dedup() } func generateIsoAreas( ctx context.Context, tx *sql.Tx, feedback Feedback, tree *mesh.STRTree, heights mesh.ClassBreaks, id int64, ) error { feedback.Info("Generate iso areas") total := time.Now() defer func() { feedback.Info("Generating iso areas took %s.", time.Since(total)) }() box := mesh.Box2D{ X1: tree.Min().X, Y1: tree.Min().Y, X2: tree.Max().X, Y2: tree.Max().Y, } raster := mesh.NewRaster(box, isoCellSize) raster.Rasterize(tree.Value) areas := raster.Trace(heights) return storeAreas( ctx, tx, feedback, areas, tree.EPSG(), heights, id) } func storeAreas( ctx context.Context, tx *sql.Tx, feedback Feedback, areas []wkb.MultiPolygonGeom, epsg uint32, heights mesh.ClassBreaks, id int64, ) error { feedback.Info("Store 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() var size int for i, a := range areas { if len(a) == 0 { continue } wkb := a.AsWKB() size += len(wkb) if _, err := stmt.ExecContext( ctx, id, heights[i], epsg, wkb, contourTolerance, ); err != nil { return err } } feedback.Info("Transferred WKB size: %.2fMB.", float64(size)/(1024*1024)) return nil }