Re: [Matplotlib-users] matplotlib/basemap plot geo data on background map

[Jeff W. <js...@fa...>]

  On 9/10/10 11:06 AM, izzybitsie wrote:
> Hi,
> I'm new to matplotlib and I'm looking for an easy way to plot geographical
> data on a background map: bkgmap.png
> http://old.nabble.com/file/p29679002/bkgmap.png
> So far I only found out about warpimage() to do this but only part of
> bkgmap.png comes up in the output image.  I think this is because this image
> has no pixels covering all the world.
> http://old.nabble.com/file/p29679002/partialbkg_polygon.png
>
> Any idea on how to insert this image as background?
> THANKS
>
> Code: the polygon displays in right position even though background doesn't
> show OK (tested with map,lat/lon lines drawn too)
>
> import sys
> import Image, ImageDraw  # PIL
> from matplotlib.patches import Polygon
> from mpl_toolkits.basemap import Basemap
> import matplotlib.image as mpimg
> import matplotlib.pyplot as plt
> import numpy as np
>
> lat0=48
> lon0=13
> lllon=-15
> lllat=20
> urlon=73
> urlat=57
> map =
> Basemap(projection='stere',lat_0=lat0,lon_0=lon0,llcrnrlon=lllon,llcrnrlat=lllat,urcrnrlon=urlon,urcrnrlat=urlat,
>                resolution='c',area_thresh=1000.)
> map.warpimage(image='bkgmap.png',scale=None,ax=plt.gca())
> # points
> lat = [50.,55.,45.,40.,50.]
> lon = [-20.,-10.,10.,-10.,-20.]
> x0,y0 = map(lon[0],lat[0])
> x1,y1 = map(lon[1],lat[1])
> x2,y2 = map(lon[2],lat[2])
> x3,y3 = map(lon[3],lat[3])
> x4,y4 = map(lon[4],lat[4])
> plt.gca().add_patch(Polygon([(x0,y0),(x1,y1),(x2,y2),(x3,y3),(x4,y4)],fill=1,facecolor='red',edgecolor='black'))
> plt.show()

Here is an example of how to make this work if the map projection on the 
image is exactly the same as the map projection in Basemap:

import sys
import Image, ImageDraw  # PIL
from matplotlib.patches import Polygon
from mpl_toolkits.basemap import Basemap
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.image import pil_to_array

lat0=48
lon0=13
lllon=-15
lllat=20
urlon=73
urlat=57
map =\
Basemap(projection='stere',lat_0=lat0,lon_0=lon0,llcrnrlon=lllon,llcrnrlat=lllat,urcrnrlon=urlon,urcrnrlat=urlat,
               resolution='l')
#map.warpimage(image='bkgmap.png',scale=None,ax=plt.gca())
pilImage = Image.open('bkgmap.png')
rgba = pil_to_array(pilImage)
im = map.imshow(rgba)
# points
lat = [50.,55.,45.,40.,50.]
lon = [-20.,-10.,10.,-10.,-20.]
x0,y0 = map(lon[0],lat[0])
x1,y1 = map(lon[1],lat[1])
x2,y2 = map(lon[2],lat[2])
x3,y3 = map(lon[3],lat[3])
x4,y4 = map(lon[4],lat[4])
plt.gca().add_patch(Polygon([(x0,y0),(x1,y1),(x2,y2),(x3,y3),(x4,y4)],fill=1,facecolor='red',edgecolor='black'))
map.drawmapboundary()
map.drawcoastlines()
plt.show()

Note that the coastlines drawn by drawcoastlines don't match the 
coastlines in the image (attached png).  This means the Basemap 
projection doesn't match the image projection.  Unless you can make them 
match, you'll have to write your own transformations.  Of course, you 
can't write your own transformation unless you know exactly what the map 
projection of the image is, so the simplest approach will be to make the 
Basemap projection match the image.

-Jeff

-- 
Jeffrey S. Whitaker         Phone  : (303)497-6313
Meteorologist               FAX    : (303)497-6449
NOAA/OAR/PSD  R/PSD1        Email  : Jef...@no...
325 Broadway                Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web    : http://tinyurl.com/5telg