C++ libproj, swap coordinates

Question

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?

Answer 1

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; }

Answer 2

TLDR;
PROJ4 format is not perfect, get JSON projection string with proj_as_projjson and check the key ["target_crs"]["coordinate_system"]["axis"][0]["direction"].

const char* projInfoJson = proj_as_projjson(projContext_m, projEPSG_m, NULL);

json::parse(projInfoJson)["target_crs"]["coordinate_system"]["axis"][0]["direction"] == "north";

---

After a few day debugging and documentation reading, it became clear that the PROJ4 format doesn't provide all the parameters.

That is the reason why cs2cs produced swapped results when using with only EPSGs, or with a projection string.

Also, a few EPSGs (including 2180) for historic/tradition reasons outputs the results in other than 'Easting first' manner.
Have a look at epsg's website and check the 'Coordinate System' part (emphasis mine):

Cartesian 2D CS. Axes: northing, easting (x,y). Orientations: north, east. UoM: m.

Unfortunately, there is no way to tell this information from a PROJ4 type projection string. There are other formats that provide more information.

So the solution is, get the JSON projection string, extract the axis info, and use a swap when the stars align.
Here is the code that does this:

#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);
// ----- here comes the change
    bool swapNeeded = false;
    if (projEPSG_m != nullptr)
    {   
        const char* projInfoJson = proj_as_projjson(projContext_m, projEPSG_m, NULL);
        nlohmann::json projInfo = nlohmann::json::parse(projInfoJson);
        swapNeeded = projInfo["target_crs"]["coordinate_system"]["axis"][0]["direction"] == "north";
    }

    // 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
    if (swapNeeded) {
        std::swap(result.xy.x, result.xy.y);
    }

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

I used the nlohmann::json library for JSON parsing, but one might use other library to read those few characters.