Opened 17 years ago
Closed 16 years ago
#1790 closed defect (fixed)
after creating qix file no features are read from 3D point shapefile
Reported by: | bartvde | Owned by: | Mateusz Łoskot |
---|---|---|---|
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 (2)
Change History (8)
comment:1 by , 17 years ago
Cc: | added |
---|---|
Component: | default → OGR_SF |
Keywords: | shape qix added |
Milestone: | → 1.4.3 |
Owner: | changed from | to
Version: | unspecified → svn-trunk |
comment:3 by , 17 years ago
Priority: | normal → highest |
---|---|
Status: | new → assigned |
Unfortunately, not yet. I'm including this ticket in my today tasks.
by , 16 years ago
Attachment: | layer-and-index-vis.png added |
---|
Visualization of the test shapefile and spatial index. Everything looks OK but doesn't work OK. Hmm...
comment:4 by , 16 years ago
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)
comment:5 by , 16 years ago
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.
comment:6 by , 16 years ago
Keywords: | shptree added |
---|---|
Resolution: | → fixed |
Status: | assigned → closed |
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.
Mateusz,
Could you try to reproduce and track this down fairly soon?