1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
|
from __future__ import (absolute_import, division, print_function)
"""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 mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
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 = plt.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=plt.cm.gray, interpolation='nearest')
plt.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(np.arange(10,76,5), labels=[0,0,1,0], color=overlay_color)
m.drawparallels(np.arange(-90,90,5), labels=[1,0,0,0], color=overlay_color)
# add timestamp and save
fig = plt.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))
plt.title('Meteosat Geostationary Satellite Image - Portion of Full Earth',y=1.05,fontsize=12)
plt.show()
#fig.savefig(plot_name)
#print 'Plot saved to %s' % (plot_name)
|
