Generating custom projection in pyproj

Question

I have little experience with pyproj, but I want to achieve the following task:

I have a trajectory of a robot, given in lat/long coordinates which I want to convert to x/y Coordinates. In aerospace sciences, we used the NED coordinate system for this, therefore I want to transfer a bunch of lat/long coordinates to local x/y coordinates whereby the first x/y coordinate shall be at 0/0,(which is typically just the start of the trajectory).

So I have been trying to use the PyProj library, which sucessfully converts lat/long to x/y, however the coordinate system in which the result is given is that of england, however ideally I want the 0/0-point to be my first lat/long point. So I though perhaps it would be possible to define a custom coordinate system originating at a specific coordinate system, but couldn't find anything on google.

Therefore: does anyone know how I can define a NED coordinate System at a specific lat/long combination?

Here is my code so far:

import pyproj as proj

crs_wgs = proj.Proj(init='epsg:4326')  # assuming you're using WGS84 geographic
    crs_bng = proj.Proj(init='epsg:27700')  # use a locally appropriate projected CRS
    # then cast your geographic coordinate pair to the projected system
    x, y = proj.transform(crs_wgs, crs_bng, s['gps_lng'], s['gps_lat'])

Answer 1

You can do this but I'm not sure its a good thing. I'd work in EPSG:27700 because all my other map data is most likely to have this, and if I lose the projection info for something with a custom CRS then guessing it is very very hard.

But here's how to do it:

Suppose we start at 3W, 58N. What's that in EPSG:27700?

>>> proj.transform(crs_wgs, crs_bng, -3,58)
(340989.6785735455, 901635.9303086236)

The underlying description of EPSG:27700 can be looked up at SpatialReference https://spatialreference.org/ref/epsg/27700/proj4/ and is:

>>> bng = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs"
>>> 

Note the x_0 and y_0 offset parameters. You need to adjust these offsets by your start coordinate so that your start coordinate projects to (0,0). I'm rounding it to integers here but for precision use the full decimal value:

>>> -100000-901635
-1001635
>>> 400000-340989
59011
>>> 

Now create a projection string like EPSG:27700 but with those new offsets:

>>> cust = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=59011 +y_0=-1001635 +ellps=airy +datum=OSGB36 +units=m +no_defs"

and a corresponding pyproj object:

>>> crs_cust = proj.Proj(cust) 

If we got this right then 3W 58N should project to (0,0):

>>> proj.transform(crs_wgs, crs_cust, -3,58)
(0.6785735454977839, 0.930308623588644)

which is close enough (these are in metres) and I think the remaining fraction of a metre is due to my integer rounding. Try it with the full decimal and you should get (0,0), or very very close to it.