Coordinate transformation in arcpy

16743
4
Jump to solution
05-07-2015 02:32 AM
ChrisHills
Occasional Contributor

As part of a script I need to convert coordinates from one SRS to another. This is the code I have so far, but I am missing the bit in the middle. I have done some extensive googling but I was not able to find a solution so I would appreciate some assistance.

import arcpy

inputSRS = 'Projected Coordinate Systems/National Grids/Europe/British National Grid' # Projected
outputSRS = 'Geographic Coordinate Systems/World/WGS 1984' # Geographic
trsOut = 'OSGB_1936_To_WGS_1984_Petroleum'

srIn = arcpy.SpatialReference(inputSRS)
srOut = arcpy.SpatialReference(outputSRS)

inX = 210000
inY = 310000

# Code goes here

print outX, outY # "-4.81042685256 52.6543521582"
0 Kudos
1 Solution

Accepted Solutions
JaiSiva1
New Contributor III

Hi Chris,

  try the following code :

import arcpy
inputSRS = 'Projected Coordinate Systems/National Grids/Europe/British National Grid' # Projected  
outputSRS = 'Geographic Coordinate Systems/World/WGS 1984' # Geographic  
trsOut = 'OSGB_1936_To_WGS_1984_Petroleum'  
srIn = arcpy.SpatialReference(inputSRS)    
srOut = arcpy.SpatialReference(outputSRS)  
pt = arcpy.Point()
ptGeoms = []
pt.X = 210000  
pt.Y = 310000  
ptGeoms.append(arcpy.PointGeometry(pt,srIn ))
arcpy.CopyFeatures_management(ptGeoms, 'base')
arcpy.Project_management('base',trsOut, srOut)

It will first create a Feature named 'base' with input SRS and then creates your output feature with output SRS.

You can set the environment to save the output in a specific location.

View solution in original post

4 Replies
JaiSiva1
New Contributor III

Hi Chris,

  try the following code :

import arcpy
inputSRS = 'Projected Coordinate Systems/National Grids/Europe/British National Grid' # Projected  
outputSRS = 'Geographic Coordinate Systems/World/WGS 1984' # Geographic  
trsOut = 'OSGB_1936_To_WGS_1984_Petroleum'  
srIn = arcpy.SpatialReference(inputSRS)    
srOut = arcpy.SpatialReference(outputSRS)  
pt = arcpy.Point()
ptGeoms = []
pt.X = 210000  
pt.Y = 310000  
ptGeoms.append(arcpy.PointGeometry(pt,srIn ))
arcpy.CopyFeatures_management(ptGeoms, 'base')
arcpy.Project_management('base',trsOut, srOut)

It will first create a Feature named 'base' with input SRS and then creates your output feature with output SRS.

You can set the environment to save the output in a specific location.

JaiSiva1
New Contributor III

additionally you can cross check by :

with arcpy.da.SearchCursor("SR2",['SHAPE@X','SHAPE@Y']) as cursor:
    for row in cursor:
        print row[0]
        print row[1]

It will print :

  -4.810426852

  52.654352152

ChrisHills
Occasional Contributor

Jai, thank you very much this is exactly what I wanted!

Edit: Is it possible to do this without saving to a file? I will be doing this on thousands of coordinates so this is undesirable as it will place increased load on the disk, and slow down the process.

0 Kudos
curtvprice
MVP Esteemed Contributor

Is it possible to do this without saving to a file?

Yes, you can use the arcpy shape method .projectAs(). This can be tricky but should work great!

projectAs (spatial_reference, {transformation_name})

Projects a geometry and optionally applies a geotransformation.

To project, the geometry needs to have a spatial reference, and not have an UnknownCoordinateSystem. The new spatial reference system passed to the method defines the output coordinate system. If either spatial reference is unknown the coordinates will not be changed. The Z- and measure values are not changed by the ProjectAs method.

Point—Help | ArcGIS for Desktop

PointGeometry—Help | ArcGIS for Desktop

I also highly recommend using WKID numbers instead of those long coordinate system names which are easy to misspell. A very easy place to look them up is http://spatialreference.org.

Note this is untested code:

import arcpy  
inputSRS = arcpy.SpatialReference(7405)  # British National Grid
outputSRS = arcpy.SpatialReference(4326) # GCS WGS84
gt = 'OSGB_1936_To_WGS_1984_Petroleum'    
pt = arcpy.Point()  
pt.X = 210000    
pt.Y = 310000    
print "Input XY: {} {}".format(pt.X, pt.Y)
ptgeo = arcpy.PointGeometry(pt, inputSRS)
ptgeo1 = ptgeo.projectAs(outputSRS, gt)
pt1 = ptgeo1.lastPoint
print "Output XY: {} {}".format(pt1.X, pt1.Y)