Predicates should work on Geometry not SpatialObject

[VladimirAlexiev]

Predicates are defined with domain=range=SpatialObject. That is an abstract class that includes Feature and Geometry.

Imagine a database of features and respective geometries, where 1k pairs satisfy sfWithin.

What is ?x geo:sfWithin ?y supposed to return?

• 1k Geometry pairs , or
• 1k Feature pairs, or
• 1k <Geometry, Feature> pairs, or
• ...
• 4k pairs <Geometry|Feature, Geometry|Feature>

Because of the union nature of SpatialObject, the last bullet is suggested.. but that's absurd.

And there's another problem: a Feature can have multiple Geometries:

• "equivalent" but using different SRS
• or even non-equivalent: centroid vs bbox vs detailed; time-dependent; at high/low tide, etc

Say you have F1.G1 (centroid), F1.G2 (polygon), F2.G1 (polygon). And say that F1.G1 is within F2.G1 but F1.G2 is not within F2.G1. What should "F1 within F2" return?

Unfortunately the standard gives no examples of predicate use. Nor does it say something like:
"If a Feature has only one Geometry then it can be used in or returned by spatial predicates"

It poses equivalence of predicates and functions, but how can that be if predicates work on SpatialObject and functions only on geomLiteral?

I suggest to change domain=range=Geometry.

[ashleysommer]

Because of the union nature of SpatialObject, the last bullet is suggested.. but that's absurd.

I agree. If the geo predicates are used as magic predicates to generate new solutions, then this is a problem, because it would be expected to generate 4k pairs, when we really want 1k Geometries.

Unfortunately the standard gives no examples of predicate use. Nor does it say something like: "If a Feature has only one Geometry then it can be used in or returned by spatial predicates"

This is because the geo predicates were not designed to be used as magic predicates that can dynamically produce new solutions.
They are designed to be used as pre-calculated assertions. The user designing the dataset can calculate geometry relationships (using whatever offline methods), then assert these relationships in the data, eg geom:1 geo:sfWithin geom:2 (geometry within geometry). Or bio:occurrence1 geo:sfWithin sites:field2 (feature within feature).

That is why the domain is Geometry|SpatialObject. It doesn't matter if those features have multiple geometries, the user will indicate the relationships based on their own business rules.

It poses equivalence of predicates and functions, but how can that be if predicates work on SpatialObject and functions only on geomLiteral?

This is one I'm not sure about. The predicates and the functions are not equal, because the functions only work on geomLiterals (probably for the reasons pointed out above), and they have strict rules about how they are implemented.

[VladimirAlexiev]

But @ashleysommer a "real geosparql implementation" doesn't rely on offline custom materialisation processes. A "real" implementation uses a geo index (eg quadtree) to quickly find candidate solutions, whether none, one or both args are unbound.

The reason is that there are too many relations to precompute: up to 10×N^2.

• Say in a concrete case you only need sfWithin and there are maybe 2×N (a few large polygons, the rest are points, contained in a couple of the polygons)
• but the totality of relations are 10-15 kinds and I believe you'll quickly get quadratic complexity
• when and what process decides to precompute another kind?
• and how do you handle updates?

I challenge you to prove that the spec suggests or demands precomputation.

@situx or @nicholascar please weigh in

[VladimirAlexiev]

Eg Qlever uses precomputing (materialization) at https://qlever.dev/osm-planet/.
This query picks two features connected by ogc:sfContains and shows their representation:

curl "https://qlever.dev/api/osm-planet?query=PREFIX+osmkey%3A+%3Chttps%3A%2F%2Fwww.openstreetmap.org%2Fwiki%2FKey%3A%3E%0APREFIX+osmrel%3A+%3Chttps%3A%2F%2Fwww.openstreetmap.org%2Frelation%2F%3E%0APREFIX+ogc%3A+%3Chttp%3A%2F%2Fwww.opengis.net%2Frdf%23%3E%0APREFIX+geo%3A+%3Chttp%3A%2F%2Fwww.opengis.net%2Font%2Fgeosparql%23%3E%0Aconstruct+%7B%3Fx+%3Fp+%3Fy%7D+WHERE+%7B%0A++values+%3Fx+%7B%0A++++%3Chttps%3A%2F%2Fwww.openstreetmap.org%2Fnode%2F10000263025%3E+%3Chttps%3A%2F%2Fosm2rdf.cs.uni-freiburg.de%2Frdf%2Fgeom%23osmnode_tagged_10000263025%3E%0A++++%3Chttps%3A%2F%2Fwww.openstreetmap.org%2Fnode%2F9451734825%3E+%3Chttps%3A%2F%2Fosm2rdf.cs.uni-freiburg.de%2Frdf%2Fgeom%23osmnode_tagged_9451734825%3E%0A++%7D%0A++%3Fx+%3Fp+%3Fy%0A%7D"|riot --formatted ttl -syntax ttl

Looks like this (shortened):

<https://www.openstreetmap.org/node/9451734825>
        a       <https://www.openstreetmap.org/node>;
        <http://www.opengis.net/ont/geosparql#hasGeometry>
                <https://osm2rdf.cs.uni-freiburg.de/rdf/geom#osmnode_tagged_9451734825>;
        <http://www.opengis.net/rdf#sfContains>
                <https://www.openstreetmap.org/node/10000263025>;
        <http://www.opengis.net/rdf#sfCovers>
                <https://www.openstreetmap.org/node/10000263025>;
        <http://www.opengis.net/rdf#sfEquals>
                <https://www.openstreetmap.org/node/10000263025>;
        <http://www.opengis.net/rdf#sfIntersects>
                <https://www.openstreetmap.org/node/10000263025> , ...;
        <https://www.openstreetmap.org/wiki/Key:addr:city>
                "South Windsor";
        <https://www.openstreetmap.org/wiki/Key:addr:housenumber>
                "46";
        <https://www.openstreetmap.org/wiki/Key:addr:postcode>
                "06074";
        <https://www.openstreetmap.org/wiki/Key:addr:state>
                "CT";
        <https://www.openstreetmap.org/wiki/Key:addr:street>
                "Kennedy Road";
        <https://www.openstreetmap.org/wiki/Key:addr:unit>
                "3" .

<https://www.openstreetmap.org/node/10000263025>
        a       <https://www.openstreetmap.org/node>;
        <http://www.opengis.net/ont/geosparql#hasGeometry>
                <https://osm2rdf.cs.uni-freiburg.de/rdf/geom#osmnode_tagged_10000263025>;
        <http://www.opengis.net/rdf#sfContains>
                <https://www.openstreetmap.org/node/9451734825>;
        <http://www.opengis.net/rdf#sfCovers>
                <https://www.openstreetmap.org/node/9451734825>;
        <http://www.opengis.net/rdf#sfEquals>
                <https://www.openstreetmap.org/node/9451734825>;
        <http://www.opengis.net/rdf#sfIntersects>
                <https://www.openstreetmap.org/node/9451734825> , ...
        <https://www.openstreetmap.org/wiki/Key:addr:city>
                "South Windsor";
        <https://www.openstreetmap.org/wiki/Key:addr:housenumber>
                "46";
        <https://www.openstreetmap.org/wiki/Key:addr:postcode>
                "06074";
        <https://www.openstreetmap.org/wiki/Key:addr:state>
                "CT";
        <https://www.openstreetmap.org/wiki/Key:addr:street>
                "Kennedy Road" .

<https://osm2rdf.cs.uni-freiburg.de/rdf/geom#osmnode_tagged_9451734825>
        <http://www.opengis.net/ont/geosparql#asWKT>
                "POINT(-72.585525 41.854633)"^^<http://www.opengis.net/ont/geosparql#wktLiteral> .

<https://osm2rdf.cs.uni-freiburg.de/rdf/geom#osmnode_tagged_10000263025>
        <http://www.opengis.net/ont/geosparql#asWKT>
                "POINT(-72.585525 41.854634)"^^<http://www.opengis.net/ont/geosparql#wktLiteral> .

So spatial relations are between Features not Geometries (like @ashleysommer suggests).
But Qlever use the ogc namespace (eg ogc:sfContains) not the standard geo namespace:

PREFIX ogc: <http://www.opengis.net/rdf#> # materialized
PREFIX geo: <http://www.opengis.net/ont/geosparql#> # standard

And that's a good thing because if their dump used the geo namespace and you load it to a repo that supports dynamic magic predicates, there'll be a lot fo confusion.

---

Now let's do some counts:

• https://qlever.dev/osm-planet/uG2e30
  Number of features=geometries (they don't use the GeoSPARQL classes):

SELECT (count(*) as ?c) WHERE {
  ?x geo:hasGeometry ?y
}
# 11,384,152,648

• https://qlever.dev/osm-planet/sm6CuK
  Number of materialized topological relations:

?p	?count_p	?percent	comment
sfIntersects	33,187,249,918	49.82
sfCovers	15,224,253,618	22.86	non-standard: ad-freiburg/osm2rdf#132
sfContains	14,982,122,627	22.49	inverse: sfWithin
sfTouches	2,182,091,666	3.28
sfCrosses	643,745,207	0.97
sfOverlaps	205,302,326	0.31
sfEquals	183,895,584	0.28
sfDisjoint	missing

The total is about 65B so 6x the number of geometries not quadratic.
However, they lack :sfDisjoint, and that is quadratic.
Any pair that's not related by the above relations must necessarily be sfDisjoint.

Now, surely sfDisjoint should not be used with entirely free variables.
But I can see a reasonable case eg "Given a region, find its disjoint decomposition into subregions",
I could write it something like this:

select * {
  bind(... as ?region)
  ?subregion geo:sfWithin ?region
  filter not exists {
    ?subregion2 geo:sfWithin ?region
    filter (!exists {?subregion2 geo:sfDisjoint ?subregion}
      && !exists {?subregion2 geo:sfWithin ?subregion})
  }
}

[mperry455]

The predicates for topological relations were not intended to be used as magic predicates. As @ashleysommer suggested, the intention was for them to be used to assert topological relations that could then be used for logical reasoning to infer additional topological relations.

The topological query functions are designed to be used with a spatial index to quickly find solutions that satisfy a spatial constraint.

[VladimirAlexiev]

@mperry455 I think you got it backward!

• Boolean functions cannot generate solutions, they can only check. You can't invoke a boolean function with any UNDEF params because then it will simply return UNDEF.

select (geof:sfWithin(?x,?y) as ?z) {}
# UNDEF

• You can only use boolean functions in FILTER. Again, if any of the args is unbound, it will return UNDEF and the filter will fail:

select * {
  bind("POINT(...)"^^geo:wktLiteral as ?y)
  filter(geof:sfWithin(?x,?y))
}
# no rows

• So you must somehow enumerate all candidate geometries and then filter them: that is quadratic complexity!
  Say to find caves in cities, you can use a query like this.
  If you have a million caves and a million cities, the Cartesian product will be 10^12 and you're screwed:

select * {
  ?x a geo:Feature; dct:type <feat/cave>; geo:hasGeometry/geo:asWKT ?xWkt.
  ?y a geo:Feature; dct:type <feat/city>; geo:hasGeometry/geo:asWKT ?yWkt.
  filter(geof:sfWithin(?xwkt,?yWkt))
}

• Furthermore, function arguments are geoLiteral. So even if they COULD generate candidate solutions, how will you then relate geoLiterals to Geometries to Features?

In contrast, magic predicates can generate candidate solutions.
If you rewrite the query this way and there are 2 caves per city on average, then it will return 2M results: a manageable affair.

select * {
  ?y a geo:Feature; dct:type <feat/city>.
  ?x geo:sfWithin ?y.
  ?x a geo:Feature; dct:type <feat/cave>.
}

---

Thinking of my original question "What is ?x geo:sfWithin ?y supposed to return?":
13. Query Rewrite Extension does define such combinatory logic (this is RIF):

Forall ?so1 ?so2 ?g1 ?g2 ?g1Serial ?g2Serial (
    ?so1[ogc:relation->?so2] :- Or (
        And (
            # feature - feature rule
            ?so1[geo:hasDefaultGeometry->?g1]
                ?so2[geo:hasDefaultGeometry->?g2]
                ?g1[ogc:asGeomLiteral->?g1Serial]
                ?g2[ogc:asGeomLiteral->?g2Serial]
                External(ogc:function(?g1Serial, ?g2Serial))
        )
        And (
            # feature - geometry rule
            ?so1[geo:hasDefaultGeometry->?g1]
                ?g1[ogc:asGeomLiteral->?g1Serial]
                ?so2[ogc:asGeomLiteral->?g2Serial]
                External(ogc:function(?g1Serial, ?g2Serial))
        )
        And (
            # geometry - feature rule
            ?so1[ogc:asGeomLiteral->?g1Serial]
                ?so2[geo:hasDefaultGeometry->?g2]
                ?g2[ogc:asGeomLiteral->?g2Serial]
                External(ogc:function(?g1Serial, ?g2Serial))
        )
        And (
            # geometry - geometry rule
            ?so1[ogc:asGeomLiteral->?g1Serial]
                ?so2[ogc:asGeomLiteral->?g2Serial]
                External(ogc:function(?g1Serial, ?g2Serial))
        )
    )
)

It's obvious to me that such a rule can only work with bound arguments ?so1, ?so2.

[ashleysommer]

@mperry455 is correct, (and he would know, he wrote the spec), the GeoSPARQL predicates are intended to be used for asserting known topological relationships between geospatial features.
The intention is to then use those relationships to infer or derive further topological relationships using SPARQL mechanisms.

Eg,

feature:A a geo:Feature .
feature:B a geo:Feature .
feature:C a geo:Feature .

# We _know_ feature B is within feature A.
feature:B geo:sfWithin feature:A .
# And we _know_ feature C is within feature B.
feature:C geo:sfWithin feature:B .

# The sparql query can allow us to discover that 
# feature C is within feature A.

SELECT ?f WHERE {
  feature:C sfWithin+ ?f
}

a "real geosparql implementation" doesn't rely on offline custom materialisation processes.

Say in a concrete case you only need sfWithin and there are maybe 2×N (a few large polygons, the rest are points, contained in a couple of the polygons)
but the totality of relations are 10-15 kinds and I believe you'll quickly get quadratic complexity

I was the lead software engineer on the LOC-I project in Australia in 2018, and yes that is exactly what we did. We took in four large geospatial datasets, computed the topological relationships between a) geospatial features within each dataset, and b) top-level geospatial features between datasets.

We then loaded these into a triplestore, and used SPARQL to demonstrate the ability to derive all topological relationships between all of the geospatial features.

Back then there was no concept of "magic" predicates, which could calculate the topological relationships between features based on their Geometry literals, and produce new solutions at runtime. And the GeoSPARQL v1.0 spec was certainly not written with that in mind.

Your argument in this thread - that the predicates should work on Geometry not SpatialObject - is valid, and should be considered. But I'm giving you the context for why it is that way. The design was for the topological relationships betwen geospatial features to be properties of the Features themselves, not the geometry of the feature.