Converting projected coordinates to lat/lon using pyproj fails
If you examine the answer of afalciano in Converting projected coordinates to lat/lon using Python?
1) you define the two projections
# original projection
p = pyproj.Proj("+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-105 +k=90 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs")
# resulting projection, WGS84, long, lat
outProj =pyproj.Proj(init='epsg:4326')
2) conversion
x1,y1 = [-1902530.61073866, 3422503.38926134]
lon,lat = pyproj.transform(p,outProj,x1,y1)
print lon, lat
(104.06918350995736, 53.539892485824495)
3) with many points, simply use a for loop
points = np.array([[-1902530.61073866, -1897767.61073866, -1893004.61073866], [3422503.38926134, 3427266.38926134, 3432029.38926134]])
for x, y in zip(*points):
print pyproj.transform(p,outProj,x,y)
(104.06918350995736, 53.539892485824495)
(103.97445324515407, 53.52377036854651)
(103.87981286777422, 53.507556732352896)
Thanks for the comments. However, I figured out that the problem was related to a missed call of meshgrid. By doing this I was able to run my code and generate the right lon/lat arrays.
x=np.array(dset.variables['x'][:])
y=np.array(dset.variables['y'][:])
xv, yv = np.meshgrid(x, y)
p = pyproj.Proj("+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-105 +k=90
+x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs")
lons, lats = p(xv, yv, inverse=True)
I think it is merely equivalent to what @gene was saying.