3

I'm trying to convert pairs of x,y coordinates to latitude/longitude using python. The projection is "NAD_1983_HARN_StatePlane_Oregon_North_FIPS_3601_Feet_Intl" and a sample x,y pair is 7625659, 714714. So far, I have found this question and have tried using the code in the first solution:

from pyproj import Proj, transform

inProj = Proj(init='epsg:????')
outProj = Proj(init='epsg:????')
x1,y1 = 7625659,714714
x2,y2 = transform(inProj,outProj,x1,y1)
print x2,y2

However, I'm having trouble finding an epsg code for the location. I came across this but I can't seem to find the epsg code for that area. All I can see is the sr-org code. Also, it seems strange that the location on the map is not in Oregon. I was able to convert it (the output looked reasonable) using the tool at http://beta.ngs.noaa.gov/gtkweb/ using datum=NAD83, projection for conversion=SPC, Norting=714714, Easting=7625659, Units=International feet and SPC Zone=OR N-3601. How should I convert this projection to latitude and longitude?

Edit: I tried using the tool at http://prj2epsg.org/search which gave me 2913 (I also tried 3646 per BradHards suggestion, it came up the same) for inProj and 4326 (BradHards) for outProj, but it's giving me a location (-62.0472510133 32.7243167795) which is in the middle of the ocean, not in Portland. The lat/lon from http://beta.ngs.noaa.gov/gtkweb/ is (45.605018016, -122.752592697) which is in Portland (Sorry, I can't comment on the site...).

2
  • Have a look at this, spatialreference.org/ref/epsg/4152
    – Fezter
    Nov 16, 2016 at 3:02
  • The outProj is probably meant to be epsg:4326 if you want WGS-84 longitude / latitude. It might be something else if you want a different ellipsoid. Perhaps epsg:3646 for the input.
    – BradHards
    Nov 16, 2016 at 3:08

1 Answer 1

4

I needed to add preserve_units = True and the epsg codes provided by BradHards.

from pyproj import Proj, transform

inProj = Proj(init='epsg:3646', preserve_units = True)
outProj = Proj(init='epsg:4326')
x1,y1 = 7625659,714714
x2,y2 = transform(inProj,outProj,x1,y1)
print x2,y2

This returns lat/lon (-122.752652268, 45.6050168495). Strangely enough http://beta.ngs.noaa.gov/gtkweb/ gives me lon/lat (-122.752592697, 45.605018016) which is a couple of feet over. I'm going to assume that's a round off error somewhere and it shouldn't make too much of a difference.

1

Your Answer

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

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