How to create a projection from a CRS string using pyproj?

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 :()
``````

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.

• 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 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. 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. 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()
p1.ImportFromEPSG(4326)

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)
``````