Ticket #1790 (closed defect: fixed)

Opened 10 months ago

Last modified 8 months ago

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

Reported by: bartvde Assigned to: mloskot
Priority: highest Milestone: 1.4.3
Component: OGR_SF Version: svn-trunk
Severity: normal Keywords: shape qix shptree
Cc: warmerdam

Description

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

Attachments

layer-and-index-vis.png (340.7 kB) - added by mloskot on 10/25/07 00:11:31.
Visualization of the test shapefile and spatial index. Everything looks OK but doesn't work OK. Hmm…
shptreevis-qgis.png (4.1 kB) - added by mloskot on 10/25/07 00:38:55.
More clear shptree visualization

Change History

09/04/07 10:14:58 changed by warmerdam

  • cc set to warmerdam.
  • component changed from default to OGR_SF.
  • owner changed from warmerdam to mloskot.
  • version changed from unspecified to svn-trunk.
  • milestone set to 1.4.3.
  • keywords set to shape qix.

Mateusz,

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

09/18/07 05:22:14 changed by bartvde

Any news on this one?

09/18/07 05:51:41 changed by mloskot

  • priority changed from normal to highest.
  • status changed from new to assigned.

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

10/25/07 00:11:31 changed by mloskot

  • attachment layer-and-index-vis.png added.

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

10/25/07 00:38:55 changed by mloskot

  • attachment shptreevis-qgis.png added.

More clear shptree visualization

10/25/07 00:45:45 changed by 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)

10/25/07 01:09:17 changed by 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.

10/27/07 21:49:28 changed by mloskot

  • keywords changed from shape qix to shape qix shptree.
  • status changed from assigned to closed.
  • resolution set to fixed.

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.