Auto-select suitable UTM Zone based on Grid Intersection

Question

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 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?

Answer 1

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

Answer 2

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:4326'),
    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)
    c.write(
        {'geometry': geometry.mapping(geom_proj)})

Answer 3

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

https://pyproj4.github.io/pyproj/stable/examples.html#find-utm-crs-by-latitude-and-longitude

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",
    area_of_interest=AreaOfInterest(
        west_lon_degree=-93.581543,
        south_lat_degree=42.032974,
        east_lon_degree=-93.581543,
        north_lat_degree=42.032974,
    ),
)
utm_crs = CRS.from_epsg(utm_crs_list[0].code)

This is what geopandas uses internally for the estimate_utm_crs method.