c++coordinate-transformationproj

C++ libproj, swap coordinates


I have a C++ program running on linux for years. It uses some coordinate transformation, we use the libproj (version 6.3.1) for this purpose.

Interestingly, one client gave the input parameter epsg:2180 then everything has blown up.

It seems like that for a few epsg transformations the output is swapped. Here are the commands I used to prove it:

First, get the projection string:

$ projinfo epsg:2180

PROJ.4 string:
+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs

If I use it with a sample coodinate (LatLon format) like below, I got the result as Easting first. (Basically the command is cs2cs epsg:4326 +to <projection string here> -f %.4f)

echo "55.12578267 16.77291767" | cs2cs epsg:4326 +to +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs -f %.4f

Result that everyone expects:

358037.8085 809219.1725 0.0000

However, if I use the simple command:

echo "55.12578267 16.77291767" | cs2cs epsg:4326 epsg:2180 -f %.4f

Then the result is flipped (customer not happy):

809219.1725 358037.8085 0.0000

Problem doesn't stop here, in C++, I use the same library, with the same problem. The question is, how can I make sure that the Easting will be the first parameter?

I've already tried to use proj_normalize_for_visualization, as the documentation suggested that it would do a swap on the coordinates, but in reality it doesn't swap the result.

#include <cstdlib>
#include <cstdio>
#include <proj.h>
int main(int argc, char** argv)
{
    PJ_CONTEXT* projContext_m = proj_context_create();
    
    // hardcode epsg (in reality it is user defined)
    int epsg_p = 2180;
    
    char _buff[100];
    snprintf(_buff, sizeof (_buff), "EPSG:%d", epsg_p);
    PJ* projEPSG_m = proj_create_crs_to_crs(projContext_m,
                                            "EPSG:4326",
                                            _buff, NULL);

    if (projEPSG_m != nullptr)
    {   
        proj_normalize_for_visualization(projContext_m, projEPSG_m);
    }
    
    // fake input (in reality it is user defined)
    PJ_COORD input_p;
    input_p.xy.x = 55.12578267;
    input_p.xy.y = 16.77291767;
    
    PJ_COORD result = proj_trans(projEPSG_m, PJ_FWD, input_p);    // input_p is the coordinate we need to transform

    // we assume `result.x` is the Easting
    printf("Easting: %f, Northing %f", result.xy.x, result.xy.y);
    
    return 0;
}

Output:

Easting: 809219.172518, Northing 358037.808510

Expected:

Easting: 358037.808510, Northing 809219.172518

What should I do to enforce result.x is always an Easting?


Solution

  • I had the same problem and I solved it with proj_normalize_for_visualization().

      auto p1 = proj_create(ctx, name.c_str()); // create crs
      if (!p1) throw Err() << "Can't create libproj crs: " << name;
    
      // always north up
      auto p2 = proj_normalize_for_visualization(ctx, p1);
      if (p2) { proj_destroy(p1); p1 = p2; }
    
      // promote to 3d
      auto p3 = proj_crs_promote_to_3D(ctx, nullptr, p1);
      if (p3) { proj_destroy(p1); p1 = p3; }