queryForGeoJSONFeaturesInTable throws an error for some geopackage files

m-schuetz commented 4 years ago

Version: 2.1.0

First of all, thanks for this comprehensive geopackage library!

I recently ran into a problem with some geopackage files. Calling queryForGeoJSONFeaturesInTable results in the following error:

Uncaught (in promise) TypeError: Cannot read property 'x' of null
    at geocentricToGeodetic (geopackage.js:58208)
    at datum_transform (geopackage.js:58393)
    at transform (geopackage.js:58506)
    at transformer (geopackage.js:58545)
    at proj4$1 (geopackage.js:58584)
    at BoundingBox.projectBoundingBox (geopackage.js:1264)
    at FeatureTableIndex.queryWithBoundingBox (geopackage.js:5724)
    at FeatureDao.queryForGeoJSONIndexedFeaturesWithBoundingBox (geopackage.js:13493)
    at module.exports.GeoPackage.queryForGeoJSONFeaturesInTable (geopackage.js:14250)
    at resolver (GeoPackageLoader.js:94)

The issue is, that point ends up being null in geocentricToGeodetic.

The point becomes null during datum_transform after calling geodeticToGeocentric.

And the reason is because the latitude is not in polar coordinates:

This is the location in my code where the error occurs: https://github.com/potree/potree/blob/72748db38fd46a98c4d0800a70b1fe2570365152/src/loader/GeoPackageLoader.js#L92

Note that this looks a bit different than the official examples and the GeoPackage Viewer, because the version: 2.1.0 release seems to have a different API (no getFeatureRow() function, for example). What I'd like to accomplish is to load 3D data of all stored features, so that I can create 3D models out of them.

Here is an example that works: 3D Scene: http://mschuetz.potree.org/geopackage/examples/geopackage.html Geopackage file: http://mschuetz.potree.org/geopackage/examples/morro_bay_shp/gpkg/geopackage.gpkg

And this is the geopackage that causes an error when loading with geopackage-js: 3D Scene: http://mschuetz.potree.org/geopackage/not_working/gpgp.html Geopackage file: http://mschuetz.potree.org/geopackage/not_working/test2d_potree17_epsg2056.gpkg

Help on how to fix this issue or on how to correctly load 3D geometry data for features would be highly appreciated. Thanks!

danielbarela commented 4 years ago

I will check into this and get back to you soon.

danielbarela commented 4 years ago

OK, so the problem is that the queryForGeoJSONFeaturesInTable takes a bounding box but does not specify that this bounding box should be a bounding box in 4326. So, when you got the bounding box from the FeatureDao it was in the EPSG:2056 projection. So when you pass that bounding box in to the queryForGeoJSONFeaturesInTable function, it ends up passing it to the FeatureDao.queryForGeoJSONIndexedFeaturesWithBoundingBox method. This method expects that the bounding box that was passed in should be in 4326. I believe that what you are doing is perfectly reasonable and so I am going to add in the ability to pass in the projection that your bounding box is in, that way the query will work properly.

In the meantime you can call this: boundingBox.projectBoundingBox(dao.projection, 'EPSG:4326') on line 92 of GeoPackageLoader just after you get the bounding box from the feature dao. When you pass that in, it should all work since that bounding box will now be in 4326.

Let me know if this fixes your issue, but I will add the ability to pass in the projection, or maybe I will add the projection that the bounding box is in, to the BoundingBox that comes back from the featureDao.getBoundingBox method. Either way, what you are doing is totally valid and I will ensure that it works in the future.

I will leave this open, in case there are other issues that come up due to the feature data being not in 4326.

m-schuetz commented 4 years ago

I will leave this open, in case there are other issues that come up due to the feature data being not in 4326.

Thanks for that! The workaround with projectBoundingBox works and I can now load the data!

There seems to be another issue where the feature geometry doesn't align with the map. For example, in the GeoPackage Viewer, the feature data from this geopackage file is displayed with a mismatch like this:

And if I load it in Potree, the mismatch looks like this:

I've actually been able to work around this issue in Potree by using the projection function from geopackage-js, rather than proj4.js, so this is how it's supposed to look:

The issue manifests like this:

My 3D coordinate system is the same as the point cloud coordinate system, in this case EPSG:2056.
geopackage.js, specifically data.queryForGeoJSONFeaturesInTable(table, boundingBox) returns coordinates in WGS84
I use proj4.js to convert the features from WGS84 to scene coordinates, using this proj4 string from spatialreference.org.
This works fine for http://mschuetz.potree.org/geopackage/examples/geopackage.html where the geopackage itself seems to be WGS84 to begin with
For the missmatched geopackage file, proj4 and geopackage's featureDao.projection functions return slightly different output coordinates for the same input coordinates.

m-schuetz commented 4 years ago

You can reproduce this projection discrepancy in the example page above if you set a breakpoint inside GeoPackageLoader.js line 92, let boundingBox = dao.getBoundingBox();

Then execute this code in the dev tools:

{
    let epsg_4326 = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs";
    let epsg_2056 = "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs";

    proj4.defs("WGS84", epsg_4326);
    proj4.defs("EPSG:2056", epsg_2056);

    let transform = proj4("WGS84", "EPSG:2056");

    let [long, lat] = [6.939912664130361, 46.99803776018691];

    console.log(dao.projection.forward([long, lat]))
    console.log(transform.forward([long, lat]))

}

output:

[2562000.000046756, 1205127.1312170895]
[2562067.9954865887, 1205341.1478612714]

Proj4 strings from

danielbarela commented 4 years ago

To create the proj4js projection, we will pass the WKT defined in the definition column of the gpkg_spatial_ref_sys table to proj4. Except in the case of 4326. So, it appears that the issue is, the definition in your GeoPackage is:
PROJCS["CH1903+ / LV95",GEOGCS["CH1903+",DATUM["CH1903+",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],TOWGS84[674.374,15.056,405.346,0,0,0,0],AUTHORITY["EPSG","6150"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4150"]],PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],PARAMETER["latitude_of_center",46.95240555555556],PARAMETER["longitude_of_center",7.439583333333333],PARAMETER["azimuth",90],PARAMETER["rectified_grid_angle",90],PARAMETER["scale_factor",1],PARAMETER["false_easting",2600000],PARAMETER["false_northing",1200000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","2056"]] which is very slightly different from (https://spatialreference.org/ref/epsg/ch1903-lv95/ogcwkt/):
PROJCS["CH1903+ / LV95",GEOGCS["CH1903+",DATUM["CH1903",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],TOWGS84[674.374,15.056,405.346,0,0,0,0],AUTHORITY["EPSG","6150"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4150"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Hotine_Oblique_Mercator"],PARAMETER["latitude_of_center",46.95240555555556],PARAMETER["longitude_of_center",7.439583333333333],PARAMETER["azimuth",90],PARAMETER["rectified_grid_angle",90],PARAMETER["scale_factor",1],PARAMETER["false_easting",2600000],PARAMETER["false_northing",1200000],AUTHORITY["EPSG","2056"],AXIS["Y",EAST],AXIS["X",NORTH]]

When I modify your GeoPackage file and change the definition in the srs table to the above, the features are drawn in the correct location.

ngageoint / geopackage-js

queryForGeoJSONFeaturesInTable throws an error for some geopackage files #139