tile_gall_peters_projection.py

#!/usr/bin/env python

import math
import json
import pyproj

R = 6371000


def tile_gall_peters_projection():


    src_proj = pyproj.Proj('+proj=cea +lon_0=0 +lat_ts=45 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs')
    dst_proj = pyproj.Proj(init='EPSG:4326')

    
    x = lambda lng: (R * math.pi * lng) / (180 * math.sqrt(2))
    y = lambda lat: R * math.sqrt(2) * math.sin(lat)

    xmin, xmax = x(-180), x(180)
    ymin, ymax = y(-90), y(90)

    nx, ny = 20, 20

    x_range = xmax - xmin
    y_range = ymax - ymin
    
    x_size = x_range / nx
    y_size = y_range / ny

    gj = {
        'type': 'FeatureCollection',
        'features': [
        ]
    }

    for j in range(ny):
        for i in range(nx):

            x0 = xmin + i * x_size
            x1 = xmin + (i + 1) * x_size

            y0 = ymin + j * y_size
            y1 = ymin + (j + 1) * y_size

            s0, t0 = pyproj.transform(src_proj, dst_proj, x0, y0)
            s1, t1 = pyproj.transform(src_proj, dst_proj, x1, y1)

            gj['features'].append({
                'type': 'Feature',
                'properties': {},
                'geometry': {
                    'type': 'Polygon',
                    'coordinates': [
                        [
                            [s0, t0],
                            [s0, t1],
                            [s1, t1],
                            [s1, t0],
                            [s0, t0]
                        ]
                    ]
                }
            })
            # print s0, t0, s1, t1
    
    print json.dumps(gj, indent=2)


if __name__ == '__main__':
    tile_gall_peters_projection()