4

I'm trying to convert US census data (NAD83) to WGS84. I know the difference between the two is quite minor, but trying to apply the transformation using pyproj does not do anything. e.g.:

from pyproj import Proj, transform
p1 = Proj(init='epsg:4269')
p2 = Proj(init='epsg:4326')
transform(p1, p2 , -80.00001, 40.00001)

yields:

(-80.00001, 40.00001000000001)

Any idea why? Shouldn't I expect small shifts of ~1 meter?

Zohar Bar-Yehuda
  • 141
  • 1
  • 2
  • 2
    It will depend on what pyproj is using (or not) for the transformation between the 2 systems. It might be using 0,0,0 in which case--no offsets. I think you might have to create the Proj object with a proj.4 string that includes +towgs84 and parameters rather than using the WKIDs. – mkennedy Nov 29 '18 at 00:45

2 Answers2

3

Until you project your point(s), you won't know where they lie in one system in relation to themselves in the other system.

WGS84 and NAD83 are not, themselves projections, but coordinate systems differing by the datum (or 'representations of the spheroidal earth, defined mathematically') that they use.

Imagine two points which have almost no horizontally spacing, but one very high up in the air and one close to the ground. In this case you have no shift in meters because you have no flat surface, or single plane, to travel across in order to measure. Your output set of points is essentially completely coincident with the first set so yes, you should have no practical shift, but you must project each set to calculate the difference.

There is a very good explanation here in the answer to Difference Between Projection and Datum that I think addresses your question further.

I hope this was helpful.

Pictory
  • 194
  • 10
  • 1
    I see what you mean. So in order to "transform" between the two, I can either:
    1. Go through a metric projection (e.g. NAD -> UTM -> WGS)
    2. Use a transformation matrix such as those referenced here (WGS_1984_(ITRF08)_To_NAD_1983_2011).
     – Zohar Bar-Yehuda Nov 28 '18 at 16:34
  • I think the second approach will serve you better. Your reference looks like exactly what you're after. You may have also seen this post in your searches, more to your initial question regarding pyproj. – Pictory Nov 28 '18 at 18:27
1

As suggested in @mkennedy's comment pyproj is using no offsets (towgs84=0,0,0) in the transformation between NAD83 and WGS84.

p1 = Proj(init='epsg:4269')
print(p1.definition)

'proj=longlat datum=NAD83 no_defs ellps=GRS80 towgs84=0,0,0'

As also suggested in @mkennedy's comment, you might have to create the Proj object with a proj.4 string that includes +towgs84 and parameters rather than using the WKIDs.

From @mkennedy's answer to the question "PROJ4 / PostGIS transformations between WGS84 and NAD83 transformations in Alaska" I derived the following:

from pyproj import Proj, transform
p2 = Proj(init='epsg:4326')

no offsets


p1 = Proj(init='epsg:4269')


method: coordinate frame


p1cf = Proj(proj='latlong', ellps='GRS80', datum='NAD83', towgs84='-0.9956,1.9013,0.5215,-0.025915,-0.009426,-0.0011599,-0.00062')


method: position vector


p1pv = Proj(proj='latlong', ellps='GRS80', datum='NAD83', towgs84='-0.9956,1.9013,0.5215,0.025915,0.009426,0.0011599,-0.00062')


print(transform(p1, p2 , -80.00001, 40.00001))
print(transform(p1cf, p2 , -80.00001, 40.00001))
print(transform(p1pv, p2 , -80.00001, 40.00001))

Output:

(-80.00001, 40.00001)
(-80.00001904501643, 40.00003299028925)
(-80.00001618601023, 40.00001788467985)
user2856
  • 65,736
  • 6
  • 115
  • 196