Converting projected coordinates to lat/lon using pyproj fails

Question

I have been trying to work on this issue on myself mainly referring to this post (Converting projected coordinates to lat/lon using Python?) but still didn't succeed.

I have some 2-D data function of X, Y projected coordinates (units meters) but no lat/lon information.

Given that I have already the Proj4 string in the file I'm trying to use PyProj to convert these coordinate into lat/lon with a simple script.

import pyproj
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")
lon, lat = p(x, y, inverse=True)

Where x and y are the variables taken directly from the file (they are 1-D with 1121 and 881 elements, respectively). Unfortunately I'm getting a dimensional error given that x and y have not the same dimension (Buffer lengths not the same).

So my question is... Why should they have the same size? I thought that you could have any rectangular subset into a projection. Am I missing something?

• mkennedy: What's +k=90 doing? Try taking that out.

Answer 1

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)

• Joseph Sheedy: Iterating over the array is slow if the array is large. transform() accepts numpy arrays so transformed_points = pyproj.transform(p, outProj, *points) is faster.

Answer 2

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.

Answer 3

The issue may lie in the fact of the format of your input:

Input coordinates can be given as python arrays, lists/tuples, scalars or numpy/Numeric/numarray arrays. Optimized for objects that support the Python buffer protocol (regular python and numpy array objects).

Python interface to PROJ.4 library

Performs cartographic transformations and geodetic computations. The Proj class can convert from geographic (longitude,latitude) to native map projection (x,y) coordinates and vice versa, or from one map projection coordinate system directly to another.

The Geod class can perform forward and inverse geodetic, or Great Circle, computations. The forward computation involves determining latitude, longitude and back azimuth of a terminus point given the latitude and longitude of an initial point, plus azimuth and distance. The inverse computation involves determining the forward and back azimuths and distance given the latitudes and longitudes of an initial and terminus point.

git repository where you may access the most up-to-date source

• Droid: But they are normal arrays. For example x is array([-1902530.61073866, -1897767.61073866, -1893004.61073866, ..., 3422503.38926134, 3427266.38926134, 3432029.38926134]) .