The end goal is to update an ongoing shapefile with missing well data. As field agents find missing wells, they can update a spreadsheet with GPS coords (Decimal Degrees). Once a month or so I can run the Python code in ArcMAP and update the shapefile with new points. I am trying to insert the GPS coordinates as points into a shapefile in ArcMAP 10.7 using Python 2.7.
Everything appears to be working correctly, the data table and attributes are correct, the newly added points appear to be correct relative to one another. The issue is when I try to update an existing file. The newly added points are way off somewhere else. In the example here, I tried to add the exact same GPS coords (in decimal degrees) as an existing point. You can see the first entry, the Y_Coord and X_Coord for Well_1 is the same as 1. However the new point is way below the other points.
Immediately, I am thinking this is a coordinate issue of some sort, and searches on the problem indicate a coordinate issue as well but I have a million-times checked to make sure the shapefile, a boundary layer (visual check only), and dataframe are all in the correct geographic and projection coordinate systems. GCS being used is GCS_WGS_1984 (WKID: 4326) and PCS is WGS_1984_UTM_Zone_16N (WKID: 32616).
I have double-checked, triple, a million times checked to make sure the coordinates systems are correct. I have tried to specifically assign it in Python (although I am not sure I am calling it correctly.
Questions:
Am I defining the spatial reference in Python correctly? I feel like I am just calling the object but not actually initiating it.
If not, is there a way to specifically tell Python to apply this in the correct GCS/PCS?
import arcpy
#one coordinate only for testing
coords = [(1, (43.20779, -89.78513))]
# link to shapefile
shapefile = r'M:\ARC\gis_users\Missing_Wells\Missing_Wells.shp'
# Define the spatial reference (GCS_WGS_1984 = 4326, PCS_WGS_1984_UTM16N = 32616)
arcpy.SpatialReference(32616)
#use SHAPE@XY syntax to add point features to a point feature class
with arcpy.da.InsertCursor(shapefile, ['Well_ID', 'Y_Coord', 'X_Coord', 'SHAPE@XY']) as cursor:
#Insert new rows that include the well ID and a x,y coordinate pair
for Well_ID, (lat, lon) in coords:
cursor.insertRow((Well_ID, lat, lon, (lat, lon)))
print('added')
...
sys.exit()
p.s. yes, I know there are tools in ArcGIS Online that can accomplish this goal much easier and in a more timely manner. Our Org will be moving over to ArcGIS Pro next year. We just don't have the budget for it at the moment.
Solved! Go to Solution.
Thanks so much for the help @DavidPike ! This really led me to the right path! Back to fundamentals and quit trying to overspecify!
SOLUTION:
TlDR: Set GCS to WGS84. That's it. No PCS. Then you can enter the data as decimal degrees.
Long version: I was over-specifying when I specified a PCS. As David pointed out, ArcMAP then interpreted the inputs as cartesian coordinates. My SHAPE@XY was then overwriting the coordinates with what I specified it to be so it LOOKED like the coordinates were correct in the attribute table. However, if the point was examined the decimal degree coordinates were actually close to zero. Solution was to just specify the GCS.
Much appreciated again @DavidPike!!
coords = [(1, (43.20779, -89.78513))]
longitude should be X, latitude is Y
Thank you. I know the lat and long is flipped which is why they are flipped in the code as well. I've tried it both ways (lat/long, long/lat) this just happen to be the last time I tried it and just left it. As long as the code is assigning it correctly, which it is.
Your coordinate system is UTM Zone 16. SHAPE@XY will be expecting a tuple of UTM cartesian coordinates.
Can you check the coordinates of your new point in ArcMap - you'll likely see it about 100m from the UTM origin.
I'd also run a search cursor to print out the SHAPE@XY values and you'll be able to see what's going on.
Thanks! This is the first concrete step I've had for a while!
I see what you're saying here. How do I correct this? I tried switching the projection to Web Mercartor (WKID: 3857) but that didn't fix it either. What projection should I be using if I want the input data as decimal degrees? Or should I not be using the SHAPE@XY at all and be using a different tool?
You need to change the projection of your shapefile/feature class to WGS84 (EPSG 4326) then the SHAPE@XY will be a tuple of decimal lat lon.
Run the 'project' tool.
https://pro.arcgis.com/en/pro-app/latest/tool-reference/data-management/project.htm
Obviously this will mean your data is now all in WGS84. There is an alternative where you could run a geometry operation or pass the coords into a few geoprocessing tools as intermediate steps to get the x an y in UTM for you to plug in. But certainly easier to just run with the data in WGS84 if you can get away with it.
Thanks so much for the help @DavidPike ! This really led me to the right path! Back to fundamentals and quit trying to overspecify!
SOLUTION:
TlDR: Set GCS to WGS84. That's it. No PCS. Then you can enter the data as decimal degrees.
Long version: I was over-specifying when I specified a PCS. As David pointed out, ArcMAP then interpreted the inputs as cartesian coordinates. My SHAPE@XY was then overwriting the coordinates with what I specified it to be so it LOOKED like the coordinates were correct in the attribute table. However, if the point was examined the decimal degree coordinates were actually close to zero. Solution was to just specify the GCS.
Much appreciated again @DavidPike!!