python - Auto-select suitable UTM Zone based on Grid Intersection

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

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/questions/40132542/get-a-cartesian-projection-accurate-around-a-lat-lng-pair
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
#Get Bounds of user selected area
bound_points_shapely = geometry.MultiPoint([(limLngMin, limLatMin), (limLngMax, limLatMax)])
get lat/lng from center
input_lon_center = bound_points_shapely.centroid.coords[0][0] #True centroid (coords may be multipoint)
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/questions/127427/transforming-shapely-polygon-and-multipolygon-objects

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',
}
with fiona.open('projectedShapefile.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)})