1

I have a table of points in my database that I'm showing in OpenLayers in 2 different ways:

  1. Retrieving via API as GeoJSON and then display in OL using GeoJSON format (white point)
  2. Retrieving via GeoServer which serves a vector tile layer and displays in OL using MVT format (pink point)

This is the result I get for the same point:

enter image description here The map is in EPSG:3857

The data in the database is in the CRS NZGD2000 (EPSG:4167) but with a shift in longitude of -160 degrees so a point has the coordinates POINT (10.9909915167 -45.0688542). --- (note the longitude of 10.9909915167 instead of 170.9909915167.)

To deal with the shifted longitudes we have defined a custom CRS (we named it EPSG:1) that basically copies the specs of EPSG:4167 (NZGD2000) and adds 160 to the Prime meridian, like so:

Definition of the CRS in OpenLayers using proj4:

proj4.defs("EPSG:1", "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs +pm=160");

Definition of the CRS in GeoServer (added to the bottom of users_projections/epsg.properties file):

1=GEOGCS["NZGD2000_160_Shift", DATUM["New_Zealand_Geodetic_Datum_2000_160_shift", SPHEROID["GRS 1980",6378137,298.257222101, AUTHORITY["EPSG","7019"]], TOWGS84[0,0,0,0,0,0,0], AUTHORITY["EPSG","6167"]], PRIMEM["Custom",160], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]]]

To display the GeoJSON data, we use the following code to tell OpenLayers the data is coming in EPSG:1 and to project it to EPSG:3857, so OL uses the proj4 definition we declared above to do the transformation.

const source = new VectorSource({
      features: new GeoJSON().readFeatures(features, { dataProjection: "EPSG:1", featureProjection: "EPSG:3857" }),
    });
    layer.setSource(source);

To display the data coming from the GeoServer vector tile we use the following piece of code retrieving the data in EPSG:900913, so the transformation is done by GeoServer with the epsg we added to its projections file.

const source = new VectorTileSource({
      url: `${geoserverURL}/geoserver/gwc/service/tms/1.0.0/myworkspace:marks@EPSG%3A900913@pbf/{z}/{x}/{-y}.pbf`,
      format: new MVT(),
      wrapX: true,
      tileGrid: createXYZ({ maxZoom: 19 }),
    });

These are the FlatCoordinates (in EPSG:3857) we get when clicking on the map and do getGeometry():

projected by OpenLayers (white point): [19034630.105876002, -5632367.691823326]

projected by GeoServer (pink point): [19034630.10393626, -5632367.698059561]

What makes us believe that the point projected by OpenLayers is the right one, is that when we transformed the coordinates of the database POINT (10.9909915167 -45.0688542) from EPSG:1 to EPSG:3857, this was the result 19034630.1059;-5632367.69182 (similar to the OpenLayers one).

For some reason, GeoServer is not transforming it properly!

Just one more evidence, when accessing the GeoServer tile that retrieves that point I get this geometry:

http:///geoserver/gwc/service/tms/1.0.0/myworkspace:marks@EPSG%3A900913@pbf/19/511167/188457.geojson

enter image description here

Not sure if the coordinates showing only 2 decimal places means the pbf that OpenLayers receives also has this number of decimal places...

I wonder what might be causing this offset of millimetres and if there is a good way to debug this issue.

I wish I could see the coordinates with more decimal places.

Vince
  • 20,017
  • 15
  • 45
  • 64
Falcoa
  • 335
  • 2
  • 9
  • 1
    If I was concerned about mm accuracy I would not be using EPSG:3857 to display my map, see https://gis.stackexchange.com/a/8674/79 for a discussion of why you probably can't believe those lat/lon values anyway – Ian Turton Aug 19 '21 at 07:24
  • What would you use then? I would need the features to be matching the position as they come from the same DB source. Thanks for the link to the discussion. – Falcoa Aug 19 '21 at 08:59
  • 2
    I would also avoid GeoJSON as there is an unnecessary convertion to WGS84 and back again inherent in the format. If mm accuracy is important then store and display the coordinates in EPSG:4167 (as they are collected?) – Ian Turton Aug 19 '21 at 10:01
  • 1
    If you aren't using a scanning electron microscope for collection, it seems unlikely that you really need nine decimal places of precision in meters. Most geodata is unlikely to be accurate to meters, much less centimeters. – Vince Aug 19 '21 at 13:12
  • This is topographical data (surveying), so yes, is accurate to mm. – Falcoa Aug 19 '21 at 20:35
  • As stated in the question, in this case, the challenge is that data at identical positions are coming from two sources, and the sources don't line up. Agree they don't need to be accurate to 9dp of a metre, but we can see from the example coordinates the mismatch between the coordinates is about 8mm - so 2dp with rounding. – nstillwell Aug 19 '21 at 20:49

0 Answers0