RTreeWithin not working

[Sacha]

Hello,

I have recently updated from using SpatiaLite 2.3, initially to 2.4
and now trying 3.0. I was using the mbrCache with 2.3 but have now
updated my databases to use spatialindex and query using RTreeWithin/
RTreeIntersects. Unfortunately not all my spatial queries that
previously worked with the MbrWithin are working using RTreeWithin.

One example test I have effectively queries a UTM zone (although my
coordinates are stored in lat/long). There are a large number of
polygons well within the zone (verified the coordinates in the
spatialite gui), but RTreeWithin query for that domain returns
nothing. If I edit the sql query slight to use the old form:

SELECT pkid FROM idx_objectstore_geom WHERE
xmin >= x1
and ymin >= y1
and xmax <= x2
and ymax <= y2

(x1, y1 and x2, y2 are the south west and north east corners of my
query bounding box)

I get the results as expected, whereas the following returns nothing:

SELECT pkid FROM idx_objectstore_geom WHERE RTreeWithin(x1, y1, x2,
y2)

I am not sure if I am misunderstanding something about the R*tree
methods, but assumed that under the hood it essentially did the same
as the old method.

I have also seen some issues related to floating-point accuracy as
discussed elsewhere, but I can probably get round this by adjusting
the way I query (although these didn't seem to be an issue with 2.3
and MbrCache).

Thanks in advance for any thoughts on this.

I am working in C++ on Windows using the statically linked library
with MSVC.

[Roman]

Sacha,

this looks like the problem I described in previous thread (http://groups.google.com/group/spatialite-users/browse_thread/thread/93ccddb4f148256e). Try to add some margin to requested bounding box.

RTreeWithin(x1-margin, y1-margin, x2+margin,y2+margin)

where margin for lat/long dataset should be 0.00001.

Roman

[a.fu...@lqt.it]

> SELECT pkid FROM idx_objectstore_geom WHERE
> xmin >= x1
> and ymin >= y1
> and xmax <= x2
> and ymax <= y2
>

> SELECT pkid FROM idx_objectstore_geom
> WHERE RTreeWithin(x1, y1, x2, y2)
>

not at all the same: they are not equivalent.
because Within implies "geom2 is *fully* within
geom1". So, if some geometry actually crosses the
BBOX will be excluded by this spatial filter.
http://www.postgis.org/docs/ST_Within.html

the equivalent spatial filter is RTreeContains:
this will extract both geometries fully
contained within the BBOX *and* geometries
crossing the BBOX itself
in other words, when accessing the R*Tree
as a Spatial Index (quick BBOX filtering)
you are always expected to use RTreeContains.

bye Sandro

[Roman]

So, if some geometry actually crosses the
BBOX will be excluded by this spatial filter.
http://www.postgis.org/docs/ST_Within.html

the equivalent spatial filter is RTreeContains:
this will extract both geometries fully
contained within the BBOX *and* geometries
crossing the BBOX itself
in other words, when accessing the R*Tree
as a Spatial Index (quick BBOX filtering)
you are always expected to use RTreeContains.

Really? I think you meant RTreeIntersects. Requested BBOX is the second parameter (geometry B), so RTreeContains extracts geometries whose MBR contains requested BBOX and does not extract geometries inside requested BBOX.

Roman

[Sacha]

Thanks for the replies.

Roman, the bounding box expansion is not an issue as in my current
test I am querying (-5.999999,-80,-1e-006,84) with an example mbr in
my index table being (-3.802227, 51.624760, -3.801922, 51.624966).
This should be well within the query!

On Aug 8, 5:36 pm, a.furi...@lqt.it wrote:
> not at all the same: they are not equivalent.
> because Within implies "geom2 is *fully* within
> geom1". So, if some geometry actually crosses the
> BBOX will be excluded by this spatial filter.http://www.postgis.org/docs/ST_Within.html

Sandro., I do want to query geometries completely within the bbox. I
use the contains method when I query a point and want all geometries
that contain that point. In this case I simply want all the geometries
completely within the bbox as the values above, which in my case
should return hunderds of geometries. The definition from the
spatialite documentation of RTreeWithin says:

"Retrieves from an R*Tree Spatial Index any entity whose MBR falls
within the rectangle identified by extreme points x1 y1 and x2 y2"

which is exactily what I want in this case.

[Roman]

Sacha,

now I see it, your query should be

SELECT pkid FROM idx_objectstore_geom WHERE pkid MATCH RTreeWithin(x1, y1, x2,y2)

For explanation, read this: http://groups.google.com/group/spatialite-users/browse_thread/thread/a8fdbe35e44ce8eb

[Sacha]

Yes. Sorry, I realised i missed the "pkid MATCH" part in my initial
post, but that is exactly what I have been doing.

[Sacha]

Having done some further investigations today I had some interesting
findings. My original database has geometries in several cities around
the world, spread over several different UTM zones. When doing my
query in zone 30 I get no results. I gradually cut the geometries in
other UTM zones out of the database. Each time there were still no
results returned, until I removed geometries from all the other UTM
zones. In other words, my query works when my database only has
geometries within the UTM zone I am querying.

I also tried querying individual UTM zones for the other geometry data
and got no results. Perhaps there is a problem the way the bounding
box being queried with is initially being checked against the largest
bounding box of all spatial data in the database? It is probably worth
doing some testing with smaller bounding box queries.

Any thoughts on what the issue could be?

[Dan Putler]

Hi, This should be a simple question. I'm trying to do a join between two tables that in some case is many to one, and the geometry is in the "one" table. To make this more concrete, From these two tables I want to create a third. The structure of the two original tables are: Table1   ID: Integer   Name: Text   Type: Text Table2   ID: Integer  Geometry: Blob geometry column in my case Table2 is from a line string shapefile (road segments), and Table1 from a dbf. The ID and geometries are unique in Table2, but some of the values of the ID are not unique in Table1, and Table1 has more records than Table2. The code I'm using to create the third table (table3) is: CREATE TABLE table3(   ID INTEGER,  Name TEXT,  Type TEXT ); SELECT AddGeometryColumn( 'table3', 'Geometry', 4269, 'XY'); INSERT INTO table3   SELECT a.ID, a.NAME, a.TYPE, ST_LineFromWKB(b.Geometry, 4269)   FROM table1 as a, table2 as b   WHERE a.ID = b.ID; I'm getting no errors from the spatialite command line console. When I read the new table into QGIS, the program indicates I'm getting the number of records in the file table I expect, but the layer does not display in QGIS (Table2 does display correctly). Looking at the metadata for the new data in QGIS indicates that the new table has no extent information (which is not true for Table2). My real question is whether I'm treating the geometry column correctly. I've looked at the spatialite manual and the cookbook, and haven't come across an example exactly like this one, so I'm concerned I may be error, which is causing the QGIS issues. Dan

[a.fu...@lqt.it]

Hi Dan,

You cannot use ST_LineFromWKB() this way; this is
wrong and return an invalid Geometry (i.e. NULL).
WKB *is not* the internal geometry format used
by SpatiaLite.

Just to be pedantic, this is valid:
ST_LineFromWKB(ST_AsBinary(b.Geometry, 4269))
- using ST_AsBinary to convert from "internal"
� to WKB
- then using ST_LineFromWKB to get back again
� the same "internal" geometry from WKB
Obviously useless: you already have the "internal"
geometry. there is no need at all to convert back
and forward twice using WKB

And I notice a second potential flaw: in your
"table3" there is no PRIMARY KEY.
Here is your SQL rewritten in a better way:

CREATE TABLE table3(
� pk_id INTEGER PRIMARY KEY AUTOINCREMENT,
� ID INTEGER,

� Name TEXT,
� Type TEXT
);

SELECT AddGeometryColumn( 'table3', 'Geometry', 4269, 'XY');

INSERT INTO table3 (pk_id, ID, Name, Type)
� SELECT NULL, a.ID, a.NAME, a.TYPE, b.Geometry
� FROM table1 as a, table2 as b
� WHERE a.ID = b.ID;

bye Sandro

[a.fu...@lqt.it]

Hi Sasha,

I'm unable to understand the exact terms of
your problem.

- I understand you surely have some worldwide dataset:
and I easily imagine this is SRID=4326, WGS84, long/lat

- then you are attempting to retrieve any data located
within a single UTM zone: and this probably implies
some appropriate CRS transformation.
I strongly suspect this is someway related to your
problem: but I need some SQL snippet in order to
give you any further evaluation.

bye Sandro

[Sacha]

Hi Sandro,

On Aug 10, 10:12 am, a.furi...@lqt.it wrote:
> - I understand you surely have some worldwide dataset:
>  and I easily imagine this is SRID=4326, WGS84, long/lat

Correct.

> - then you are attempting to retrieve any data located
>  within a single UTM zone: and this probably implies
>  some appropriate CRS transformation.
>  I strongly suspect this is someway related to your
>  problem: but I need some SQL snippet in order to
>  give you any further evaluation.

Here are the queries used to create my table (it's somewhat
complicated due to some bespoke data structures so I add further
columns required after these queries):

CREATE TABLE objectstore(key INTEGER NOT NULL PRIMARY KEY);
SELECT AddGeometryColumn('objectstore'', 'geom', 4326, 'POLYGON',
'XYZ');
SELECT CreateSpatialIndex('objectstore', 'geom')
(.....add other columns....)

I then insert my database entires using the C API to create geometry
blobs of polygons. When doing my spatial query I use the following:

SELECT key,subtype FROM objectstore WHERE ROWID IN (SELECT rowid FROM
idx_objectstore_geom WHERE pkid MATCH
RTreeWithin(-5.999999,-80,-1e-006,84))

As I mentioned in my previous post, this returns nothing for my
worldwide dataset. If I remove data located outside the above queried
bounding box, then I get all the geometries returned as expected.

I can successfully view my full database in QGIS so I am sure that all
the expected geometries are succesfully added to the database and are
valid.

[a.furieri]

Roman, you are absolutely right.
I actually intended *RTreeIntersects*
sorry for this typo, causing further confusion :-(

bye Sandro

[Sacha]

RTreeIntersects works fine in place of RTreeWithin for my above query,
returning all the expected results. This is fine in this case where I
have a very large bounding box query, but I want to be able to perform
queries where only geometries completely within the bounding box are
retrieved and any overlapping/intersecting geometries are excluded.

[a.furieri]

Hi Sasha,

I've done some testing using GeoNames (about 7 millions long/lat
points);
and the following queries surely work as expected.
Probably there is some bug affecting the rarely used RTreeWithin:
the most usual way (and thorougly tested) is using RTreeIntersects.

SELECT Count(*)
FROM "GeoNames"
WHERE ROWID IN (
SELECT pkid FROM idx_GeoNames_Geometry
WHERE pkid
MATCH RTreeIntersects(6,-90,12,90));

SELECT Count(*)
FROM "GeoNames"
WHERE ROWID IN (
SELECT pkid FROM idx_GeoNames_Geometry
WHERE pkid
MATCH RTreeIntersects(-180,-90,-174,90));

bye sandro

[Sacha]

Thanks for the response Sandro.

Will you be raising a defect for the RTreeWithin? I need to use
RTreeWithin as well as RTreeIntersects and RTreeContains.

[Dan Putler]

Hi Sandro,

Sorry for starting a new thread, I seemed to have stepped onto another thread with my first message.

Using the INSERT statement you sent:

INSERT INTO table3 (pk_id, ID, Name, Type)
SELECT NULL, a.ID, a.NAME, a.TYPE, b.Geometry
FROM table1 as a, table2 as b
WHERE a.ID = b.ID;

generates the error message "SQL Error: 5 values for 4 columns"

The natural change to make to the query is:

INSERT INTO table3 (pk_id, ID, Name, Type, Geometry)

SELECT NULL, a.ID, a.NAME, a.TYPE, b.Geometry
FROM table1 as a, table2 as b
WHERE a.ID = b.ID;

This results in the error message: "SQL Error: table3.Geometry violates Geometry constraint [geom-type or SRID not allowed]"

Given this failure, I tried:

INSERT INTO table3 (pk_id, ID, Name, Type, Geometry)
SELECT NULL, a.ID, a.NAME, a.TYPE, ST_LineFromWKB(ST_AsBinary(b.Geometry, 4269))

FROM table1 as a, table2 as b
WHERE a.ID = b.ID;

which again results in the Geometry constraint error. If it helps, I can send along the original files and the code I can ran in an off-list email.

Dan

[a.furieri]

Hi Dan,

just to recapitulate, attempting to avoid any
further typo this time :-)

INSERT INTO table3 (pk_id, ID, Name, Type, Geometry)
SELECT NULL, a.ID, a.NAME, a.TYPE, b.Geometry
FROM table1 as a, table2 as b
WHERE a.ID = b.ID;

(this is fine, mine previous was obviously wrong)

?? SELECT AddGeometryColumn('table3', 'Geometry', 4269, 'XY');
this is indisputably wrong (missing geometry class).
I suppose you actually intended:

SELECT AddGeometryColumn('table3', 'Geometry', 4269,

'LINESTRING', 'XY');

please note: this assumes that b.Geometry always
contains 2D [XY] LINESTRING SRID=4269
if such assumption fails to be verified you'll
surely get a Geometry constraint violation.
performing a preliminary sanity check is very simple:

SELECT DISTINCT ST_GeometryType(geometry),
CoordDimension(geometry), ST_Srid(geometry)
FROM table2;

bye Sandro

[Dan Putler]

Hi Sandro,

Thanks for your time and very quick response, not to mention all the effort you have put into SpatiaLite. The sanity check was a critical point. I incorrectly assumed that the importer would determine the SRID of a shapefile using the *.prj component of the set, but I now know that the SRID needs to be explicitly specified on import.

Perhaps a "recipe" along the lines of what I was trying to do would be a nice addition to the SpatiaLite cookbook when it gets revised.

Dan

--- On Wed, 8/10/11, a.furieri <a.fu...@lqt.it> wrote:

> --
> You received this message because you are subscribed to the
> Google Groups "SpatiaLite Users" group.
> To post to this group, send email to spatiali...@googlegroups.com.
> To unsubscribe from this group, send email to
> spatialite-use...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/spatialite-users?hl=en.
>
>