Reprojecting list of latitude/longitude values in .csv to points using ArcPy

Question

I am creating a tool using Python to convert .csv files into points. However, the points were collected in EPSG: 2274 and not WGS84 (EPSG:4326), so in order to create the points I have to convert the lat/long values into WGS84 for them to plot correctly I believe. I have created a script that does this and it works, but the points are not placed correctly geographically. They end up all clustered and in Arkansas (highlighted points) instead of Tennessee (blue pins).

Is there something I need to change in the code to project these values into points?

import csv
import arcpy
import ogr, osr

#set input variables
spreadsheet = "C:\\Monument_Surveyors\\ControlPointInformation_Test_oldconfiguration.csv"
fc = r"C:\Monument_Surveyors\Monuments\Monuments_copying.gdb\Monuments_STAGE_editing_SPCS"

# Spatial Reference System
inputEPSG = 4326
outputEPSG = 2274

#open csv file
with open(spreadsheet, "r") as monument_sheet:
    #create reader
    csvReader = csv.reader(monument_sheet)
    #call header
    header = next(csvReader)

#grab index position by using field name
#not necessary if index is already known
    latIndex = header.index("Latitude:")
    longIndex = header.index("Longitude")

#create empty list for string coordinates
    coords = []

#parse through rows and call data in each row index
    for row in csvReader:
        lat = float(row[latIndex])
        long = float(row[longIndex])
        coords.append([lat, long])
        
print("Coords: " + str(coords))

# Spatial Reference System
inputEPSG = 2274
outputEPSG = 4326

#convert string coordinates to float coordinates
coordsFloat = []
for coord in coords:
    for i in coord:
        coordsFloat.append(float(i))
xypoints = [coordsFloat[x:x+2] for x in range(0, len(coordsFloat),2)]

projCoords = []
projXcoords = []
projYcoords = []
latList = []
longList = []
i = 0
while i < len(xypoints):
    #define lat and long index within each coordinate pair
    pointY = xypoints[i][0]
    pointX = -abs(xypoints[i][1])
    
    # create a geometry from coordinates
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(pointX, pointY)
            
    # create coordinate transformation
    inSpatialRef = osr.SpatialReference()
    inSpatialRef.ImportFromEPSG(inputEPSG)
            
    outSpatialRef = osr.SpatialReference()
    outSpatialRef.ImportFromEPSG(outputEPSG)
            
    coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)
        
    # transform point
    point.Transform(coordTransform)
    
    #define lat and long
    lat = point.GetY()
    long = point.GetX()
    latList.append(lat)
    longList.append(long)
    #print("long: " + str(long))
    #print("GetX: " + str(point.GetX()))
    
    #concatonate list element-wise
    projCoords = [None]*(len(latList)+len(longList))
    projCoords[::2] = latList
    projCoords[1::2] = longList
    #print("latList: " + str(latList))
    #print("longList: " + str(longList))
    
    i+=1

#group coordinates
projCoordsTuple = zip(projCoords[0::2], projCoords[1::2])
    
#convert tuple to a list within a list
projCoordsList = []
for row in projCoordsTuple:
    latIndex = row[0]
    longIndex = row[1]
    projCoordsList.append([longIndex, latIndex])
print(projCoordsList)   


try:
    #create cursor
    cur = arcpy.da.InsertCursor(fc, ("SHAPE@X", "SHAPE@Y"))
    #iterate through rows to insert point
    for i in projCoordsList:
        cur.insertRow(i)
    
    print ("feature created")

except:
    print(arcpy.GetMessages())
    print ("feature not created")

finally:
    del cur

EDIT: I would like to update this specific script as a simple Append will not do the trick - I need to do some additional field calculations that do not exist in the .csv spreadsheet.

Answer 1

You could do this in a few lines using arcpy.management.XYTableToPoint (or arcpy.management.MakeXYEventLayer) then arcpy.management.Project. Much simpler. The initial feature class (in EPSG:2274) can be an in-memory FC to speed up processing.

in_csv = "test.csv"
out_fc = "test.shp"  # or "path/to/file.gdb/somefc"
x_field, y_field = "X", "Y"

arcpy.management.XYTableToPoint(in_csv , "memory/testcsv", x_field, y_field, coordinate_system=arcpy.SpatialReference(2274))
arcpy.management.Project("memory/testcsv", out_fc , arcpy.SpatialReference(4326))

Some demo points in EPSG:2274:

X,Y,id,
1718965.95131462,590521.042027273,1
1817376.99676817,592536.415227254,2
1807030.13644896,517453.796786448,3
1722344.26104266,456262.307004273,4

Result:

And if you need the InsertCursor because you are appending to an existing FeatureClass with existing features:

in_csv = "test.csv"
out_fc = "test.shp"  # or "path/to/file.gdb/somefc"
x_field, y_field = "X", "Y"

arcpy.management.XYTableToPoint(in_csv , "memory/testcsv", x_field, y_field, coordinate_system=arcpy.SpatialReference(2274))
arcpy.management.Project("memory/testcsv", "memory/testfc" , arcpy.SpatialReference(4326))
arcpy.management.Append("memory/testfc", out_fc)  # may need to specify schema_type="NO_TEST" and field_mapping

Answer 2

If your points are showing up in the wrong place, have you verified that your OGR coordinate transformation is working correctly?

coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

Since you are doing the coordinate transformation outside of ArcGIS, it is very unlikely that ArcGIS is responsible for the location being wrong in the map you are viewing.

Using OGR to do the coordinate transformation is unnecessary since you are already using ArcPy and have access to ArcGIS code to do the transformation. Why bring more Python packages into the mix if you don't need to.

Instead of attempting to modify your script, especially since I don't think using OGR for coordinate transformation is necessary, I will provide a solution using ArcPy insert cursor that you can adapt to your situation.

import arcpy
import csv

# set input variables
spreadsheet =       # path to CSV file containing coordinates with "X,Y,ID" header
fc = #              # path to pre-existing feature class in WGS 84 spatial reference

# EPSG codes for CSV coordinates and feature class
epsg_in = 2274
epsg_out = 4326

inSpatialRef = arcpy.SpatialReference(epsg_in)
outSpatialRef = arcpy.SpatialReference(epsg_out)
datum_transformation = arcpy.ListTransformations(inSpatialRef, outSpatialRef)[0]

with open(spreadsheet, newline='') as file:
    with arcpy.da.InsertCursor(fc, "SHAPE@", datum_transformation) as cur:
        csvfile = csv.reader(file)
        next(csvfile)
        for x, y, oid in csvfile:
            pointGeom = arcpy.PointGeometry(arcpy.Point(float(x), float(y)), inSpatialRef)
            cur.insertRow([pointGeom])
 
 

The InsertCursor can do the datum transformation for you, so there is no need to do it ahead of time. Although I created a PointGeometry first, it may be the cursor could be defined with ["SHAPE@X", "SHAPE@Y"] instead "SHAPE@" and just have the coordinates passed in separately.

Answer 3

Finally found a method that works...

import csv
import arcpy

#set input variables
spreadsheet = "C:\\Monument_Surveyors\\ControlPointInformation_Test_oldconfiguration.csv"
fc = r"C:\Monument_Surveyors\Monuments\Monuments_copying.gdb\Monuments_STAGE_editing_SPCS"

# Spatial Reference System
inputEPSG = 2274
outputEPSG = 4326

#open csv file
with open(spreadsheet, "r") as monument_sheet:
    #create reader
    csvReader = csv.reader(monument_sheet)
    #call header
    header = next(csvReader)

#grab index position by using field name
#not necessary if index is already known
    latIndex = header.index("Latitude:")
    longIndex = header.index("Longitude")
    elevIndex = header.index("Elev.:")

#create empty list for string coordinates
    coords = []

#parse through rows and call data in each row index
    for row in csvReader:
        lat = float(row[latIndex])
        long = float(row[longIndex])
        elev = float(row[elevIndex])
        coords.append([lat, long, elev])

#convert string coordinates to float coordinates
coordsFloat = []
for coord in coords:
    for i in coord:
        coordsFloat.append(float(i))
xypoints = [coordsFloat[x:x+3] for x in range(0, len(coordsFloat),3)]
print(xypoints)

try:
    #create cursor
    cur = arcpy.da.InsertCursor (fc, ["SHAPE@Y", "SHAPE@X", "SHAPE@Z"])
    for i in xypoints:
        cur.insertRow(i)
    print("feature created")

except Exception as e:
    print(e)
    print(arcpy.GetMessages())
    print ("feature not created")

finally:
    del cur