0

Calculate UTM Zone (Cartography) assigns a UTM zone to each feature within a feature class as a field value. The value returned is a 'UTM spatial reference string', which is a long string that looks something like this:

PROJCS["GCS North American 1983 UTM Zone 10S (Calculated)",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-123.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

I'd like to use this information to create a spatial reference object in the projection indicated in the UTM spatial reference string. I'm not sure how to go about this, however. It looks like parsing out a value in between the brackets after PROJCS would get me close. However, given the different format of the string (spaces vs underscores, etc.) I'm not sure how reliably I'd be able to get it to work anyways. Is there a better way to get a spatial reference object that each feature falls within?

Maybe there's a better way to get UTM zone per feature, but I would like the code to stand alone and not require a UTM Zone feature class to function.

  • 1
    The CRS WKT should work fine, but it doesn't match the corresponding NAD83 UTM zone in the software. It's even calculating the latitude band ("S") instead of the hemisphere. If your data's all in the northern hemisphere, you could extract the zone number, then concatenate it with 269xx, so 26910 is NAD83 UTM 10 North. – mkennedy Aug 26 '15 at 22:40
2

Riffing on @mkennedy's comment above, here's a small script to parse and create a spatial reference if the datum is NAD83.

import arcpy

inputString = r'PROJCS["GCS North American 1983 UTM Zone 10S (Calculated)",GEOGCS[ETCETCETC]'
editString = inputString.replace(" ", "_").upper()

if "NORTH_AMERICAN_1983" in editString and "UTM_ZONE" in editString:
  utmIndex = editString.find("UTM_ZONE_")
  utmZone = editString[(utmIndex + 9):(utmIndex + 11)]
  if utmZone[-1:].isalpha:
    utmZone = "2690" + utmZone[:1]
  else:
    utmZone = "269" + utmZone
  sr = arcpy.SpatialReference(int(utmZone))
  print sr.name

NAD_1983_UTM_Zone_10N

  • Unfortunately, Calculate UTM zone doesn't always return a value with NORTH AMERICAN 1983. Run on a different data set than the one in my op: PROJCS["GCS WGS 1984 UTM Zone 10S (Calculated)",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],'... This is why the code in my answer hunts down three-letter words and checks if the first two are numeric for its zone matching. – Emil Brundage Aug 27 '15 at 19:05
1

I've come up with a solution that's maybe not the prettiest but it works. Below is code that will parse out the needed information from a feature class field after the Calculate UTM Zone tool is run, and then use it to create the spatial reference object.

>>> cursor = arcpy.da.SearchCursor ("test", "utm_zone")
>>> for utmStr, in cursor:
    words = utmStr.split (" ")
    for word in words:
        if len (word) != 3:
            continue
        zone = word[:2].lstrip ("0")
        if not zone.isdigit():
            continue
        grid = word[-1]
        if grid in "XWVUTSRQPN":
            hem = "N"
        else:
            hem = "S"
        srName = "WGS 1984 UTM Zone {0}{1}".format (zone, hem)
        sr = arcpy.SpatialReference (srName)
        print sr.name
        break

WGS_1984_UTM_Zone_7N
WGS_1984_UTM_Zone_8N
WGS_1984_UTM_Zone_9N
WGS_1984_UTM_Zone_10N
WGS_1984_UTM_Zone_11N
WGS_1984_UTM_Zone_12N
WGS_1984_UTM_Zone_13N
WGS_1984_UTM_Zone_14N
>>>

If a more concise solution is available I'd like to see it.

0

The arcpy SpatialReference object has the loadFromString method, which allows for the updating of a spatial reference object from WKT.

>>> import arcpy
>>> sr = arcpy.SpatialReference ()
>>> utmWkt = 'PROJCS["GCS North American 1983 UTM Zone 10S (Calculated)",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-123.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]'
>>> sr.loadFromString (utmWkt)
>>> sr.name
u'GCS North American 1983 UTM Zone 10S (Calculated)'

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.