Ticket #1790 (closed defect: fixed)

Opened 4 years ago

Last modified 4 years ago

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

Reported by: bartvde Owned by: 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 Download (340.7 KB) - added by mloskot 4 years ago.
Visualization of the test shapefile and spatial index. Everything looks OK but doesn't work OK. Hmm…
shptreevis-qgis.png Download (4.1 KB) - added by mloskot 4 years ago.
More clear shptree visualization

Change History

Changed 4 years ago by warmerdam

  • cc warmerdam added
  • 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 shape qix added

Mateusz,

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

Changed 4 years ago by bartvde

Any news on this one?

Changed 4 years ago 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.

Changed 4 years ago by mloskot

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

Changed 4 years ago by mloskot

More clear shptree visualization

Changed 4 years ago 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)

Changed 4 years ago 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.

Changed 4 years ago by mloskot

  • keywords shptree added
  • 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.

Note: See TracTickets for help on using tickets.