topic UTM Conversion in Python in Python Questions

UTM Conversion in Python

JohanL — Fri, 11 Nov 2011 13:13:37 GMT

Hi Everyone,

I was wondering if anyone could give me any ideas/thoughts/links to get myself familiar UTM Conversion from lat long. I am using a lot of DEM data (i.e. SRTM or AGDEM) who provide the data in tif format with the name representing the lat long in WGS84 (center?).

i.e. ASTMGTM2_N37W110 = North037 and W-110 degrees that according to http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html conversion tool would be in UTM zone 12.

If I could apply an automated process similar to that conversion tool on the downloaded datasets (using the name in lat/long)then I would simply need to apply the standard project raster tool in ArcGIS to apply the actual projection of the data. I realize that one could define the projection in the project raster manually and then use a batch process but I was wondering if a fully automated solution in python is possible?

Sincerely,
Bjorn

Re: UTM Conversion in Python

JakeSkinner — Sat, 11 Dec 2021 17:38:00 GMT

Hi Bjorn,

You could automate the procedure of finding which zone each raster falls in. Each install of ArcGIS includes reference data in the directory 'C:\Program Files (x86)\ArcGIS\Desktop10.0\Reference Systems'. There is a UTM.shp located here and you can use this to determine each zone for each TIFF. Here is an example on how to do this:

import arcpy
from arcpy import env
env.overwriteOutput = True
prjFile = r"C:\Program Files (x86)\ArcGIS\Desktop10.0\Coordinate Systems\Geographic Coordinate Systems\World\WGS 1984.prj"
spatRef = arcpy.SpatialReference(prjFile)

# list each raster and extract coordinates
env.workspace = r"C:\temp\python\rasters"
lstRasters = arcpy.ListRasters("*", "TIF")
for ras in lstRasters:
    rasName = ras.split('.')[0]
    lat_lon = rasName.split('_')[1]
    if 'N' in lat_lon:
        lat = int(lat_lon[1:3])
    else:
        lat = int(lat_lon[1:3]) * -1
    if 'W' in lat_lon:
        lon = int(lat_lon[4:]) * -1
    else:
        lon = int(lat_lon[4:])
    
    # Create a point based on the coordinate values
    point = arcpy.Point(lon, lat)

    UTMfc = r"C:\TEMP\Python\UTM.shp"
    
    # create an empty point feature class with a field called 'Name'
    env.workspace = r"C:\temp\python\test.gdb"    
    arcpy.CreateFeatureclass_management(env.workspace, "LatLon", "POINT", "", "", "", spatRef)
    arcpy.AddField_management("LatLon", "Name", "TEXT")

    # Populate the feature class with the coordinate value and the raster name
    rows = arcpy.InsertCursor("LatLon")
    row = rows.newRow()
    row.Name = rasName
    row.Shape = point
    rows.insertRow(row)

    del row, rows

    # Perform an intersect to find which zone the raster falls in
    arcpy.Intersect_analysis(["LatLon", UTMfc], "Intersect", "", "", "POINT")
    rows = arcpy.SearchCursor("Intersect")
    for row in rows:
        print str(row.Name) + " is in zone UTM" + str(row.Zone)

    arcpy.Delete_management("LatLon")
    arcpy.Delete_management("Intersect")

You could probably append some additional code onto this to project the raster based on the zone it falls within.

Re: UTM Conversion in Python

ArnaudTemme — Wed, 08 Aug 2012 16:59:30 GMT

Hello colleagues,

Thank you for this code, it's partly what I was looking for too. Large amounts of ASTER to project to an appropriate UTM zone. I can get the zone using your sample code above, but that is not all, it seems.

To create the out_coordinate_system for use in arcpy.ProjectRaster_management , a huge string with the coordinate system's properties needs to be made. Is there a way to do that, based on the now-known zone?

Re: UTM Conversion in Python

KimOllivier — Sat, 11 Dec 2021 17:38:03 GMT

Yes, don't bother with the string. Create a SpatialReference object and set the Factory Code or the Name.

sr = arcpy.SpatialReference()
sr.factoryCode = 32759 # case sensitive
sr.create() # essential step that is non-intuitive!
print sr.name # to check

At 10.1 the library of strings for standard projections has disappeared! They are now built in to ArcGIS.
The factorycode is the most reliable but hardest to find, look up the projected_coordinate_systems.PDF in the root ArcGIS documentation directory. I can never remember what name property is required to be set to make it work.

Re: UTM Conversion in Python

ArnaudTemme — Tue, 14 Aug 2012 21:34:53 GMT

Thanks a lot, this is brilliant! I will be able to make it work now.