Get a cartesian projection accurate around a lat, lng pair

Question

I have list of lat, lng fixes in WGS84, on which I would like to do calculations like distance measurements between point, polygons etc... To do this I was planning on using shapely but then, I need to convert it to a cartesian space, which is only locally accurate.

The problem I have is that my location fixes can come from all over the world, so if I use a fixed projection optimized for my area, I'll introduce errors in other parts of the world. Is it possible to define my own cartesian projection centered around a location pair, depending on the mean position of the current list of locations? The location fixes that I need to do calculations on will always be close to each other, but different lists of location fixes can be spread all over the world.

For example: say I get 5 fixes on which I need to do calculations. Then, I'd like to define a projection which is accurate in the neighbourhood of those fixes, since those fixes wil always be within a few km of each other. When I get the next 5 fixes, which can be on a totally different part of the world, I'd like to define a projection optimized for those locationfixes.

How would I tackle this? It seems like using pyproj (which uses proj.4 if I understood well) is a good idea, but I can't make sense of the string that is required to initialize a projection, as pasted below. Can someone help me?

local_proj = pyproj.Proj(r'+proj=tmerc +lat_0=51.178425 +lon_0=3.561298 +ellps=GRS80 +units=meters')

Answer 1

UPDATE: There is now an excellent little python package utm that should be used rather than the very general snippet below.

Consider using the local UTM zone instead of a cartesian system. The UTM zone works easily via shapely and pyproj and gives you a locally accurate projection system that will work for distance queries:

def convert_wgs_to_utm(lon, lat):
    utm_band = str((math.floor((lon + 180) / 6 ) % 60) + 1)
    if len(utm_band) == 1:
        utm_band = '0'+utm_band
    if lat >= 0:
        epsg_code = '326' + utm_band
    else:
        epsg_code = '327' + utm_band
    return epsg_code

# setup your projections
utm_code = convert_wgs_to_utm(input_lon, input_lat)
crs_wgs = proj.Proj(init='epsg:4326') # assuming you're using WGS84 geographic
crs_utm = proj.Proj(init='epsg:{0}'.format(utm_code))

# then cast your geographic coordinates to the projected system, e.g.
x, y = proj.transform(crs_wgs, crs_utm, input_lon, input_lat)

# proceed with your calculations using shapely...

See this related question on stackoverflow: Determining UTM zone (to convert) from longitude/latitude