4

I'm trying to convert between topocentric coordinates and geographic coordinates with pyproj.

First, I made a topocentric CRS using the proj4 string from the example here. Then, I made a geodetic CRS using EPSG:4326. The Transformer gave me issues when I tried to make the conversion. Code here:

from pyproj import CRS
from pyproj import Transformer

crs_4326 = CRS.from_epsg(4326)
crs_topo = CRS.from_proj4(f"+proj=topocentric +ellps=WGS84 +lon_0=5 +lat_0=55 +h_0=200")
Transformer.from_crs(crs_4326, crs_topo, always_xy=True)

I get this error: ProjError: Input is not a transformation.: (Internal Proj Error: pipeline: Pipeline: Mismatched units between step 1 and 2)

How can I use a Transformer to convert to or from a topocentric system?

1

2 Answers 2

2

I came across this question trying to solve the same problem in C. Here's what I came up with as a solution to my problem, which may be similar to yours.

#include <stdio.h>
#include <stdlib.h>

#include <proj.h>

int main(void)
{
    PJ_CONTEXT *C;
    PJ *P;
    PJ_COORD a, b;

    C = proj_context_create();
    
    P = proj_create(C, "proj=pipeline step proj=cart ellps=WGS84 step proj=topocentric ellps=WGS84 lon_0=-74 lat_0=40 h_0=0");
    a = proj_coord(-74.0, 40.1, 0, 2020);

    a.lpzt.lam = proj_torad(a.lpzt.lam);
    a.lpzt.phi = proj_torad(a.lpzt.phi);

    b = proj_trans(P, PJ_FWD, a);
    
    printf("east %fmeters, north %fmeters, height %fmeters\n", b.enu.e, b.enu.n, b.enu.u);
    
    proj_destroy(P);
    proj_context_destroy(C);
}

Compiles with gcc test.c -lproj, outputs: east 0.000000meters, north 11103.553727meters, height -9.689709meters. Note the lack of + in the proj string, and the conversion from degrees to radians.

I gleaned this by creating a debug build of proj 9.0.0, and debugging the cct program. The manual gives this as a command line example:

echo 2.12955 53.80939444 73 2020 | cct -d 3 +proj=pipeline +step +proj=cart +ellps=WGS84 +step +proj=topocentric +ellps=WGS84 +lon_0=5 +lat_0=55 +h_0=200

and in debugging this I found the proj string created from the command line parameters, and the conversion from degrees to radians. Et voila!

The lat/lon I've used in my example C code are, roughly New York City. The coordinate a lat/lon are exactly 0.1 degrees north of the topocentric origin I put in the proj string, and reassuringly the northing that is output is +11.1km (which is about right).

2
  • Ahhh thanks a ton, I did not realize the C/C++ runtime interface doesn't like the preceding '+' for args like it does at command line.
    – mr3
    Commented Apr 21, 2023 at 19:47
  • @mr3 no worries. I remember going nuts trying to work this out, the documentation left me floundering, confused and dazed. Definitely one worth putting up here!
    – bazza
    Commented Apr 22, 2023 at 9:08
1

This might work for you (from @snowman2 's link):

from pyproj import CRS
from pyproj import Transformer


t = Transformer.from_pipeline("+proj=pipeline +step +proj=cart +ellps=WGS84 +step +proj=topocentric +ellps=WGS84 +lon_0=5 +lat_0=55 +h_0=200")
# Check a point close to the specified origin
t.transform(5.0001, 55.0001, 201)
# (6.399598230444809, 11.132708495215015, 0.9999870808796674)

The two-step pipeline here is also described in this section at the bottom of your link: https://proj.org/operations/conversions/topocentric.html#topocentric-origin-described-as-geographic-coordinates

I'm using pyproj==3.3.1, proj==8.0.1.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.