"""geos_demo_2.py

This script shows how to plot data onto the Geostationary Satellite projection
when the data is from a portion of the full Earth image. The script assumes that
the data is already contained in a regular grid in the geos projection and that
the corner points of the data grid are known in lat-long.

Dependencies: Matplotlib, Basemap toolkit, Python Imaging Library

"""
from PIL import Image
from matplotlib.toolkits.basemap import Basemap
from pylab import cm, arange, colorbar, text, savefig, gcf, clim, title, show
from matplotlib.image import pil_to_array

plot_name = 'geos_demo.png'
overlay_color = 'black'

# read in jpeg image to rgb array
pilImage = Image.open('200706041200-msg-ch01-SAfrica.jpg')
#data = asarray(pilImage)
data = pil_to_array(pilImage)
data = data[:, :, 0] # get data from first channel in the image

# define data region and projection parameters
ll_lon = 9.74
ll_lat = -35.55
ur_lon = 48.45
ur_lat = 0.2
lon_0 = 0.0
satellite_height = 35785831.0

# create Basemap instance for a Geostationary projection.
m = Basemap(projection='geos', lon_0=lon_0, satellite_height=satellite_height,
            resolution='l', llcrnrlon=ll_lon, llcrnrlat=ll_lat, urcrnrlon=ur_lon, urcrnrlat=ur_lat)
# add data
m.imshow(data, cmap=cm.gray, interpolation='nearest')
# draw coastlines.
m.drawcoastlines(linewidth=0.5, color=overlay_color)
m.drawcountries(linewidth=0.5, color=overlay_color)
# can't label meridians on bottom, because labels would
# be outside map projection region.
m.drawmeridians(arange(10,76,5), labels=[0,0,1,0], color=overlay_color)
m.drawparallels(arange(-90,90,5), labels=[1,0,0,0], color=overlay_color)
# add a colobar
clim(0, 255)
colorbar()
# add timestamp and save
fig = gcf()
fig.text(x=0.275, y=0.025, s=u'Meteosat-9 VIS 0.6 channel - 12:00 UTC 04/06/2007\n    \N{COPYRIGHT SIGN} EUMETSAT 2007',
            horizontalalignment='left',
            verticalalignment='bottom',
            fontsize=10,
            fontweight='bold',
            bbox=dict(facecolor='gray', alpha=0.25, pad=15))
fig.set_size_inches((8, 6))
title('Meteosat Geostationary Satellite Image - Portion of Full Earth',y=1.05,fontsize=12)

show()
#fig.savefig(plot_name)
#print 'Plot saved to %s' % (plot_name)

"""geos_demo_2.py

This script shows how to plot data onto the Geostationary Satellite projection
when the data is from a portion of the full Earth image. The script assumes that
the data is already contained in a regular grid in the geos projection and that
the corner points of the data grid are known in lat-long.

Dependencies: Matplotlib, Basemap toolkit, Python Imaging Library

"""
from PIL import Image
from matplotlib.toolkits.basemap import Basemap
from pylab import figure, cm, arange, colorbar, text, savefig, gcf, clim, title, show
from matplotlib.image import pil_to_array

plot_name = 'geos_demo.png'
overlay_color = 'black'

# read in jpeg image to rgb array
pilImage = Image.open('200706041200-msg-ch01-SAfrica.jpg')
#data = asarray(pilImage)
data = pil_to_array(pilImage)
data = data[:, :, 0] # get data from first channel in the image

# define data region and projection parameters
ll_lon = 9.74
ll_lat = -35.55
ur_lon = 48.45
ur_lat = 0.2
lon_0 = 0.0
satellite_height = 35785831.0

fig = figure(figsize=(7,7))
ax = fig.add_axes((0.1,0.1,0.8,0.8))
# create Basemap instance for a Geostationary projection.
m = Basemap(projection='geos', lon_0=lon_0, satellite_height=satellite_height,
            resolution='l', llcrnrlon=ll_lon, llcrnrlat=ll_lat, urcrnrlon=ur_lon, urcrnrlat=ur_lat)
# add data
m.imshow(data, cmap=cm.gray, interpolation='nearest')
clim(0, 255)
# draw coastlines.
m.drawcoastlines(linewidth=0.5, color=overlay_color)
m.drawcountries(linewidth=0.5, color=overlay_color)
# can't label meridians on bottom, because labels would
# be outside map projection region.
m.drawmeridians(arange(10,76,5), labels=[0,0,1,0], color=overlay_color)
m.drawparallels(arange(-90,90,5), labels=[1,0,0,0], color=overlay_color)
# add a colobar
#colorbar()
# add timestamp and save
fig = gcf()
fig.text(x=0.275, y=0.025, s=u'Meteosat-9 VIS 0.6 channel - 12:00 UTC 04/06/2007\n    \N{COPYRIGHT SIGN} EUMETSAT 2007',
            horizontalalignment='left',
            verticalalignment='bottom',
            fontsize=10,
            fontweight='bold',
            bbox=dict(facecolor='gray', alpha=0.25, pad=15))
fig.set_size_inches((8, 6))
title('Meteosat Geostationary Satellite Image - Portion of Full Earth',y=1.05,fontsize=12)

show()
#fig.savefig(plot_name)
#print 'Plot saved to %s' % (plot_name)

remove colorbar