3

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

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? – nickponline Sep 23 '15 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. – nickponline Sep 23 '15 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 Sep 23 '15 at 19:28
1

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 Sep 24 '15 at 3:53
  • That's true. There's nothing really "exact" about NAD27, unless you happen to be standing in Meades Ranch. – Mintx Sep 24 '15 at 15:39
0

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)

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.