I was wondering if the following approach a) would be considered valid and b) someone already did this in python.

I have a python script where users can choose an area of analysis (anywhere on the world). The tool will select data based on this rectangle (data in WGS1984 lat/lon) and process it. The result is a shapefile (from Shapely) that can be mapped in QGIS or Arcgis (etc.). My current problem is: the users of my script don't have any background in GIS. Therefore, the tool should do as much automatically as is possible. As a further automation step, I would like to output the resulting Shapefile in a UTM-Projected Coordinate System suitable for the chosen area.

  • users will always choose a pretty small area (local scale), therefore, it is very unlikely that an area stretches significantly across two UTM-Zones
  • Minimal Distance/Area distortions are important, but only to some degree - the users will perform some calculations that require consistent distances (e.g. GetisOrd Statistic) UTM-Grid UTM-Grid on Wikimedia I think I should be able to select the best UTM-Zone by intersecting with the UTM-Grid:

Is there any python package that will detect the best UTM-Zone given an lat/lng point (and perhaps an area extent)? Are there any concerns regarding my approach?


3 Answers 3


There's the utm package, a bidirectional UTM-WGS84 converter for python. It's able to convert a (latitude, longitude) tuple into an UTM coordinate:

import utm
utm.from_latlon(51.2, 7.5)
>>> (395201.3103811303, 5673135.241182375, 32, 'U')
  • Many thanks Antonio! So I can get the UTM Grid Coordinate "32U" based on a lat/lng pair. How would I select the matching projected Coordinate System (e.g. by getting the SRID))? For example using: project = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'), pyproj.Proj(init='epsg:26913')) # destination coordinate system
    – Alex
    Commented Jan 29, 2018 at 9:59
  • I think I need the SRID to assign a projection system to my shapefile (example: Fiona): with fiona.open(shapeout, 'w', crs=from_epsg(3996), driver='ESRI Shapefile', schema=myschema) as output: [...]
    – Alex
    Commented Jan 29, 2018 at 10:05
  • Found the answer: Here's a (very similar) question and an answer that provides code for extracting the matching EPSG Code (SRID) based on UTM-Coordinates from Lat/Lng
    – Alex
    Commented Jan 29, 2018 at 10:18

Alright, the answer from Antonio above is definitely right and pointed me in the correct direction. Here is the complete code:

# convert_wgs_to_utm function, see https://stackoverflow.com/a/40140326/4556479
def convert_wgs_to_utm(lon: float, lat: float):
    """Based on lat and lng, return best utm epsg-code"""
    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
        return epsg_code
    epsg_code = '327' + utm_band
    return epsg_code
# Get Bounds of user selected area
bound_points_shapely = geometry.MultiPoint(
    [(limLngMin, limLatMin), (limLngMax, limLatMax)])
# get lat/lng from center (below: True centroid - since coords may be multipoint)
input_lon_center = bound_points_shapely.centroid.coords[0][0] 
input_lat_center = bound_points_shapely.centroid.coords[0][1]

# Get best UTM Zone SRID/EPSG Code for center coordinate pair
utm_code = convert_wgs_to_utm(
    input_lon_center, input_lat_center)
# define project_from_to with iteration
# see https://gis.stackexchange.com/a/127432/33092
project = lambda x, y: pyproj.transform(
    pyproj.Proj(init='epsg:{0}'.format(utm_code)), x, y)
# define schema
schema = {
    'geometry': 'Polygon',
# write shapefile
with fiona.open(
    'projected_shapefile.shp', mode='w', encoding='utf-8',
    driver='ESRI Shapefile', schema=schema,
    crs=from_epsg(utm_code)) as c:
    # for each shapely geometry:
    geom_proj = transform(project, geom)
        {'geometry': geometry.mapping(geom_proj)})
  • 1
    Super helpful. Thank you for adding the epsg code calculation.
    – Ian
    Commented Nov 29, 2021 at 17:03

pyproj can assist with this using the database query functionality. This is simpler to use and more robust as it queries the area of use of the CRS, which handles some strangeness in some parts of the world.

Step 1: Convert point to geographic

Step 2: make query


from pyproj import CRS
from pyproj.aoi import AreaOfInterest
from pyproj.database import query_utm_crs_info

utm_crs_list = query_utm_crs_info(
    datum_name="WGS 84",
utm_crs = CRS.from_epsg(utm_crs_list[0].code)

This is what geopandas uses internally for the estimate_utm_crs method.

  • Thanks for the answer. So if it returns a list, you simply use the first epsg entry returned?
    – Alex
    Commented Feb 11, 2022 at 16:17
  • 1
    Yes, that should be sufficient. If there are more than one returned, then you are on a boundary and it is a bit arbitrary which one you pick.
    – snowman2
    Commented Feb 11, 2022 at 17:22

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.