Adding a custom grid to the proj library

[gabrielwaibel]

Is there a manual how to add custom grid for vertical transformations to the PROJ library.

Currently I am using EGM08 when converting from WGS84 to a local coordinate system. I would like to change it to DHHN2016 (EPSG: 7837). It seems like that this grid is not available. However, a .tif file is available for download here.

Is there an easy way to use this file for the vertical conversion?

[kbevers]

Make sure the grid is discoverable by PROJ, for instance by placing it in the $PROJ_DATA folder, and set up the transformation you need. Here's a very basic, generic example:

echo 12 34 45.67 | PROJ_DATA=/path/to/grid/folder cct +proj=vgridshift +grids=your_grid.tif

Or you can point PROJ towards the specific grid:

echo 12 34 45.67 | cct +proj=vgridshift +grids=./path/to/your_grid.tif

[gabrielwaibel]

Thanks for the fast answer. That works.

But I would also like to do cs2cs transforms and use it in c++ via proj_create_crs_to_crs. How do I do that? Somehow add it here grid_alternatives.sql, and probably a few other places?

[kbevers]

It depends on what you are trying to achieve. It's likely that you need to add custom transformations between existing CRS's as well as modifying the grid alternatives. There's a rather extensive example in data/sql/nkg.sql

[jjimenezshaw]

DHHN2016 grid file is in PROJ-data for at least one year Look at cdn.proj.org

[gabrielwaibel]

Thanks I saw it. I used version 8.2, but I will now use 9.3. It still does not work as with the egm geoid models.

Do I need to insert it into a custom WGS84 to DHHN2016 Transformation into coordinate_operation in proj.db?

[jjimenezshaw]

Thanks I saw it. I used version 8.2, but I will now use 9.3. It still does not work as with the egm geoid models.

Do I need to insert it into a custom WGS84 to DHHN2016 Transformation into coordinate_operation in proj.db?

Could you please write an example of what you are doing and what you expect?

Btw, that type of question is better at the mailing list, not as GitHub issue: https://lists.osgeo.org/mailman/listinfo/proj

[gabrielwaibel]

I will write to the mail list after I approved.

In order to clarify my problem. I have the following code:

// Workspace
#include "conversion/Converter.hpp"

namespace conversion {

Converter::Converter(std::string sourceCrs, std::string targetCrs)
    : sourceCrs_(std::move(sourceCrs)), targetCrs_(std::move(targetCrs)) {

  PJ_INFO info = proj_info();
  std::cout << "PROJ version: " << info.version << std::endl;
  std::cout << "PROJ release date: " << info.release << std::e   ndl;
  std::cout << "PROJ database version: " << info.searchpath << std::endl;

  context_ = proj_context_create();
}

bool Converter::setup() {
  if (!crsCodeToWkt2(sourceCrs_, sourceWkt2_)) {
    std::cout << "Failed to convert source CRS to WKT2" << std::endl;
    return false;
  }

  if (!crsCodeToWkt2(targetCrs_, targetWkt2_)) {
    std::cout << "Failed to convert target CRS to WKT2" << std::endl;
    return false;
  }

  return true;
}

// Function to extract codes (either EPSG or ESRI) and generate WKT2 string
bool Converter::crsCodeToWkt2(const std::string& inputCode, std::string& outputWkt2) {
  if (context_ == nullptr) {
    std::cout << "Failed to create PROJ context" << std::endl;
    return false;
  }

  // This function takes one or multiple EPSG or ESRI codes and returns a compound CRS
  PJ* crs = proj_create(context_, inputCode.c_str());
  if (crs == nullptr) {
    std::cout << "Failed to create compound CRS" << std::endl;
    return false;
  }

  // Convert the compound CRS to WKT2 format
  const char* wkt2 = proj_as_wkt(context_, crs, PJ_WKT2_2019, nullptr);
  if (wkt2 == nullptr) {
    std::cout << "Failed to convert compound CRS to WKT2" << std::endl;
    proj_destroy(crs);
    return false;
  }

  // Store the WKT2 string
  outputWkt2.assign(wkt2);  // Make a copy of the WKT2 string
 std::cout << "WKT2 string: " << outputWkt2 << std::endl;

  // Clean up
  proj_destroy(crs);

  return true;
}

bool Converter::convertToCartesian(double latitude, double longitude, double altitude, double& easting, double& northing,
                                                   double& up) {
  // Create context once at the start
  if (context_ == nullptr) {
    std::cout << "Failed to create PROJ context" << std::endl;
    return false;
  }

  // Create CRS from source WKT2 string
  PJ* crsSource = proj_create(context_, sourceWkt2_.c_str());
  if (crsSource == nullptr) {
    std::cout << "Failed to create CRS from source WKT2" << std::endl;
    return false;
  }

  // Create CRS from target WKT2 string
  PJ* crsTarget = proj_create(context_, targetWkt2_.c_str());
  if (crsTarget == nullptr) {
    std::cout << "Failed to create CRS from target WKT2" << std::endl;
    proj_destroy(crsSource);
    return false;
  }

  // Create CRS to CRS transformation
  std::cout << std::fixed << std::setprecision(std::numeric_limits<double>::digits10) << "Source CRS: " << sourceCrs_
                                             << " Target CRS: " << targetCrs_ << std::endl;
  std::cout << std::fixed << std::setprecision(std::numeric_limits<double>::digits10) << "Latitude: " << latitude
                                             << " Longitude: " << longitude << " Altitude: " << altitude << std::endl;
  // PJ* P = proj_create_crs_to_crs(context_, sourceCrs_.c_str(), targetCrs_.c_str(), /* or EPSG:32632 */
  //                            nullptr);
  PJ* P = proj_create_crs_to_crs_from_pj(context_, crsSource, crsTarget, nullptr, nullptr);
  if (P == nullptr) {
    std::cout << "Failed to create CRS to CRS" << std::endl;
    proj_destroy(crsSource);
    proj_destroy(crsTarget);
    return false;
  }

  /* This will ensure that the order of coordinates for the input CRS */
  /* will be longitude, latitude, whereas EPSG:4326 mandates latitude, */
  /* longitude */
  PJ* norm = proj_normalize_for_visualization(context_, P);
  if (norm == nullptr) {
    std::cout << "Failed to normalize transformation" << std::endl;
    proj_destroy(P);
    proj_destroy(crsSource);
    proj_destroy(crsTarget);
    return false;
  }
  proj_destroy(P);
  P = norm;

  // Prepare input coordinates
  PJ_COORD a = proj_coord(longitude, latitude, altitude, 0);

  /* transform to UTM zone 32, then back to geographical */
  PJ_COORD b = proj_trans(P, PJ_FWD, a);
  easting = b.enu.e;
  northing = b.enu.n;
  up = b.enu.u;
  std::cout << std::fixed << std::setprecision(std::numeric_limits<double>::digits10) << "easting: " << b.enu.e
                                             << ", northing: " << b.enu.n << ", altitude: " << b.enu.u << std::endl;
  // b = proj_trans(P, PJ_INV, b);
  // std::cout << std::fixed << std::setprecision(std::numeric_limits<double>::digits10) << "longitude: " << b.lpz.lam
  //                                            << ", latitude: " << b.lpz.phi << ", altitude: " << b.lpz.z << std::endl;

  /* Clean up */
  proj_destroy(P);
  proj_destroy(crsSource);
  proj_destroy(crsTarget);

  return true;
}

}  // namespace conversion

I am currently using the following coordinate reference systems (CRS) for my data transformations:

• Source CRS: EPSG:4979
• Target CRS: EPSG:31467 + EPSG:3855 (EGM08)

This configuration works perfectly fine. However, I encountered an issue when I switched the target CRS to use the DHHN2016 geoid model:

New Target CRS: EPSG:31467 + EPSG:7837 (DHHN2016) After making this change, the altitude values in my conversions remain unchanged. Specifically, the EPSG:7837 geoid model does not appear to affect the altitude values during the transformation.

Question:

What could be causing the altitude values to remain unchanged with the DHHN2016 model? I suspect that I may not be using the correct geoid model for this transformation. Any guidance on how to resolve this issue would be greatly appreciated.

[jjimenezshaw]

Since 9.4.0 this works:

echo 50 9 0 | PROJ_NETWORK=ON  ./cs2cs EPSG:4979 EPSG:31467+7837
5540407.24  3500074.96 -47.87