Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

GeoTIF files for UK don't yield expected results when used in GDAL #2360

Closed
dfileccia opened this issue Sep 23, 2020 · 16 comments
Closed

GeoTIF files for UK don't yield expected results when used in GDAL #2360

dfileccia opened this issue Sep 23, 2020 · 16 comments

Comments

@dfileccia
Copy link

Since I have not found a place to ask questions:

I am using PROJ.4 7.1.1 and GDAL 3.1.3 trying to read the updated proj-data files to get the orthometric values for the UK since the toolkit provided by the national survey doesn't play well with our system. According to the README file there should only be an error of 0.00076m error from the source files but I am seeing 1m+ errors.

With other files (Australia, US, etc), I am seeing <= 1cm error which we would consider acceptable since we are a non-survey application.

My code to read the files is pretty simple for my test:

GDALAllRegister ();

GDALDataset* dset = (GDALDataset*) GDALOpen ("uk_os_OSGM15_GB.tif", GA_ReadOnly);
if (!dset) return;

const int nx = dset->GetRasterXSize ();
const int ny = dset->GetRasterYSize ();

double adfTransform[6];
dset->GetGeoTransform (adfTransform);

GDALRasterBand* band = dset->GetRasterBand (1);
if (!band) return;

while (read test data from file) {
   // Compute the row and column to read
   const int col = (lon - adfTransform[0]) / adfTransform[1];
   const int row = (adfTransform[3] - lat) / adfTransform[5];

   // Read a 2x2 cell and interpolate in-between
   float gsep[4] { 999, 999, 999, 999 };
   if (band->RasterIO (GF_Read, col, row, 2, 2, gsep, 2, 2, GDT_Float32, 0, 0) == CE_None) {
      // Compute the slopes used in the interpolation
      const double mx = (lon - adfTransform[0] - col * adfTransform[1]) / adfTransform[1];
      const double my = (lat - adfTransform[3] + row * adfTransform[5]) / adfTransform[5];

      const double ortho = doInterpolate (gsep, mx, my);

      print the lat, lon, ortho height
   }
}

GDALClose (dset);

The doInterpolate is a long standing function used for all our geoid formats we support so I know it is OK. There are no good examples on reading the data from the file just what I can deduce from GDAL.

Is there a way to read these files using the PROJ api?

Thanks for any help you can get me.
@rouault
Copy link
Member

rouault commented Sep 23, 2020

Is there a way to read these files using the PROJ api?

Sure you can perform the bilinar interpolation PROJ does with someling like

echo 2 55 0 | cct proj=vgridshift +grids=uk_os_OSGM15_GB.tif +multiplier=1

Output:

       2.0000        55.0000       42.4159           inf

CC'ing @and-marsh in case he had remarks on how he deduced error of the reprojected grid from the original one

@dfileccia
Copy link
Author

I need something from the PROJ C/C++ api. The string is very useful though. Thanks for the quick reply.

@rouault
Copy link
Member

rouault commented Sep 23, 2020

from the PROJ C/C++ api

just study the source code of cct (src/apps/cct.cpp). It uses proj_create() + proj_trans()

@kbevers
Copy link
Member

kbevers commented Sep 23, 2020

Since I have not found a place to ask questions

See https://proj.org/community/channels.html - the mailing list is the preferred option. As for the rest of your inquiry you seem to be well covered already.

@dfileccia
Copy link
Author

Thanks!

I have my test program now as:

   PJ_CONTEXT* ctx = proj_context_create ();
   int ce = proj_context_errno (ctx);
   if (ce != 0) {
      qDebug () << proj_errno_string (ce);
   }
   else {
      QDir geoidDir (QApplication::applicationDirPath ());
      geoidDir.cd ("../../Geoids");

      QFileInfo fi (geoidDir, "proj/uk_os_OSGM15_GB.tif");

      QByteArray x = QString::asprintf ("+proj=vgridshift +grids=%s +multiplier=1", fi.absoluteFilePath ().toUtf8 ().constData ()).toUtf8 ();
      qDebug () << "DGF: x=" << x;
      PJ* P = proj_create (ctx, x.constData ());
      int pe = proj_errno (P);
      if (pe != 0) {
         qDebug () << proj_errno_string (pe);
      }
      else {
         PJ_COORD src;
         src.lpz.lam = lon * D2R;
         src.lpz.phi = lat * D2R;
         src.lpz.z   = 0;

         PJ_COORD dst = proj_trans (P, PJ_FWD, src);
         pe = proj_errno (P);
         if (pe != 0) {
            qDebug () << proj_errno_string (pe);
         }

         proj_destroy (P);

         return dst.lpz.z;
      }

      proj_context_destroy (ctx);
   }

   return 999;

I get the same results as when I read the file myself so I am left to assume the data is not quite correct or there is still something wrong with the LAT/LON pair. I am assuming WGS84 as the coordinate.

@dfileccia
Copy link
Author

Thanks for info on the mailing list. Would you know if the file is accurate? I downloaded the entire data set from the proj.org site?

@rouault
Copy link
Member

rouault commented Sep 23, 2020

I am assuming WGS84 as the coordinate.

If you look at the metadata of the grid (gdalinfo), you can see this is actually ETRS89 (EPSG:4937)

Thanks for info on the mailing list. Would you know if the file is accurate? I downloaded the entire data set from the proj.org site?

This is from a contributor. We don't guarantee anything. It might eat your children etc. The lineage of the file is indicated at https://github.com/OSGeo/PROJ-data/tree/master/uk_os#united-kingdom-osgm15-height-odn-height---etrs89-ellipsoidal-heights and the conversion script used to generate uk_os_OSGM15_GB.tif at https://github.com/OSGeo/PROJ-data/blob/master/uk_os/OSGM15_GB.sh

@and-marsh
Copy link
Contributor

Hi @dfileccia,
You can look at my PR where we discussed the accuracy of this geoid OSGeo/PROJ-data#25. To get the accuracy and RMSE I used test data from the Ordnance Survey survey developer pack compared to transformation EPSG:4258 -> EPSG:27700+5701 with gdaltransform CLI utility.

Recomened transformation pipeline for EPSG:4258 -> EPSG:27700+5701 is:


unknown id, ETRS89 to ODN height (2) + Inverse of OSGB 1936 to ETRS89 (2) + British National Grid, 0.038 m, UK - Great Britain mainland onshore

PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=uk_os_OSGM15_GB.tif +multiplier=1 +step +inv +proj=hgridshift +grids=uk_os_OSTN15_NTv2_OSGBtoETRS.tif +step +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy

Here geoid height is applied to the coordinates in ETRS89, not WGS84. As I know difference between these CRS can be up to 1m in Europe.

@dfileccia
Copy link
Author

dfileccia commented Sep 24, 2020 via email

@and-marsh
Copy link
Contributor

That proj-string I've posted is output of projinfo utility, try it out, it helps to understand the sequence of transformation and possible variants.

I think I should mention that each geoid file has its own expected CRS of input coordinates. For OSGM15 it's ETRS89, so for results to be correct you should apply geoid to coordinates in ETRS89. By the way, there are two geoids in Australia and they also should be applied not to WGS84, but to GDA97 or GDA2020.

In the transformation above we have next steps:

  1. geoid height correction applied to coordinates in ETRS89
  2. Datum conversion for horizontal coordinates from ETRS89 to OSGB1936 by horizontal grid
  3. applying projection to the coordinates in OSGB1936

I think you may not be interested in last two steps, but you should notice that there are not coordinates in WGS84 at all.

You can't measure accurate enough in WGS84 since it's aligned to the Earth center, when ETRS89 - to the continent. And the continent is moving relative to the center of the Earth all the time. The question is where do you get your source coordinates from? For surveying there are so called NTRIP providers that provide base point with coordinates in the continental CRS. They are different for each country, in Europe it mostly ETRS89.

@dfileccia
Copy link
Author

dfileccia commented Sep 25, 2020 via email

@and-marsh
Copy link
Contributor

@dfileccia, try to compare your results with ones from the Ordnance Survey developer pack. You can skip horizontal coordinates transformation and check only if heights are matching. At least we'll find out if you're using the same transformation as me.

@dfileccia
Copy link
Author

dfileccia commented Sep 25, 2020 via email

@and-marsh
Copy link
Contributor

and-marsh commented Sep 25, 2020 •
edited
Loading

I've tried to convert your test points, look at my results:

>> echo -2.11771100 57.15709900 100 | ./cct +proj=vgridshift +grids=uk_os_OSGM15_GB.tif +multiplier=1 +inv
 -2.1177110000   57.1570990000       50.2294           inf

Close enough to required 50.230. Therefore I suppose the geoid is correct and your WGS84 points are close enough to ETRS89 to don't convert them. So, yes, you're doing something wrong =)

@dfileccia
Copy link
Author

dfileccia commented Sep 25, 2020 via email

@dfileccia
Copy link
Author

dfileccia commented Sep 25, 2020 via email

@rouault rouault closed this as completed Sep 25, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
4 participants
@rouault @dfileccia @kbevers @and-marsh