Skip to content

Commit

Permalink
Merge pull request #8 from PDOK/tests-and-logging
Browse files Browse the repository at this point in the history 
Tests and logging
  • Loading branch information
@WouterVisscher
WouterVisscher authored May 27, 2021
2 parents 560e26a + 35d1454 commit 0b417eb
Show file tree
Hide file tree
Showing 2 changed files with 221 additions and 83 deletions.
183 changes: 100 additions & 83 deletions sieve.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
}

Expand All @@ -102,17 +102,16 @@ 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
}

// getSourceTableInfo collects the source table information
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

Expand All @@ -122,7 +121,7 @@ func getSourceTableInfo(h *gpkg.Handle) []table {
var srsID int
err := rows.Scan(&t.name, &t.gcolumn, &gtype, &srsID)
if err != nil {
log.Fatal(err)
log.Fatalf("error ready the source table information: %s", err)
}

t.columns = getTableColumns(h, t.name)
Expand All @@ -131,31 +130,22 @@ func getSourceTableInfo(h *gpkg.Handle) []table {

tables = append(tables, t)
}
defer rows.Close()
return tables
}

// getSpatialReferenceSystem extracts this based on the given SRS id
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
}

Expand All @@ -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
}

Expand All @@ -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{
Expand All @@ -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
Expand All @@ -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() {
Expand All @@ -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{}
Expand All @@ -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:
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -325,69 +314,40 @@ 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
}
}

if ext == nil {
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 {
Expand All @@ -398,33 +358,53 @@ 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()

tables := getSourceTableInfo(srcHandle)

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
Expand All @@ -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)

Expand All @@ -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)
Expand All @@ -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]
Expand Down
Loading

0 comments on commit 0b417eb

Please # to comment.