0

I have a set of coordinates in epsg:3857:

webmerc_coords = (-8956562.6 5368787.8)

Is there a way to determine what the appropriate local NAD83 UTM projection should be for this point?

I'm specifically looking for a Python-based (GeoPandas/pyproj etc.) or even better, JavaScript-based (proj4js) implementation of a function that can do this.

The coords given above, if plugged into this hypothetical function, would return 'EPSG:26917' as that's the NAD83 grid those EPSG:3857 coords most appropriately 'belong to'.

Looking at https://epsg.io/26917 for example, there are center coordinates listed, as well as bounds for this projection. I think what I'm looking for is essentially a function that will find the projection whose center coordinates are closest to MY coordinates or am I wrong in my logic somewhere here and asking for something impossible?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
  • 1
    There are many coordinate transformation libraries like Proj, Proj4js, geotools, and others. Do some preliminary studies and come back with more focused question. – user30184 Feb 10 '22 at 13:21
  • 1: https://gis.stackexchange.com/a/247881, 2: https://gis.stackexchange.com/a/93275, 3: https://gis.stackexchange.com/a/390565 – Kadir Şahbaz Feb 10 '22 at 13:33
  • You are missing what my question actually is. I know how to reproject coordinates from projection A to projection B using any of the above mentioned libraries. I'm not asking how to do that. What I'm after/wondering about is if there is a library/built in function that can find the most relevant "projection B" given projection A, where projection A is epsg:3857 and projection B is some NAD83 Zone xx. – user11705556 Feb 10 '22 at 14:00
  • Something like this https://gis.stackexchange.com/questions/269518/auto-select-suitable-utm-zone-based-on-grid-intersection? – user30184 Feb 10 '22 at 14:14
  • That's looking promsing @user30184 ! I'll try and implement and report back. Thanks! – user11705556 Feb 10 '22 at 15:26

2 Answers2

2

Solution was here: Auto-select suitable UTM Zone based on Grid Intersection

Essentially all I needed was the function:

def convert_wgs_to_utm(lon: float, lat: float):
    """Based on lat and lng, return best utm epsg-code"""
    utm_band = str((np.floor((lon + 180) / 6 ) % 60) + 1)
    print(utm_band)
    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

To get that to NAD83, you never have negative latitudes so dump the if statement, and replace the '327' (for WGS) with '269':

def convert_wgs_to_utmNAD83(lon: float, lat: float):
    """Based on lat and lng, return best utm epsg-code"""
    utm_band = str((np.floor((lon + 180) / 6 ) % 60) + 1)
    print(utm_band)
    if len(utm_band) == 1:
        utm_band = '0'+utm_band
    epsg_code = '269' + utm_band
    return epsg_code

Thanks, User30184 for pointing me to that post :-)

2

pyproj can assist with this using the database query functionality.

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="NAD83", 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)

https://pyproj4.github.io/pyproj/stable/api/database.html#pyproj.database.query_crs_info

snowman2
  • 7,321
  • 12
  • 29
  • 54
  • Awesome that's another great sol'n. Thanks – user11705556 Feb 11 '22 at 11:59
  • I definitely recommend it as there can be some wierdness in some regions of the world when determining the utm zone. This handles it because it queries the areas of use for the CRS. – snowman2 Feb 11 '22 at 14:36