Project lat/long with proj.4

project_data.py

'''
Functions for plotting map data in different projections supported by proj.4
'''
import json
import subprocess
from pprint import pprint

import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry.polygon import Polygon
from descartes import PolygonPatch

PROJ = '/Users/kevers/dev/proj.4/src/proj'

GRATICULE_INTERVAL = 10
COASTLINE_DATA = 'data/coastline.geojson'
LAND_DATA = 'data/land.geojson'

PROJ_STR = '+proj=omerc +lonc=135 +lat_1=35 +lat_2=75'
GRATICULE_WIDTH = 20
N_POINTS = 100

LAT_MIN = -90
LAT_MAX = 90
LON_MIN = -180
LON_MAX = 180

lon_overlap = LON_MAX
lat_overlap = LAT_MAX

def project(coordinates, proj_string, in_radians=False):
    '''
    Project geographical coordinates

    Input:
    ------
        coordinates:        numpy ndarray of size (N,2) and type double.
                            longitude, latitude
        proj_string:        Definition of output projection

    Out:
    ----
        numpy ndarray with shape (N,2) with N pairs of 2D
        coordinates (x, y).
    '''

    # proj expects binary input to be in radians
    if not in_radians:
        coordinates = np.deg2rad(coordinates)

    # set up cmd call. -b for binary in/out
    args = [PROJ, '-b']
    args.extend(proj_string.split(' '))

    proc = subprocess.Popen(args, stdin=subprocess.PIPE, stdout=subprocess.PIPE)
    stdout, _ = proc.communicate(coordinates.tobytes())

    out = np.frombuffer(stdout, dtype=np.double)
    return np.reshape(out, coordinates.shape)


with open(COASTLINE_DATA) as data:
    coastline = json.load(data)
'''
with open(LAND_DATA)as data:
    land = json.load(data)
'''

# build graticule
graticule = []
for lon in xrange(LON_MIN, LON_MAX+1, GRATICULE_WIDTH):
    coords = np.zeros((N_POINTS, 2))
    coords[:,0] = lon
    coords[:,1] = np.linspace(LAT_MIN, LAT_MAX, N_POINTS)
    graticule.append(coords)

for lat in xrange(LAT_MIN, LAT_MAX+1, GRATICULE_WIDTH):
    coords = np.zeros((N_POINTS, 2))
    coords[:,0] = np.linspace(LON_MIN, LON_MAX, N_POINTS)
    coords[:,1] = lat
    graticule.append(coords)


# Plot stuff
fig = plt.figure()
axes = fig.add_subplot(111)

for feature in coastline['features']:
    C = np.array(feature['geometry']['coordinates'])

    if np.any(C[:,0] > 180.0):
        C[C[:,0] > 180.0, 0] = 180.0

    if np.any(C[:,0] < -180.0):
        C[C[:,0] < -180.0, 0] = -180.0

    if np.any(C[:,1] > 90.0):
        C[C[:,1] > 90.0, 1] = 90.0

    if np.any(C[:,1] < -90.0):
        C[C[:,1] < -90.0, 1] = -90.0

    coords = project(C, PROJ_STR)
    x, y = zip(*coords)
    plt.plot(x, y, '-k', linewidth=0.5)

'''
for feature in land['features']:
    if len(feature['geometry']['coordinates']):
        continue
    coords = project(feature['geometry']['coordinates'], PROJ_STR)
    poly = Polygon(coords[0])
    patch = PolygonPatch(poly, fc='#cc00cc')
    axes.add_patch(patch)
'''

graticule = project(graticule, PROJ_STR)
for feature in graticule:
    x, y = zip(*feature)
    plt.plot(x, y, '-k', linewidth=0.2)

axes.axis('off')
font = {'family': 'serif',
        'color': 'black',
        'style': 'italic',
        'size': 12}
plt.suptitle(PROJ_STR, fontdict=font)

plt.autoscale(tight=True)
plt.savefig('out.png', dpi=300)