SpatialReference in Loop

[EricGreen]

Hey all,

I am trying to set the projection of a temp shapefile so I can add a length field.  I am using the factory code to set it.  I have a huge list of shapefiles in (GCS) and their corresponding ideal factory codes for State Plane CS.  I am trying to loop through each and set the projection, but it only works the first time.  The second time I get this error:

Traceback (most recent call last):
  File "C:\Users\egreen\Desktop\BigFiles KEEP\CreatePROX_FINAL.py", line 67, in <module>
    sr = arcpy.SpatialReference(ProjCode)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\mixins.py", line 927, in __init__
    self._arc_object.createFromFile(item)
RuntimeError: ERROR 999999: Error executing function.
the input is not a geographic or projected coordinate system

The error is misleading since it is not a bad factory code (or WKID).  The second shapefile has the same factory code.  I think the error is due to the variable sr already existing.  I'm curious if I need to destroy sr somehow (I don't see a method for that though).  Or am I misunderstanding this.

Here's my code:

import arcpy
import xml.etree.cElementTree as ET
import os
import xlrd
from arcpy import env
import os.path
import shutil
import zipfile
import glob

#set paths to documents
myPath = os.getcwd() + "\\ScriptFiles"
excel_path = myPath + "\\Inputs.xls"
wb = xlrd.open_workbook(excel_path)

# Load XLRD Excel Reader
sheetname = wb.sheet_names() #Read for XCL Sheet names
sh1 = wb.sheet_by_index(0) #Login

#flag for first row
Skip = True

### Set Geoprocessing environments
arcpy.env.scratchWorkspace = myPath + "\\Scratch.gdb"
arcpy.env.workspace = myPath + "\\PROXCreate.gdb"

Coordinate_System = "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]]" # provide a default value if unspecified
# Set overwrite option
arcpy.env.overwriteOutput = True

tmpRoads = "tmpRoads"
tmpBoundaries = "tmpBoundaries"

#Loop through Excel
for rownum in range(sh1.nrows):

    rows = sh1.row_values(rownum)
    
    #skip first row
    if Skip == True:
        Skip=False
    else:
        CountyName = rows[0]
        StateName = rows[2]
        InputRoadShpfl = rows[10]
        InputBoundaryShpfl = rows[11]
        ProjName = rows[9]
        ProjCode = rows[12]

        if arcpy.Exists(tmpRoads):
            arcpy.Delete_management(tmpRoads)
            
        if arcpy.Exists(tmpBoundaries):
            arcpy.Delete_management(tmpBoundaries)

        if (arcpy.Exists(myPath + "\\TigerFiles\\" + InputRoadShpfl)) and (myPath + "\\TigerFiles\\" + InputBoundaryShpfl):

            print "Executing shapefiles:" + InputBoundaryShpfl + " & " + InputRoadShpfl + "(" + StateName + "\\" + CountyName + ")"
           
            # Make a layer from the feature class
            #arcpy.MakeFeatureLayer_management(myPath + "\\TigerFiles\\" + InputRoadShpfl, tmpRoads)
            arcpy.CopyFeatures_management(myPath + "\\TigerFiles\\" + InputRoadShpfl, tmpRoads)
            arcpy.MakeFeatureLayer_management(myPath + "\\TigerFiles\\" + InputBoundaryShpfl, tmpBoundaries)

            # Process: Define Projection (3)
            #sr = arcpy.SpatialReference(102629)
            sr = arcpy.SpatialReference(ProjCode)
            #sr = arcpy.SpatialReference(ProjName)
            #sr.factoryCode = ProjCode
            sr.create()
            arcpy.DefineProjection_management(tmpRoads, sr)

[NobbirAhmed]

Add a print statement to check type of ProjCode as follows:

ProjCode = rows[12]
print type(ProjCode)

Most probably the type is not an "int" - in that case convert it to an int:

ProjCode = rows[12]
ProjCode = int(ProjCode.strip())

[EricGreen]

Add a print statement to check type of ProjCode as follows:

ProjCode = rows[12]
print type(ProjCode)

Most probably the type is not an "int" - in that case convert it to an int:

ProjCode = rows[12]
ProjCode = int(ProjCode.strip())

Interesting.  I tried what you suggested and ProjCode is defined as a float.  Which is odd again because it does run the first time and not the second (and they are both floats).  I tried converting to int but it crashes (because it is a float not a string).  So I removed the strip (ProjCode = int(ProjCode)) and it runs but I look at the results and it is clearly using some other projection.  When I run it without the int conversion the result of the first loop is using the correct coordinate system.  It seems the int conversion must be forcing the spatialreference to some generic or default system.

I still feel like the issue is somewhere in the sr object, but I can;t be sure.

Confused.

[NobbirAhmed]

I don't get the error with only this code:

#Loop through Excel
for rownum in range(sh1.nrows):

    columns = sh1.row_values(rownum)
    print "row num = ", rownum
    
    #skip first row
    if Skip == True:
        Skip=False
        
    else:
        ProjCode = columns[0]
        print ProjCode
        print type(ProjCode)  # check type for all rows
        ProjCode = int(ProjCode)

        sr = arcpy.SpatialReference(ProjCode)
        sr.create()
        print sr.name

Could you just try the same?

I'm also confused on this line (may be not related to your problem):

if (arcpy.Exists(myPath + "\\TigerFiles\\" + InputRoadShpfl)) and (myPath + "\\TigerFiles\\" + InputBoundaryShpfl):

Here, Exists checks only the first shapefile but not the second one.

[EricGreen]

First, I really appreciate your help!

I too don't get a crash with the int code.  But my results indicated that the wrong coordinate system was being used.  Although, I know that is not the case after printing the name of the spatial reference as you suggested.  Seems my problem is elsewhere.

The rest of my code simply adds a field and then calculates the length of each shape.  Unfortunately the numbers I am getting are really small.  E.g. I get 0.00000033935 for a shape that I measure to be 0.03384 miles (or 0.00000312066 for 0.296182 miles).  Not sure what unit that is in, but most of my shapes (roads) should be about a mile long.  This number makes me think it is in hemispheres!  They seem to be about 100,000 miles.  Oddly, if I run my code without the int and for only one record... it works!

if (arcpy.Exists(myPath + "\\TigerFiles\\" + InputRoadShpfl)) and arcpy.Exists(myPath + "\\TigerFiles\\" + InputBoundaryShpfl):

            print "Executing shapefiles:" + InputBoundaryShpfl + " & " + InputRoadShpfl + "(" + StateName + "\\" + CountyName + ")"
           
            # Make a layer from the feature class
            #arcpy.MakeFeatureLayer_management(myPath + "\\TigerFiles\\" + InputRoadShpfl, tmpRoads)
            arcpy.CopyFeatures_management(myPath + "\\TigerFiles\\" + InputRoadShpfl, tmpRoads)
            arcpy.MakeFeatureLayer_management(myPath + "\\TigerFiles\\" + InputBoundaryShpfl, tmpBoundaries)

            # Process: Define Projection (3)
            #sr = arcpy.SpatialReference(102629)
            print ProjCode
            sr = arcpy.SpatialReference(ProjCode)
            #sr = arcpy.SpatialReference(ProjName)
            #sr.factoryCode = ProjCode
            sr.create()
            print sr.name
            arcpy.DefineProjection_management(tmpRoads, sr)
 
            # Process: Add Field
            arcpy.AddField_management(tmpRoads, "Length", "DOUBLE", "6", "3", "", "", "NULLABLE", "NON_REQUIRED", "")

            # Process: Calculate Field
            arcpy.CalculateField_management(tmpRoads, "Length", "!shape.length@miles!", "PYTHON", "")

            # Process: Define Projection (2)
            arcpy.DefineProjection_management(tmpRoads, "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]]")

[NobbirAhmed]

One more thing. What's the coordinate systems of your data "before" you are applying Define Projection? If not "Unknown" then try using Project tool instead. Define projection does not change the coordinate system - it just updates the metadata.

And, if you can share your data then I can try to reproduce your issue on my machine and see what's going on.

[EricGreen]

Well I'll be!  The project tool does the trick!

arcpy.Project_management(tmpRoadsGCS, tmpRoads, sr)

Thanks a million!

BTW the input files I am using are straight from TigerCensus.