[GIS] Mapping x, y pixel coordinates to latitude and longitudes on the WGS84 datum

coordinate systemgeojsonprojpython

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.

Best Answer

I did not know geojson. But from the perspective of the coordinate systems this does not fit together:

1.

longs = range(-180, 181, ...)
lats = range(84, -81, ...)

2.

Proj(proj='utm', zone='56H', ellps='WGS84')

1: you use a geographic coordinat system. Valid values are

x -180 to + 180

y -90 to +90

You use WGS. WGS is (the most important part of) a geographic coordinate system (GCS). A geographic coordinate system uses a three-dimensional spherical surface to define locations on the earth. A point is referenced by its longitude and latitude values.

2: you use a projected coordinate system. Valid values ​​are e.g.: x 123456 (or 56123456 with zone identifyer) y 55123456 A projected coordinate system (PCS) is defined on a flat, two-dimensional surface. A map projection converts features between a spherical surface and a flat surface. The UTM coordinate sytem use WGS as its spherical surface. But the coordinates on the spherical surface are projected to the flat surface. If you work with UTM you work with coordinates on the flat surface.

Related Question