make_scatter

#!/usr/bin/env python
import matplotlib.pyplot as plt
import numpy as np
import psycopg2
from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Polygon

def plot_poly(lon1,lat1,lon2,lat2,out_path):
    def draw_screen_poly( lats, lons, m):
        x, y = m( lons, lats )
        xy = zip(x,y)
        poly = Polygon( xy, facecolor='red', alpha=0.4 )
        plt.gca().add_patch(poly)

    lats = [ lat1, lat2, lat2, lat1 ]
    lons = [ lon1, lon1, lon2, lon2 ]

    plt.clf()
    m = Basemap(projection='npstere',boundinglat=73,lon_0=0,resolution='l')
    m.drawcoastlines()
    m.drawmapboundary()
    draw_screen_poly( lats, lons, m )

    plt.savefig(out_path, dpi=150, bbox_inches='tight', pad_inches=0)
    #plt.show()

def plot_data(x,y,out_path,title,icor):
    plt.clf()
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    ## the data
    N=1000
    #x = np.random.randn(N)
    #y = np.random.randn(N)

    ## left panel
    ax1.scatter(x,y,color='blue',s=3,edgecolor='none')
    #ax1.set_aspect(1./ax1.get_data_ratio()) # make axes square
    plt.axis((0, 700, 0, 700))
    # Line
    ax1.plot([0,700],[0,700],color='black', alpha=0.8, linewidth=0.5)
    ax1.set_xlabel('Insitu data [km/month]', fontsize=10)
    ax1.set_ylabel('TOPAZ rean. data [km/month]', fontsize=10)
    ax1.set_title(title)
    # Axis
    ax1.set_xlim(0,700)
    ax1.set_ylim(0,700)
    # Corr coef
    ax1.text(0.03, 0.95, 'Correlation: %3.1f' % icor, horizontalalignment='left', verticalalignment='center', transform=ax1.transAxes, fontsize=10)

    """
    ## right panel
    ax2 = fig.add_subplot(122)
    props = dict(alpha=0.5, edgecolors='none' )

    handles = []
    colors = ['blue', 'green', 'magenta', 'cyan']
    for color in colors:
        x = np.random.randn(N)
        y = np.random.randn(N)
        s = np.random.randint(50,200)
        handles.append(ax2.scatter(x, y, c=color, s=s, **props))

    ax2.set_ylim([-5,11])
    ax2.set_xlim([-5,11])

    ax2.legend(handles, colors)
    ax2.grid(True)
    ax2.set_aspect(1./ax2.get_data_ratio())
    """
    #plt.show()
    plt.savefig(out_path, dpi=150, bbox_inches='tight', pad_inches=0)

def get_and_plot(lon1,lat1,lon2,lat2,title,out_path):
    conn = psycopg2.connect(database='gis', user='gisuser')
    curs = conn.cursor()
    table_name = 'comparison'

    curs.execute("""SELECT date, mag_insitu, mag_model, st_AsText(geom) FROM %s WHERE st_within ( %s.geom, st_geomfromtext('POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s))') );""" % (table_name,table_name,lon1,lat1,lon1,lat2,lon2,lat2,lon2,lat1,lon1,lat1) )

    x_insitu = np.array([])
    y_model = np.array([])
    for row in curs.fetchall():
        i_date = row[0]
        i_mag_insitu = row[1]
        x_insitu = np.append(x_insitu, float(i_mag_insitu))
        i_mag_model = row[2]
        y_model = np.append(y_model, float(i_mag_model))
        #print i_mag_insitu, i_mag_model
    conn.commit()
    conn.close()

    corr = np.corrcoef(x_insitu, y_model)
    cor = corr[0][1]
    print 'Correlation: %s' % cor
    plot_data(x_insitu, y_model, out_path, title, cor)

try:
    # FS
    get_and_plot(-30,75,0,85,'Fram Strait','plots/fs.png')
    plot_poly(-30,75,0,85,'plots/map_fs.png')

    # Arctic
    get_and_plot(-180,60,180,90,'Arctic','plots/arctic.png')
    #plot_poly(-179,60,179,90,'plots/map_arctic.png')

    # Laptev and Siberian area
    get_and_plot(120,75,16,85,'Laptev and Siberain area','plots/laptev.png')
    plot_poly(120,75,16,85,'plots/map_laptev.png')

    # Beaufort area
    get_and_plot(-150,75,180,80,'Beaufort area','plots/bea.png')
    plot_poly(-150,75,-180,80, 'plots/map_bea.png')

    # Transpolar drift area
    get_and_plot(-60,85,160,90,'Transpolar Arctic drift area','plots/transpolar.png')
    plot_poly(0,85,160,90,'plots/map_transpolar.png')
except:
    print 'Error!'