Wednesday 10 April 2019

python - Converting projected coordinates to lat/lon using pyproj fails


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?



Answer



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.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...