spatial function support

[kttary]

I was wondering if there is a plan to add spatial function to geopackage? I'm developing an offline application rely on geopackage library but miss some spatial function such as:

1. ST_Transform to perform coordinate system transformation between projection system
2. ST_Area to calculate area
3. ST_Intersect to find some feature within given bounding box
4. ST_PointOnSurface to find a point on 2D Surface
5. ST_Distance Thank you in advance

[bosborn]

1: Handled by simple-features-proj-ios, see SFPProjectionTransform transformWithGeometry method

2: See GPKGGeometryUtils computeAreaOfDegreesPath method

3: See GPKGTileBoundingBoxUtils isPoint:inBoundingBox methods and GPKGFeatureOverlayQuery queryFeaturesWithBoundingBox methods

4 & 5: Handled by simple-features-ios, see SFGeometryUtils. Also GPKGTileBoundingBoxUtils

[kttary]

Nice, though i'm stumbling at number 2.

let pt1: SFPoint = SFPoint(hasZ: false, andHasM: false, andX:  -90.000, andY: 38)
let pt2: SFPoint = SFPoint(hasZ: false, andHasM: false, andX:  -90.001, andY: 38.001)
let pt3: SFPoint = SFPoint(hasZ: false, andHasM: false, andX: -90.002, andY: 38)
let pt4: SFPoint = SFPoint(hasZ: false, andHasM: false, andX:  -90, andY: 38)
let area: Double = GPKGGeometryUtils.computeArea(ofDegreesPath: [pt1, pt2, pt3, pt4])
print(area)

and swift throw runtime error: libswiftCore.dylib`-[_SwiftNativeNSArrayBase retain]: 0x10d897700 <+0>: pushq %rbp 0x10d897701 <+1>: movq %rsp, %rbp 0x10d897704 <+4>: pushq %rbx 0x10d897705 <+5>: pushq %rax 0x10d897706 <+6>: movq %rdi, %rbx -> 0x10d897709 <+9>: callq 0x10d852e20 ; swift_retain 0x10d89770e <+14>: movq %rbx, %rax 0x10d897711 <+17>: addq $0x8, %rsp 0x10d897715 <+21>: popq %rbx 0x10d897716 <+22>: popq %rbp 0x10d897717 <+23>: retq
0x10d897718 <+24>: nopl (%rax,%rax)

[bosborn]

You found an infinite recursion bug, now fixed here.

For now you can change your call to this: let area: Double = fabs(GPKGGeometryUtils.computeSignedArea(ofDegreesPath: [pt1, pt2, pt3, pt4]))

[kttary]

Thank you, its Works. I'm continuing to number 3 with following code:

let myBB: GPKGBoundingBox = GPKGBoundingBox(minLongitude: -90.0025, andMinLatitude: 37.999, andMaxLongitude: -89.999, andMaxLatitude: 38.0025)
let featureTables : NSArray = geoPackage.getFeatureTables() as NSArray;
let featureName: String = featureTables.object(at: featureTables.index(of: "Parks")) as! String;
let featureTable: GPKGFeatureDao = geoPackage.getFeatureDao(withTableName: featureName);
let myTiles:GPKGFeatureTiles = GPKGFeatureTiles(featureDao: featureTable)
let myOverlay:GPKGFeatureOverlay = GPKGFeatureOverlay(featureTiles: myTiles)
let myQuery:GPKGFeatureOverlayQuery = GPKGFeatureOverlayQuery(featureOverlay: myOverlay)
let myRes : GPKGFeatureIndexResults = myQuery.queryFeatures(with: myBB)
print(myRes.count())

given code above, geopackage throw an error "Index Manager is not set on the Feature Tiles and is required to query indexed features". Looks like index manager of Parks is nil. Am i miss some lines of code or need to set spatial index?

[bosborn]

The FeatureOverlayQuery and FeatureTiles might actually be a little high up if you just want to query. We use those more when we draw tiles from features once or on demand, and still want to query the features. So you can either use the index manager or manually query all and filter. I've included examples of both below. Even though your bounding box is in the same projection as the stored data, I added it so you can see how to query when it is not (or all the time to be safe).

let myBB: GPKGBoundingBox = GPKGBoundingBox(minLongitude: -90.0025, andMinLatitude: 37.999, andMaxLongitude: -89.999, andMaxLatitude: 38.0025)
let featureTables : NSArray = geoPackage.getFeatureTables() as NSArray
let featureName: String = featureTables.object(at: featureTables.index(of: "Parks")) as! String
let featureTable: GPKGFeatureDao = geoPackage.getFeatureDao(withTableName: featureName)

let myBBProjection: SFPProjection = SFPProjectionFactory.projection(withAuthority: PROJ_AUTHORITY_EPSG, andIntCode: PROJ_EPSG_WORLD_GEODETIC_SYSTEM)

// Indexed query
let indexer:GPKGFeatureIndexManager = GPKGFeatureIndexManager(geoPackage: geoPackage, andFeatureDao: featureTable)
indexer.indexLocation = GPKG_FIT_GEOPACKAGE // to index in GeoPackage or GPKG_FIT_METADATA for external metadata
indexer.index()
let myRes:GPKGFeatureIndexResults = indexer.query(with: myBB, andProjection: myBBProjection)
print(myRes.count())
myRes.close()

let transform:SFPProjectionTransform = SFPProjectionTransform(from: featureTable.projection, andTo: myBBProjection)

// Unindexed query
var count = 0
let allRows:GPKGResultSet = featureTable.queryForAll()
while(allRows.moveToNext()){
    let row:GPKGFeatureRow = featureTable.getFeatureRow(allRows)
    let geomData:GPKGGeometryData! = row.getGeometry()
    if(geomData != nil){
        let geometry:SFGeometry! = geomData.geometry
        if(geometry != nil){
            var envelope:SFGeometryEnvelope! = geomData.envelope
            if(envelope == nil){
                envelope = SFGeometryEnvelopeBuilder.buildEnvelope(with: geometry)
            }
            let geometryBoundingBox:GPKGBoundingBox = GPKGBoundingBox(geometryEnvelope: envelope)
            let transformBoundingBox = geometryBoundingBox.transform(transform)
            if(GPKGTileBoundingBoxUtils.overlap(with: myBB, andBoundingBox: transformBoundingBox) != nil){
                count = count + 1
            }
        }
    }
}
print(count)
allRows.close()

[kttary]

Perfect. Thank you for the usefull snippet. For point 4, it is slighty different. ST_POINTONSURFACE return a geometry that is guarantee lies on the polygon surface while some function on SFGeometryUtil return boolean. However, i don't heavily need that function. So it is enough for now.

[bosborn]

Ah, no method matching that behavior right now. There is a centroid method you can use, but the point is the centroid and may be outside of the surface depending on the polygon.