[GIS] Convert shapefile to WKT using ArcPy

arcgis-10.4arcpycoordinate systemwell-known-text

I am attempting to convert a shapefile to WKT using ArcGIS 10.4. I have seen the question/answer posted here, Converting geometry to WKT using ArcPy?, which is helpful. However, when I follow the code and write it to the text file I create, the spatial reference is incorrect. The geometric shape is identical to the original shapefile, but the WKT maps to a different continent. I don't have experience writing to WKT, so I could be making an amateurish mistake that is causing my spatial reference to be off. 

import arcpy

# Create new empty text file
outFile = open(r"my_path.txt", "w")

# My shapefile
fc = r"my_path.shp"

# Convert to WKT and write to file
cursor = arcpy.da.SearchCursor(fc, ["SHAPEFILE@WKT"])
for row in cursor:
    print(row[0])
    outFile.write((row[0]) + "\n")
del row
del cursor
outFile.close()

The projection that I want my WKT to be in is NAD_1983_UTM_Zone_13N.

Best Answer

Firstly, "SHAPEFILE@WKT" should be "SHAPE@WKT".

Secondly, your WKT will be in exactly the same coordinate system as your shapefile. If you're getting strange results, your data is probably not in NAD 1983 UTM Zone 13N, but some other spatial reference and is being projected on the fly when you view it in ArcGIS.

To specify the output spatial reference for your search cursor, pass in a SpatialReference object.

For example:

import arcpy

# My shapefile (which is in WGS84 Geographic)
fc = r"C:/Temp/test.shp"

#I want to output WKT in NAD 1983 UTM Zone 13N
sr = arcpy.SpatialReference('NAD 1983 UTM Zone 13N')
# or could use
# sr = arcpy.SpatialReference(26913)

#Print out and compare the spatial refs
print(arcpy.Describe(fc).SpatialReference.exportToString())
print(sr.exportToString())

# Convert to WKT and print
cursor = arcpy.da.SearchCursor(fc, ["SHAPE@WKT"])
for row in cursor:
    print(row[0])
    break # just print 1st row to demonstrate

# Convert to WKT projected on the fly and print
cursor = arcpy.da.SearchCursor(fc, ["SHAPE@WKT"], spatial_reference=sr)
for row in cursor:
    print(row[0])
    break # just print 1st row to demonstrate

del row
del cursor

Output:

GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 11258999068426.2;-100000 10000;-100000 10000;8.98315284119522E-09;0.001;0.001;IsHighPrecision
PROJCS['NAD_1983_UTM_Zone_13N',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',-105.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]];-5120900 -9998100 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
MULTIPOLYGON (((-104.88510543531396 43.417242462433592, -104.48059706432281 43.667085868045774, -104.81372160513905 44.190567289328442, -105.15874345098445 43.893134663599653, -105.68222487226711 43.738469698220683, -104.88510543531396 43.417242462433592)))
MULTIPOLYGON (((509301.25108411006 4807156.9673131397, 541874.8360413569 4835028.9378186846, 514886.82021686854 4893055.9525310742, 487250.06276868912 4860015.9002533574, 445063.47957695409 4843052.0422064001, 509301.25108411006 4807156.9673131397)))
Related Question