Wednesday 27 November 2013

Map drawing

#########################################################
#  Plotting points on Lambert Conformal projection with blue marble image
#
#
#  JHPS
#  C: 26/11/2013

#######################################################

from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd


#Points
stops = pd.read_csv('stops.csv')      #This is just a csv (e.g. excel) file with columns 'NAME', 'LON', and 'LAT'

num = 0                                           # num picks the element in the stops list that you want to plot (starting from 0)
lon_c = stops['LON'][num]
lat_c = stops['LAT'][num]

#Setup
fig = plt.figure()
#map = Basemap(resolution='c',projection='ortho',lat_0=60.,lon_0=-60.)
ax = fig.add_axes([0.1,0.1,0.8,0.8])

map = Basemap(llcrnrlon=lon_c - 60,llcrnrlat=30,urcrnrlon=lon_c + 60,urcrnrlat=70.000,\
            rsphere=(6378137.00,6356752.3142),\
            resolution='l',area_thresh=1000.,projection='lcc',\
            lat_1=10.,lon_0=lon_c,ax=ax)

map.drawparallels(np.arange(-90.,91.,30.))
map.drawmeridians(np.arange(-180.,181.,60.))

parallels = np.arange(0.,80,20.)

map.drawparallels(parallels,labels=[1,0,0,1])
meridians = np.arange(10.,360.,30.)
map.drawmeridians(meridians,labels=[1,0,0,1])

map.bluemarble()

map.drawcountries(linewidth=1)


x,y = map(stops['LON'][1],stops['LAT'][1])

map.plot(x,y,'r.',markersize=20)

plt.show()



Tuesday 22 October 2013

Data assimilation: Interpolation example

########################################################
# Demonstration of interpolation for Hercule's cyborg
#  prints three plots for different input data and resulting
#  interpolations
#
#
# JHPS
# C: 22/10/2013
#######################################################

import numpy as np
from scipy.interpolate import Rbf
import matplotlib.pyplot as plt
from matplotlib import cm

### After first walk with points A, B and C
x = np.array([-1,1.5,-1])
y = np.array([1,0,-1])
z = np.array([5,0,-5])

ti = np.linspace(-2.0, 2.0, 100)
XI, YI = np.meshgrid(ti, ti)

# Interpolation
rbf = Rbf(x, y, z, epsilon=2)
ZI = rbf(XI, YI)

# Plotting
plt.figure()
n = plt.normalize(-2., 2.)
plt.subplot(1, 1, 1)
plt.pcolor(XI, YI, ZI, cmap=cm.jet)
plt.scatter(np.array([0]), np.array([0]), s=100, facecolors='none', edgecolors='k')
plt.scatter(x, y, 100, z, cmap=cm.jet)
plt.plot(np.array([0]),np.array([0]))
plt.xlim(-2, 2)
plt.ylim(-2, 2)
plt.colorbar()

labels = ['A','B','C']

for i in range(len(x)):
        plt.annotate(labels[i], xy=(x[i], y[i]),  xycoords='data',
            xytext=(-20, 10), textcoords='offset points',
            arrowprops=dict(arrowstyle="->")
            )

plt.annotate('D', xy=(0, 0),  xycoords='data',
            xytext=(-20, 10), textcoords='offset points',
            arrowprops=dict(arrowstyle="->")
            )


# Second walk with another point between A and B
x = np.array([-1,1.5,-1,0.5])
y = np.array([1,0,-1,1.5])
z = np.array([5,0,-5,3])

ti = np.linspace(-2.0, 2.0, 100)
XI, YI = np.meshgrid(ti, ti)

# Interpolation
rbf = Rbf(x, y, z, epsilon=2)
ZI = rbf(XI, YI)

# Plotting
plt.figure()
n = plt.normalize(-2., 2.)
plt.subplot(1, 1, 1)
plt.pcolor(XI, YI, ZI, cmap=cm.jet)
plt.scatter(np.array([0]), np.array([0]), s=100, facecolors='none', edgecolors='k')
plt.scatter(x, y, 100, z, cmap=cm.jet)
plt.plot(np.array([0]),np.array([0]))
plt.xlim(-2, 2)
plt.ylim(-2, 2)
plt.colorbar()

plt.annotate('New point', xy=(x[-1], y[-1]),  xycoords='data',
            xytext=(20, -10), textcoords='offset points',
            arrowprops=dict(arrowstyle="->")
            )


plt.annotate('D', xy=(0, 0),  xycoords='data',
            xytext=(-20, 10), textcoords='offset points',
            arrowprops=dict(arrowstyle="->")
            )


# Third walk with another point before A
x = np.array([-1.5,-1,1.5,-1])
y = np.array([1.5,1,0,-1])
z = np.array([-5,5,0,-5])

ti = np.linspace(-2.0, 2.0, 100)
XI, YI = np.meshgrid(ti, ti)

# Interpolation
rbf = Rbf(x, y, z, epsilon=2)
ZI = rbf(XI, YI)

# Plotting
plt.figure()
n = plt.normalize(-2., 2.)
plt.subplot(1, 1, 1)
plt.pcolor(XI, YI, ZI, cmap=cm.jet)
plt.scatter(np.array([0]), np.array([0]), s=100, facecolors='none', edgecolors='k')
plt.scatter(x, y, 100, z, cmap=cm.jet)
plt.plot(np.array([0]),np.array([0]))
plt.xlim(-2, 2)
plt.ylim(-2, 2)
plt.colorbar()


plt.annotate('a different new point', xy=(x[0], y[0]),  xycoords='data',
            xytext=(20, 0), textcoords='offset points',
            arrowprops=dict(arrowstyle="->")
            )


plt.annotate('D', xy=(0, 0),  xycoords='data',
            xytext=(-20, 10), textcoords='offset points',
            arrowprops=dict(arrowstyle="->")
            )

## show plots

plt.show()