Apologies in advance for the long question, I am fairly confused at this point and I figure if I show my thought process any flaws early on in the piece can be more easily pointed out.
I have a grid of 1000 by 1000 (origin is top left), on which I'm drawing a variety of polygons. I'm trying to convert the x, y locations of each of the vertices of these polygons programmatically to longitude and latitude for use in geojson.
Effectively, I want to map my grid around the globe.
I have run into a few problems, I suspect largely because I'm fairly inexperienced with GIS. I spent some time reading up on what I could, this explanation of map projections and coordinate systems was particularly helpful, and after reading it and other sources, I had the following plan.
With the bounds of -180 degrees to 180 degrees longitude, and 84 to -80 latitude, I should be able to convert arbitrary geographic coordinates to geocentric (x, y) and then scale those values to exist between 0 and 1000, and reverse the process to go the other way.
For the tools to do this, I decided upon pyproj, a wrapper around proj4 as
1) We use python
2) It seems fairly specific to the task at hand which is usually a good sign
From my reading, the bounds of longitude and latitude (in the utm projection, which is what I plan to use) are -180 to 180, and 84 to -80
from pyproj import proj p = Proj(proj='utm', zone='56H', ellps='WGS84') def x_and_y(lon, lat): x, y = p(lon, lat, errcheck=True) return (x, y)
I thought i'd use the following code to work out the x, y values of the lat and long boundaries, and use when scaling from 0-1000 to whatever the min and max x and ys would be.
This is where I run into a brick wall, effectively.
1) Using 0, 0 as lat and long gives me
RuntimeError: latitude or longitude exceeded limits, is this expected?
2) Actually, a lot of values give me errors, i.e., if I do the following:
longs = range(-180, 181, 12) lats = range(84, -81, -4) #print longs #print lats for lon in longs: for lat in lats: print '%s %s' % (lon, lat) print '%s %s' % x_and_y(lon, lat)
I get errors from and including -108 degrees longitude through to and including 60 degrees longitude.
At this point I realised my understanding must be fundamentally flawed and decided to turn to you.
I guess I have the following questions:
1) Am I travelling down the right path here? If my goal is to map (x, y) coordinates onto a global map, is this the way to go?
1a) If I'm not, recommendations?
2) Clearly my understanding is flawed, can anyone see and explain where and why?
3) Are there any tutorials or code samples around for this kind of thing? The pyproj documentation is fairly limited.
Finally, if anyone has recommended reading for this kind of work in general, I'd love a pointer or two on where to start.
Thanks for your time!
PS: If this is better suited for stack overflow, please let me know.