1

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')
KarelV
  • 487
  • 1
  • 5
  • 20

1 Answers1

7

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

songololo
  • 4,724
  • 5
  • 35
  • 49
  • Thanks, that is indeed what I needed! – KarelV Oct 20 '16 at 07:05
  • 1
    Warning, there are irregular zones in UTM: https://gis.stackexchange.com/questions/2561/what-was-the-rationale-for-the-non-standard-utm-zones-near-norway – bugmenot123 Jan 29 '18 at 15:06
  • Thanks for this conversion. However, this does not work in Python 2.7 for zones less than 10. convert_wgs_to_utm(-158, 71.8) gives 3264.0 when it should give 32604. Changing the first line to utm_band = str(int(math.floor((lon + 180) / 6 ) % 60) + 1).zfill(2) and then removing the first if statement results in the correct answer. – user20408 Oct 15 '18 at 15:35
  • There is now an excellent little package called `utm` that handles these conversions, though I don't know if it supports Python 2. – songololo Oct 15 '18 at 17:03