after creating qix file no features are read from 3D point shapefile

[bartvde]

Before creating a qix file with shptree I get features using the following ogrinfo command:

ogrinfo D:\bart\dtb\dtb_subset\dtb2000_sym.shp dtb2000_sym -spat 60000 400000 90000 410000

When I create the qix file with shptree, and rerun the exact same command, I get no features returned anymore:

shptree D:\bart\dtb\dtb_subset\dtb2000_sym.shp

D:\bart\FWTools1.3.6>ogrinfo D:\bart\dtb\dtb_subset\dtb2000_sym.shp dtb2000_sym -spat 60000 400000 90000 410000 INFO: Open of `D:\bart\dtb\dtb_subset\dtb2000_sym.shp'

using driver `ESRI Shapefile' successful.

Layer name: dtb2000_sym Geometry: 3D Point Feature Count: 0 Extent: (50000.000000, 390000.000000) - (99997.706000, 419998.443000) Layer SRS WKT: PROJCS["Rijksdriehoekstelsel_New",

GEOGCS["GCS_Amersfoort",

DATUM["Amersfoort",

SPHEROID["Bessel_1841",6377397.155,299.1528128]],

PRIMEM["Greenwich",0.0], UNIT["Degree",0.017453292519943295]],

PROJECTIONDouble_Stereographic?, PARAMETER["False_Easting",155000.0], PARAMETER["False_Northing",463000.0], PARAMETER["Central_Meridian",5.38763888888889], PARAMETER["Scale_Factor",0.9999079], PARAMETER["Latitude_Of_Origin",52.15616055555555], UNIT["Meter",1.0]]

CTE: String (8.0)

I believe the same problem is present in OGR and Mapserver.

I tested with latest FWTools 1.3.6 on Windows.

Dataset can be downloaded from (about 2 Mb):  http://www.osgis.nl/download/dtb_subset.zip

[warmerdam]

Mateusz,

Could you try to reproduce and track this down fairly soon?

[bartvde]

Any news on this one?

[mloskot]

Unfortunately, not yet. I'm including this ticket in my today tasks.

[mloskot]

Visualization of the test shapefile and spatial index. Everything looks OK but doesn't work OK. Hmm...

[mloskot]

More clear shptree visualization

[mloskot]

Using xmin=50000 and ymin=390000 gives following result:

D:\dev\gdal\bugs\1790>ogrinfo dtb2000_sym.shp dtb2000_sym -spat 50000 390000 90000 410000 -so
INFO: Open of `dtb2000_sym.shp'
      using driver `ESRI Shapefile' successful.

Layer name: dtb2000_sym
Geometry: 3D Point
Feature Count: 68887
Extent: (50000.000000, 390000.000000) - (99997.706000, 419998.443000)
Layer SRS WKT:
PROJCS["Rijksdriehoekstelsel_New",
    GEOGCS["GCS_Amersfoort",
        DATUM["Amersfoort",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Double_Stereographic"],
    PARAMETER["False_Easting",155000.0],
    PARAMETER["False_Northing",463000.0],
    PARAMETER["Central_Meridian",5.38763888888889],
    PARAMETER["Scale_Factor",0.9999079],
    PARAMETER["Latitude_Of_Origin",52.15616055555555],
    UNIT["Meter",1.0]]
CTE: String (8.0)

[mloskot]

Second observation, nodes in the index are nested but the code in shptree.c:738 returns numshapes equal to 0 for each node:

fread( &numshapes, 4, 1, fp );

If I use slightly different extent (50000 390000 90000 410000), fread() gets 127776 features on level with extent (50000 390000 50011 390006). It looks to me like the index is generated incorrectly. Will work on it tomorrow.

[mloskot]

I have found that this issue is directly related to the problem reported in #1594 which has been already fixed in the trunk.

All versions of shptree utility, also this one in cluded in FWTools (version <=1.4.1), generate broken index when big value of depth is estimated.

Below, I've included short script from tests (observe number of features):

• Query features without using spatial index

ogrinfo dtb2000_sym.shp dtb2000_sym -spat 60000 400000 90000 410000 -so
OGR: OGROpen(dtb2000_sym.shp/003B70F0) succeeded as ESRI Shapefile.
INFO: Open of `dtb2000_sym.shp'
      using driver `ESRI Shapefile' successful.
OGR: GetLayerCount() = 1

Layer name: dtb2000_sym
Geometry: 3D Point
Feature Count: 39100
...

• Generate index using patched shptree algorithm, through ogrinfo:

ogrinfo -sql "CREATE SPATIAL INDEX ON dtb2000_sym" dtb2000_sym.shp

• Query features using spatial index

ogrinfo dtb2000_sym.shp dtb2000_sym -spat 60000 400000 90000 410000 -so
OGR: OGROpen(dtb2000_sym.shp/003B70F0) succeeded as ESRI Shapefile.
INFO: Open of `dtb2000_sym.shp'
      using driver `ESRI Shapefile' successful.
OGR: GetLayerCount() = 1

Layer name: dtb2000_sym
Geometry: 3D Point
SHAPE: Used spatial index, got 39221 matches.
Feature Count: 39100
...

The patched shptree has been applied to:

• trunk
• branches/1.4
• official Shapelib repository

I'm closing this bug as fixed.

Bart, I apologise you've been waiting so long for it.