0

I am reading a shapefile through ArcPy. The shapefile has both geographical coordinate system and projected coordinate system. When I try to get the coordinates of a point geometry in the shapefile, the result is always the point's longitude and latitude.

How could I get its projected coordinates?

The code are as follow:

        sr = arcpy.Describe(shp).spatialReference
        print (sr.name)
        with arcpy.da.SearchCursor(shp, ["SHAPE@"])as lines:
            for line in lines:
                s = line[0].getPart(0)
                for p in s:
                    print(p.X)

the output:

WGS_1984_UTM_Zone_18N
-72.5049105
-72.505729
-72.5058427
...

It looks the shapefile has a projected coordinate system. I am wondering why the points are still recorded with its longitude and latitude.

  • 4
    A shapefile can only have one coordinate system.. are you saying the shapefile could be either geographic or projected? If you're getting geographic coordinates then the shapefile must be defined as geographic.To get your coordinates as a projected coordinate system create a spatial reference (SR = arcpy.SpatialReference(EPSG_CODE)) then the geometry .projectAs(SR).. to be any more helpful I would need to see the code where you're reading the geometries. – Michael Stimson Jan 15 at 6:59
  • Thank you all very much. I have posted the code and output. – S. SHU Jan 15 at 15:48
  • After for line in lines: and before s = line[0].getPart(0) is where you want to line.projectAs(arcpy.SpatialReference(EPSG_Code)) and the part 0 will be returned in the spatial reference given by EPSG_Code, which spatial reference you use depends on where your data is, I would suggest WGS84 UTM might be a good choice but you need to find what zone and hemisphere suits you. – Michael Stimson Jan 15 at 23:15
  • Your original coordinates are obviously in GCS, while the shapefile's coordinate system has been defined as PCS. I'd suggest re-defining the spatial reference as GCS, and then reproject the shapefile to PCS. – Barbarossa Jan 16 at 15:56
  • Thank you all. Barbarossa is right. It has been solved. – S. SHU Jan 17 at 7:28
3

You can use projectAs method of the PointGeometry class:

with arcpy.da.SearchCursor("points_wgs84","SHAPE@") as cursor:
    for row in cursor:
        print row[0].projectAs(arcpy.SpatialReference(3006)).centroid

Output:

667780,048763737 6580256,91702263 NaN NaN
667766,158110956 6580270,80767541 NaN NaN
667796,585255144 6580264,19307885 NaN NaN

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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