Generating a custom flat 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. And I want to do this using a "flat" cartesian coordinate system originating at a lat/long combination of my choice. I realize that this projection only works for a small area, but that is okay with me, as the area I am looking at is small (<10km in all directions.) Here is what I achieved so far, which is a transverse mercator projection (or so I was told).

I have used the following code.

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'])

So this question is two-fold:

1. Which Projection do I use to achieve a local flat cartesian coordinate system

2. Howto define this projection at my given lat/long combination

I was told to ask this question in a new thread, see my old question. Also a similar question has been asked, but it seemed to be too unspecific to be answered.

Answer 1

You could use an Azimuthal Equidistant projection centred on your origin point.

In PROJ4 syntax, "+proj=aeqd +lat_0=1.234 +lon_0=5.678".

That long/lat coordinate should project to (0,0) in that coordinate system, and points away from there should be in metres. You might want to add "+datum=WGS84 +units=m" but they may be default. Not sure what ellipsoid is default so maybe add an "+ellps=GRS80" or similar...

Note you can construct a projection string using this argument-based form:

>>> pargs = proj.Proj(proj="aeqd", lat_0=gps_lat, lon_0=gps_long, datum="WGS84", units="m")

which means you don't have to make sure your arguments line up with the {0} etc format replacement markers if you try doing it with format.

Answer 2

So @Spacedman pointed me in the right direction, all praise goes to him. I just want to show the complete code for my question:

import pyproj as proj
...
crs_wgs = proj.Proj(init='epsg:4326')  # assuming you're using WGS84 geographic

#Erect own local flat cartesian coordinate system
cust = proj.Proj("+proj=aeqd +lat_0={0} +lon_0={1} +datum=WGS84 +units=m".format(gps_lat_0, gps_long_0))
x, y = proj.transform(crs_wgs, cust, gps_long, gps_lat)

with gps_lat_0/gps_long_0 being the inital coordinates, and gps_lat/gps_long being the list of gps coordinates you want to transfer