I'm trying to convert between two coordinate systems using pyproj:

import pyproj

p1 = pyproj.Proj("+init=EPSG:4326")   
p2 = pyproj.Proj("""PROJCS["NAD27 / California zone II",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982138982]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["standard_parallel_1",39.83333333333334],PARAMETER["standard_parallel_2",38.33333333333334],PARAMETER["latitude_of_origin",37.66666666666666],PARAMETER["central_meridian",-122],PARAMETER["false_easting",2000000],PARAMETER["false_northing",0],UNIT["Foot_US",0.30480060960121924]]""")

p1 creates fine but p2 gives projection not found

I have corresponding points in both coordinate systems but transforming them using EPSG:26742 doesn't give the correct result:

import pyproj

p1 = pyproj.Proj("+init=EPSG:4326")
p2 = pyproj.Proj("+init=EPSG:26742")

a = (-121.330272, 38.547287, 5.661467)
b = (2191840.924977, 321431.770306, 44.012074)

print "Source point", a
print "Target point", b
print 'Converted: ', pyproj.transform(p1, p2, *a) # different :()

3 Answers 3


It would be easier if you use the proj.4 definition:

+proj=lcc +lat_1=39.83333333333334 +lat_2=38.33333333333334 +lat_0=37.66666666666666 +lon_0=-122 +x_0=609601.2192024384 +y_0=0 +datum=NAD27 +units=us-ft +no_defs

or the EPSG code 26742 NAD27 / California zone II

Using cs2cs on your test point with different parameters:

cs2cs +init=epsg:4326 +to +init=epsg:26742 -f "%%.8f" <CAIIdeg.txt >>CAIIft.txt
cs2cs +init=epsg:4326 +to +proj=lcc +lat_1=39.83333333333334 +lat_2=38.33333333333334 +lat_0=37.66666666666666 +lon_0=-122 +x_0=609601.2192024384 +y_0=0 +datum=NAD27 +units=us-ft +vunits=m +no_defs -f "%%.8f" <CAIIdeg.txt >>CAIIft.txt
cs2cs +init=epsg:4326 +to +proj=lcc +lat_1=39.83333333333334 +lat_2=38.33333333333334 +lat_0=37.66666666666666 +lon_0=-122 +x_0=609601.2192024384 +y_0=0 +ellps=clrk66 +units=us-ft +no_defs -f "%%.8f" <CAIIdeg.txt >>CAIIft.txt

I get these results:

2191846.80810027    321450.60775430 18.57432965
2191846.80810027    321450.60775430 5.66146700
2191543.35325840    321414.36473467 18.57432965

The first one uses standard EPSG codes with the internal NADCON datum transformation. The second one uses +vunits=m with the same height as WGS84 in meters. The last line uses the projection without using the NAD27 datum shift, again with the same height in feet. So it seems that proj does not use any vertical datum shift between NAD27 and WGS84.

https://trac.osgeo.org/proj/wiki/VerticalDatums might help to understand the differences.

  • 1
    Thanks, how did you get the proj.4 definition? Commented Sep 23, 2015 at 17:23
  • EPSG 26742 doesn't seem to represent the transform between two point correctly. I've updated the original post with an example. Commented Sep 23, 2015 at 17:25
  • It might depend on the NAD27 to WGS84 transformation used by your data provider. The WKT definition does not reveal that
    – AndreJ
    Commented Sep 23, 2015 at 19:28

I think the fault lies with the PhotoScan software, and not with pyproj or your code. (See this question and this post from the Agisoft forums)

Like Andre, I tried several methods to convert WGS84 values -121.330272, 38.547287 into NAD27 / California Zone 2 coordinates.

WGS84 -> Null Shift -> NAD83 -> NADCON -> NAD27 (EPSG:15851)

WGS84 -> Null Shift -> NAD83(HARN) -> NADCON -> NAD83 -> NADCON -> NAD27 (EPSG:8593)

WGS84 -> Geocentric transform -> NAD27 (EPSG: 1173)

All of these methods produced produced coordinates that were ~2ft away from each other, but ~20ft away from the XY value that came from PhotoScan. I wouldn't trust the coordinates to be very accurate if you need to switch datums until they update the software.

  • According to vdatum.noaa.gov/docs/datums.html#nad27, the NAD27 itself has a built-in error of up to 10 meters. So you can not treat the NADCON grid shift transformation as exact either. It is just the one that is widely used.
    – AndreJ
    Commented Sep 24, 2015 at 3:53
  • That's true. There's nothing really "exact" about NAD27, unless you happen to be standing in Meades Ranch.
    – Mintx
    Commented Sep 24, 2015 at 15:39

A minor addendum to @AndreJ 's answer: You can use the Python GDAL bindings to work with the projections instead of pyproj. For example:

from osgeo import osr

p1 = osr.SpatialReference()

p2 = osr.SpatialReference()
# Note, I had to modify the WKT string slightly to it was explicitly
# PROJECTION["Lambert_Conformal_Conic_2SP"]
p2.ImportFromWkt("""PROJCS["NAD27 / California zone II",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982138982]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",39.83333333333334],PARAMETER["standard_parallel_2",38.33333333333334],PARAMETER["latitude_of_origin",37.66666666666666],PARAMETER["central_meridian",-122],PARAMETER["false_easting",2000000],PARAMETER["false_northing",0],UNIT["Foot_US",0.30480060960121924]]""")

transform = osr.CoordinateTransformation(p1, p2)

transform.TransformPoint(-121.330272, 38.547287, 5.661467)

Which still gives the same answer as @AndreJ got by using the Proj4 string;

(2191543.353258395, 321414.3647346739, 18.57432964916668)

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.