3

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).

enter image description here

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.

7
  • 2
    Note: your CSV values are not longitude, latitude, they are eastings and northings, you want longitude and latitude as output.
    – user2856
    Jan 29, 2022 at 22:25
  • 5
    I also wouldn't mix arcpy and ogr. If you want to stick with your arcpy.da.InsertCursor, build a arcpy.PointGeometry(arcpy.Point(X, Y), arcpy.SpatialReference(2274)) and use the "SHAPE@" token in the InsertCursor, each PointGeometry will be reprojected automatically to match the FeatureClass which is in EPSG:4326.
    – user2856
    Jan 29, 2022 at 23:02
  • 4
    Because putting data into into an ogr object only to get it back out again and put it into an ArcGIS object is wasteful and will likely slow your code down. It also means you need to concentrate on learning two very different libraries, not ideal when you are new. Don't get me wrong, I use GDAL/OGR (and fiona, rasterio, geopandas etc...) as well as arcpy. I just don't recommend mixing them. You can do the entire operation in ogr or fiona or geopandas or arcpy, there's no need to mix and match.
    – user2856
    Jan 31, 2022 at 12:30
  • 1
    If you want to maintain the general processing flow you can still get rid of gdal-- create PointGeometry() objects from your CSV then use the point.ProjectAs() method: pro.arcgis.com/en/pro-app/2.8/arcpy/classes/pointgeometry.htm
    – mikewatt
    Jan 31, 2022 at 22:28
  • 1
    You should edit the question title. It says reprojecting with arcpy, while you are reprojecting with ogr, and it seems want to stick with it. The arcpy should be removed from the tags as well.
    – 0kcats
    Feb 4, 2022 at 9:12

3 Answers 3

10

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:

enter image description here

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
1
  • This could work, but I'm wondering if I can do this by editing my current script. I have to do some field calculations and can't simply just do an append. Jan 31, 2022 at 15:15
1

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.

-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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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