plotproj.py

'''
Plot map data in different projections supported by proj.4

Can be called either as

> python plotproj.py proj_strings.txt [out]

or

> python plotproj.py "+proj=lcc +lat_1=20" [out]

The output dir is optional. Uses current dir if not set.

Change PROJ and COASTLINE_DATA to path on your local system
before running the script.
'''

import os
import sys
import json
import subprocess

import numpy as np
import matplotlib.pyplot as plt

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

COASTLINE_DATA = 'data/coastline.geojson'

GRATICULE_WIDTH = 15
N_POINTS = 100

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

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)

def meridian(lon, lat_min, lat_max):
    '''
    Calculate meridian coordinates.
    '''
    coords = np.zeros((N_POINTS, 2))
    coords[:, 0] = lon
    coords[:, 1] = np.linspace(lat_min, lat_max, N_POINTS)
    return coords

def parallel(lat, lon_min, lon_max):
    '''
    Calculate parallell coordinates.
    '''
    coords = np.zeros((N_POINTS, 2))
    coords[:, 0] = np.linspace(lon_min, lon_max, N_POINTS)
    coords[:, 1] = lat
    return coords

def plotproj(proj_str, coastline, frame, graticule, outdir):
    '''
    Plot map.
    '''
    fig = plt.figure()
    axes = fig.add_subplot(111)

    # Plot frame
    for line in frame:
        line = project(line, proj_str)
        x, y = zip(*line)
        plt.plot(x, y, '-k')

    # Plot coastline
    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)

    # Plot graticule
    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(outdir + '/' + proj_str + '.png', dpi=300)

    # Clean up
    fig = None
    del fig
    axes = None
    del axes
    plt.close()

if __name__ == "__main__":
    '''
    '''

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

    if os.path.exists(sys.argv[1]):
        # it must be a file with proj-strings
        with open(sys.argv[1]) as projs:
            projs = projs.readlines()
    else:
        #assumes it is a single roj string
        projs = [sys.argv[1]]

    outdir = '.'
    if len(sys.argv) > 2:
        outdir = sys.argv[2]

    frame = []
    frame.append(parallel(LAT_MIN, LON_MIN, LON_MAX))
    frame.append(parallel(LAT_MAX, LON_MIN, LON_MAX))

    # build graticule
    graticule = []
    for lon in xrange(LON_MIN, LON_MAX+1, GRATICULE_WIDTH):
        graticule.append(meridian(lon, LAT_MIN, LAT_MAX))

    for lat in xrange(LAT_MIN, LAT_MAX+1, GRATICULE_WIDTH):
        graticule.append(parallel(lat, LON_MIN, LON_MAX))

    for i, proj in enumerate(projs):
        proj_str = proj.strip()
        print(i, proj_str)
        plotproj(proj_str, coastline, frame, graticule, outdir)

proj_strings.txt

+proj=aea
+proj=aeqd
+proj=airy
+proj=aitoff
+proj=alsk
+proj=apian
+proj=august
+proj=bacon
+proj=bipc
+proj=boggs
+proj=bonne +lat_1=10
+proj=calcofi
+proj=cass
+proj=cc
+proj=cea
+proj=chamb +lat_1=10 +lon_1=30 +lon_2=40
+proj=collg
+proj=comill
+proj=crast
+proj=denoy
+proj=eck1
+proj=eck2
+proj=eck3
+proj=eck4
+proj=eck5
+proj=eck6
+proj=eqc
+proj=eqdc +lat_1=55 +lat_2=60
+proj=euler +lat_1=67 +lat_2=75
+proj=etmerc
+proj=fahey
+proj=fouc
+proj=fouc_s
+proj=gall
+proj=geos +h=1
+proj=gins8
+proj=gn_sinu +m=2 +n=3
+proj=gnom
+proj=goode
+proj=gs48
+proj=gs50
+proj=hammer
+proj=hatano
+proj=healpix
+proj=rhealpix
+proj=igh
+proj=imw_p +lat_1=34 +lat_2=40
+proj=isea
+proj=kav5
+proj=kav7
+proj=krovak
+proj=labrd
+proj=laea
+proj=lagrng
+proj=larr
+proj=lask
+proj=lcc
+proj=lcca +lat_0=35
+proj=leac
+proj=lee_os
+proj=loxim
+proj=lsat +path=2 +lsat=1
+proj=mbt_s
+proj=mbt_fps
+proj=mbtfpp
+proj=mbtfpq
+proj=mbtfps
+proj=merc
+proj=mil_os
+proj=mill
+proj=misrsom +path=1
+proj=moll
+proj=murd1 +lat_1=10 +lat_2=20
+proj=murd2 +lat_1=10 +lat_2=20
+proj=murd3 +lat_1=10 +lat_2=20
+proj=natearth
+proj=natearth2
+proj=nell
+proj=nell_h
+proj=nicol
+proj=nsper +h=1
+proj=nzmg
+proj=ob_tran +o_proj=latlon +o_lon_p=20 +o_lat_p=20 +lon_0=180
+proj=ocea
+proj=oea +m=1 +n=2
+proj=omerc +lat_1=45 +lat_2=55
+proj=ortel
+proj=ortho
+proj=pconic +lat_1=55 +lat_2=75
+proj=patterson
+proj=poly
+proj=putp1
+proj=putp2
+proj=putp3
+proj=putp3p
+proj=putp4p
+proj=putp5
+proj=putp5p
+proj=putp6
+proj=putp6p
+proj=qua_aut
+proj=qsc
+proj=robin
+proj=rouss
+proj=rpoly
+proj=sinu
+proj=somerc
+proj=stere
+proj=sterea
+proj=gstmerc
+proj=tcc
+proj=tcea
+proj=tissot +lat_1=60 +lat_2=65
+proj=tmerc
+proj=tpeqd +lat_1=60 +lat_2=65
+proj=tpers +h=2
+proj=ups
+proj=urm5
+proj=urmfps +n=0.5
+proj=utm
+proj=vandg
+proj=vandg2
+proj=vandg3
+proj=vandg4
+proj=vitk1 +lat_1=45 +lat_2=55
+proj=wag1
+proj=wag2
+proj=wag3
+proj=wag4
+proj=wag5
+proj=wag6
+proj=wag7
+proj=weren
+proj=wink1
+proj=wink2
+proj=wintri

Output can be seen here: https://dl.dropboxusercontent.com/u/2758559/projections.pdf