diff --git a/sieve.go b/sieve.go index eab6652..90f1001 100644 --- a/sieve.go +++ b/sieve.go @@ -61,7 +61,7 @@ func geometryTypeFromString(geometrytype string) gpkg.GeometryType { // createSQL creates a CREATE statement on the given table and column information // used for creating feature tables in the target Geopackage func (t table) createSQL() string { - create := fmt.Sprintf(`CREATE TABLE IF NOT EXISTS %v`, t.name) + create := fmt.Sprintf(`CREATE TABLE IF NOT EXISTS "%v"`, t.name) var columnparts []string for _, column := range t.columns { columnpart := column.name + ` ` + column.ctype @@ -86,7 +86,7 @@ func (t table) selectSQL() string { for _, c := range t.columns { csql = append(csql, c.name) } - query := `SELECT ` + strings.Join(csql, `,`) + ` FROM ` + t.name + `;` + query := `SELECT ` + strings.Join(csql, `,`) + ` FROM "` + t.name + `";` return query } @@ -102,7 +102,7 @@ func (t table) insertSQL() string { } csql = append(csql, t.gcolumn) vsql = append(vsql, `?`) - query := `INSERT INTO ` + t.name + `(` + strings.Join(csql, `,`) + `) VALUES(` + strings.Join(vsql, `,`) + `)` + query := `INSERT INTO "` + t.name + `"(` + strings.Join(csql, `,`) + `) VALUES(` + strings.Join(vsql, `,`) + `)` return query } @@ -110,9 +110,8 @@ func (t table) insertSQL() string { func getSourceTableInfo(h *gpkg.Handle) []table { query := `SELECT table_name, column_name, geometry_type_name, srs_id FROM gpkg_geometry_columns;` rows, err := h.Query(query) - defer rows.Close() if err != nil { - log.Printf("err during closing rows: %v - %v", query, err) + log.Fatalf("error during closing rows: %v - %v", query, err) } var tables []table @@ -122,7 +121,7 @@ func getSourceTableInfo(h *gpkg.Handle) []table { var srsID int err := rows.Scan(&t.name, &t.gcolumn, >ype, &srsID) if err != nil { - log.Fatal(err) + log.Fatalf("error ready the source table information: %s", err) } t.columns = getTableColumns(h, t.name) @@ -131,6 +130,7 @@ func getSourceTableInfo(h *gpkg.Handle) []table { tables = append(tables, t) } + defer rows.Close() return tables } @@ -138,24 +138,14 @@ func getSourceTableInfo(h *gpkg.Handle) []table { func getSpatialReferenceSystem(h *gpkg.Handle, id int) gpkg.SpatialReferenceSystem { var srs gpkg.SpatialReferenceSystem query := `SELECT srs_name, srs_id, organization, organization_coordsys_id, definition, description FROM gpkg_spatial_ref_sys WHERE srs_id = %v;` - rows, err := h.Query(fmt.Sprintf(query, id)) - defer rows.Close() - if err != nil { - log.Printf("err during closing rows: %v - %v", query, err) - } - for rows.Next() { - var description *string - err := rows.Scan(&srs.Name, &srs.ID, &srs.Organization, &srs.OrganizationCoordsysID, &srs.Definition, &description) - if description != nil { - srs.Description = *description - } - if err != nil { - log.Fatal(err) - } - // On first hit (and only) return - return srs + row := h.QueryRow(fmt.Sprintf(query, id)) + var description *string + row.Scan(&srs.Name, &srs.ID, &srs.Organization, &srs.OrganizationCoordsysID, &srs.Definition, &description) + if description != nil { + srs.Description = *description } + return srs } @@ -164,19 +154,20 @@ func getTableColumns(h *gpkg.Handle, table string) []column { var columns []column query := `PRAGMA table_info('%v');` rows, err := h.Query(fmt.Sprintf(query, table)) - defer rows.Close() + if err != nil { - log.Printf("err during closing rows: %v - %v", query, err) + log.Fatalf("err during closing rows: %v - %v", query, err) } for rows.Next() { var column column err := rows.Scan(&column.cid, &column.name, &column.ctype, &column.notnull, &column.dfltValue, &column.pk) if err != nil { - log.Fatal(err) + log.Fatalf("error getting the column information: %s", err) } columns = append(columns, column) } + defer rows.Close() return columns } @@ -185,8 +176,7 @@ func buildTable(h *gpkg.Handle, t table) error { query := t.createSQL() _, err := h.Exec(query) if err != nil { - log.Println("err:", err) - return err + log.Fatalf("error building table in target GeoPackage: %s", err) } err = h.AddGeometryTable(gpkg.TableDescription{ @@ -201,7 +191,7 @@ func buildTable(h *gpkg.Handle, t table) error { M: gpkg.Prohibited, }) if err != nil { - log.Println("err:", err) + log.Println("error adding geometry table in target GeoPackage:", err) return err } return nil @@ -227,14 +217,13 @@ func initTargetGeopackage(h *gpkg.Handle, tables []table) error { // and decodes the WKB geometry to a geom.Polygon func readFeatures(h *gpkg.Handle, preSieve chan feature, t table) { rows, err := h.Query(t.selectSQL()) - defer rows.Close() if err != nil { - log.Printf("err during closing rows: %v", err) + log.Fatalf("err during closing rows: %s", err) } cols, err := rows.Columns() if err != nil { - log.Println("err:", err) + log.Fatalf("error reading the columns: %s", err) } for rows.Next() { @@ -245,8 +234,7 @@ func readFeatures(h *gpkg.Handle, preSieve chan feature, t table) { } if err = rows.Scan(valPtrs...); err != nil { - log.Printf("err reading row values: %v", err) - return + log.Fatalf("err reading row values: %v", err) } var f feature var c []interface{} @@ -256,7 +244,7 @@ func readFeatures(h *gpkg.Handle, preSieve chan feature, t table) { case t.gcolumn: wkbgeom, err := gpkg.DecodeGeometry(vals[i].([]byte)) if err != nil { - log.Fatal(err) + log.Fatalf("error decoding the geometry: %s", err) } f.geometry = wkbgeom.Geometry default: @@ -278,7 +266,7 @@ func readFeatures(h *gpkg.Handle, preSieve chan feature, t table) { case nil: c = append(c, v) default: - log.Printf("unexpected type for sqlite column data: %v: %T", cols[i], v) + log.Fatalf("unexpected type for sqlite column data: %v: %T", cols[i], v) } } f.columns = c @@ -290,6 +278,7 @@ func readFeatures(h *gpkg.Handle, preSieve chan feature, t table) { log.Fatal(err) } close(preSieve) + defer rows.Close() } // sieveFeatures sieves/filters the geometry against the given resolution @@ -325,61 +314,32 @@ func sieveFeatures(preSieve chan feature, postSieve chan feature, resolution flo close(postSieve) } -// multiPolygonSieve will split it self into the separated polygons that will be sieved before building a new MULTIPOLYGON -func multiPolygonSieve(mp geom.MultiPolygon, resolution float64) geom.MultiPolygon { - var sievedMultiPolygon geom.MultiPolygon - for _, p := range mp { - if sievedPolygon := polygonSieve(p, resolution); sievedPolygon != nil { - sievedMultiPolygon = append(sievedMultiPolygon) - } - } - return sievedMultiPolygon -} - -// polygonSieve will sieve a given POLYGON -func polygonSieve(p geom.Polygon, resolution float64) geom.Polygon { - minArea := resolution * resolution - if area(p) > minArea { - if len(p) > 1 { - var sievedPolygon geom.Polygon - sievedPolygon = append(sievedPolygon, p[0]) - for _, interior := range p[1:] { - if shoelace(interior) > minArea { - sievedPolygon = append(sievedPolygon, interior) - } - } - return sievedPolygon - } - return p - } - return nil -} - // writeFeatures writes the features processed by the sieveFeatures to the geopackages -func writeFeatures(postSieve chan feature, kill chan bool, h *gpkg.Handle, t table) { +func writeFeatures(postSieve chan feature, kill chan bool, h *gpkg.Handle, t table, p int) { var ext *geom.Extent - stmt, err := h.Prepare(t.insertSQL()) - if err != nil { - log.Println("err:", err) - return - } + var err error + + var features [][]interface{} for { feature, hasMore := <-postSieve if !hasMore { + writeFeaturesArray(features, h, t) + features = nil break } else { sb, err := gpkg.NewBinary(int32(t.srs.ID), feature.geometry) if err != nil { - log.Fatalln("err:", err) + log.Fatalf("Could not create a binary geometry: %s", err) } data := feature.columns data = append(data, sb) + features = append(features, data) - _, err = stmt.Exec(data...) - if err != nil { - log.Fatalln("stmt err:", err) + if len(features)%p == 0 { + writeFeaturesArray(features, h, t) + features = nil } } @@ -387,7 +347,7 @@ func writeFeatures(postSieve chan feature, kill chan bool, h *gpkg.Handle, t tab ext, err = geom.NewExtentFromGeometry(feature.geometry) if err != nil { ext = nil - log.Println("err:", err) + log.Println("Failed to create new extent:", err) continue } } else { @@ -398,24 +358,45 @@ func writeFeatures(postSieve chan feature, kill chan bool, h *gpkg.Handle, t tab kill <- true } +func writeFeaturesArray(features [][]interface{}, h *gpkg.Handle, t table) { + tx, err := h.Begin() + if err != nil { + log.Fatalf("Could not start a transaction: %s", err) + } + + stmt, err := tx.Prepare(t.insertSQL()) + if err != nil { + log.Fatalf("Could not prepare a statement: %s", err) + } + + for _, f := range features { + _, err = stmt.Exec(f...) + if err != nil { + log.Fatalf("Could a result summary from the prepared statement: %s", err) + } + } + + stmt.Close() + tx.Commit() +} + func main() { log.Println("start") sourceGeopackage := flag.String("s", "empty", "source geopackage") targetGeopackage := flag.String("t", "empty", "target geopackage") + pageSize := flag.Int("p", 1000, "target geopackage") resolution := flag.Float64("r", 0.0, "resolution for sieving") flag.Parse() srcHandle, err := gpkg.Open(*sourceGeopackage) if err != nil { - log.Println("err:", err) - return + log.Fatalf("error opening source GeoPackage: %s", err) } defer srcHandle.Close() trgHandle, err := gpkg.Open(*targetGeopackage) if err != nil { - log.Println("Open err:", err) - return + log.Fatalf("error opening target GeoPackage: %s", err) } defer trgHandle.Close() @@ -423,8 +404,7 @@ func main() { err = initTargetGeopackage(trgHandle, tables) if err != nil { - log.Println("err:", err) - return + log.Fatalf("error initialization the target GeoPackage: %s", err) } // Process the tables sequential @@ -434,7 +414,7 @@ func main() { postSieve := make(chan feature) kill := make(chan bool) - go writeFeatures(postSieve, kill, trgHandle, table) + go writeFeatures(postSieve, kill, trgHandle, table, *pageSize) go sieveFeatures(preSieve, postSieve, *resolution) go readFeatures(srcHandle, preSieve, table) @@ -450,9 +430,42 @@ func main() { log.Println("stop") } +// multiPolygonSieve will split it self into the separated polygons that will be sieved before building a new MULTIPOLYGON +func multiPolygonSieve(mp geom.MultiPolygon, resolution float64) geom.MultiPolygon { + var sievedMultiPolygon geom.MultiPolygon + for _, p := range mp { + if sievedPolygon := polygonSieve(p, resolution); sievedPolygon != nil { + sievedMultiPolygon = append(sievedMultiPolygon, sievedPolygon) + } + } + return sievedMultiPolygon +} + +// polygonSieve will sieve a given POLYGON +func polygonSieve(p geom.Polygon, resolution float64) geom.Polygon { + minArea := resolution * resolution + if area(p) > minArea { + if len(p) > 1 { + var sievedPolygon geom.Polygon + sievedPolygon = append(sievedPolygon, p[0]) + for _, interior := range p[1:] { + if shoelace(interior) > minArea { + sievedPolygon = append(sievedPolygon, interior) + } + } + return sievedPolygon + } + return p + } + return nil +} + // calculate the area of a polygon func area(geom [][][2]float64) float64 { interior := .0 + if geom == nil { + return 0. + } if len(geom) > 1 { for _, i := range geom[1:] { interior = interior + shoelace(i) @@ -464,6 +477,10 @@ func area(geom [][][2]float64) float64 { // https://en.wikipedia.org/wiki/Shoelace_formula func shoelace(pts [][2]float64) float64 { sum := 0. + if len(pts) == 0 { + return 0. + } + p0 := pts[len(pts)-1] for _, p1 := range pts { sum += p0[1]*p1[0] - p0[0]*p1[1] diff --git a/sieve_test.go b/sieve_test.go new file mode 100644 index 0000000..39266f6 --- /dev/null +++ b/sieve_test.go @@ -0,0 +1,121 @@ +package main + +import "testing" + +func TestShoelace(t *testing.T) { + var tests = []struct { + pts [][2]float64 + area float64 + }{ + // Rectangle + 0: {pts: [][2]float64{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}, area: float64(100)}, + // Triangle + 1: {pts: [][2]float64{{0, 0}, {5, 10}, {0, 10}, {0, 0}}, area: float64(25)}, + // Missing 'official closing point + 2: {pts: [][2]float64{{0, 0}, {0, 10}, {10, 10}, {10, 0}}, area: float64(100)}, + // Single point + 3: {pts: [][2]float64{{1234, 4321}}, area: float64(0.000000)}, + // No point + 4: {pts: nil, area: float64(0.000000)}, + // Empty point + 5: {pts: [][2]float64{}, area: float64(0.000000)}, + } + + for k, test := range tests { + area := shoelace(test.pts) + if area != test.area { + t.Errorf("test: %d, expected: %f \ngot: %f", k, test.area, area) + } + } +} + +func TestArea(t *testing.T) { + var tests = []struct { + geom [][][2]float64 + area float64 + }{ + // Rectangle + 0: {geom: [][][2]float64{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}}, area: float64(100)}, + // Rectangle with hole + 1: {geom: [][][2]float64{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}, {{2, 2}, {2, 8}, {8, 8}, {8, 2}, {2, 2}}}, area: float64(64)}, + // Rectangle with empty hole + 2: {geom: [][][2]float64{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}, {}}, area: float64(100)}, + // Rectangle with nil hole + 3: {geom: [][][2]float64{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}, nil}, area: float64(100)}, + // nil geometry + 4: {geom: nil, area: float64(0)}, + } + + for k, test := range tests { + area := area(test.geom) + if area != test.area { + t.Errorf("test: %d, expected: %f \ngot: %f", k, test.area, area) + } + } +} + +func TestPolygonSieve(t *testing.T) { + var tests = []struct { + geom [][][2]float64 + resolution float64 + sieved [][][2]float64 + }{ + // Lower resolution + 0: {geom: [][][2]float64{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}}, resolution: float64(9), sieved: [][][2]float64{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}}}, + // Higher resolution + 1: {geom: [][][2]float64{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}}, resolution: float64(101), sieved: nil}, + // Nil input + 2: {geom: nil, resolution: float64(1), sieved: nil}, + // Filterout donut + 3: {geom: [][][2]float64{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}, {{5, 5}, {5, 6}, {6, 6}, {6, 5}, {5, 5}}}, resolution: float64(9), sieved: [][][2]float64{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}}}, + // Donut stays + 4: {geom: [][][2]float64{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}, {{5, 5}, {5, 8.5}, {8.5, 8.5}, {8.5, 5}, {5, 5}}}, resolution: float64(3), sieved: [][][2]float64{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}, {{5, 5}, {5, 8.5}, {8.5, 8.5}, {8.5, 5}, {5, 5}}}}, + } + + for k, test := range tests { + geom := polygonSieve(test.geom, test.resolution) + if test.sieved != nil && geom != nil { + if area(geom) != area(test.sieved) { + t.Errorf("test: %d, expected: %f \ngot: %f", k, test.sieved, geom) + } + } else if test.sieved == nil && geom != nil { + t.Errorf("test: %d, expected: %f \ngot: %f", k, test.sieved, geom) + } else if test.sieved != nil && geom == nil { + t.Errorf("test: %d, expected: %f \ngot: %f", k, test.sieved, geom) + } + } +} + +func TestMultiPolygonSieve(t *testing.T) { + var tests = []struct { + geom [][][][2]float64 + resolution float64 + sieved [][][][2]float64 + }{ + // Lower single polygon resolution + 0: {geom: [][][][2]float64{{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}}}, resolution: float64(1), sieved: [][][][2]float64{{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}}}}, + // Higher single polygon resolution + 1: {geom: [][][][2]float64{{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}}}, resolution: float64(101), sieved: nil}, + // Nil input + 2: {geom: nil, resolution: float64(1), sieved: nil}, + // Low multi polygon resolution + 3: {geom: [][][][2]float64{{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}}, {{{15, 15}, {15, 20}, {20, 20}, {20, 15}, {15, 15}}}}, resolution: float64(1), sieved: [][][][2]float64{{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}}, {{{15, 15}, {15, 20}, {20, 20}, {20, 15}, {15, 15}}}}}, + // single hit on multi polygon + 4: {geom: [][][][2]float64{{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}}, {{{15, 15}, {15, 20}, {20, 20}, {20, 15}, {15, 15}}}}, resolution: float64(9), sieved: [][][][2]float64{{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}}}}, + // single hit on multi polygon + 5: {geom: [][][][2]float64{{{{0, 0}, {0, 10}, {10, 10}, {10, 0}, {0, 0}}}, {{{15, 15}, {15, 20}, {20, 20}, {20, 15}, {15, 15}}}}, resolution: float64(101), sieved: nil}, + } + + for k, test := range tests { + geom := multiPolygonSieve(test.geom, test.resolution) + if test.sieved != nil && geom != nil { + if len(geom) != len(test.sieved) { + t.Errorf("test: %d, expected: %f \ngot: %f", k, test.sieved, geom) + } + } else if test.sieved == nil && geom != nil { + t.Errorf("test: %d, expected: %f \ngot: %f", k, test.sieved, geom) + } else if test.sieved != nil && geom == nil { + t.Errorf("test: %d, expected: %f \ngot: %f", k, test.sieved, geom) + } + } +}