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

Question

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

Answer 1

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.

Answer 2

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.

Answer 3

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)