Import the libraries

In [102]:
import mpl_toolkits
from mpl_toolkits.basemap import Basemap 
import matplotlib.pyplot as plt 
import numpy as np
from osgeo import gdal

Load raster using gdal, mask NoData values. Get projection info

In [122]:
pathToRaster = r'I:\Data\anomaly//ano_DOY2002170.tif'
raster = gdal.Open(pathToRaster, gdal.GA_ReadOnly)
array = raster.GetRasterBand(1).ReadAsArray()
msk_array = np.ma.masked_equal(array, value = 65535)
print 'Raster Projection:\n', raster.GetProjection()
print 'Raster GeoTransform:\n', raster.GetGeoTransform()
Raster Projection:
GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]
Raster GeoTransform:
(-180.0, 0.25, 0.0, 90.0, 0.0, -0.25)

Quick plot using pyplot

In [4]:
plt.imshow(msk_array)
Out[4]:
<matplotlib.image.AxesImage at 0x99401f0>

Let's plot a real map using Basemap

In [126]:
map = Basemap(projection='cyl',resolution='c',lat_0=0,lon_0=0)

# Add some additional info to the map
cstl = map.drawcoastlines(linewidth=.5)
meri = map.drawmeridians(np.arange(0,360,60), linewidth=.2, labels=[1,0,0,1], labelstyle='+/-', color='grey' ) 
para = map.drawparallels(np.arange(-90,90,30), linewidth=.2, labels=[1,0,0,1], labelstyle='+/-', color='grey')
boun = map.drawmapboundary(linewidth=0.5, color='grey')

Add the data as an image

In [88]:
map = Basemap(projection='cyl',resolution='c',lat_0=0,lon_0=0)

datain = np.flipud( msk_array )
imsh = map.imshow(datain)

# Add some more info to the map
cstl = map.drawcoastlines(linewidth=.5)
meri = map.drawmeridians(np.arange(0,360,60), linewidth=.2, labels=[1,0,0,1], labelstyle='+/-', color='grey' ) 
para = map.drawparallels(np.arange(-90,90,30), linewidth=.2, labels=[1,0,0,1], labelstyle='+/-', color='grey') 

Add the same data as filled contours

In [124]:
map = Basemap(projection='cyl',resolution='c',lat_0=0,lon_0=0)

datain = np.flipud( msk_array )
ny = datain.shape[0]; nx = datain.shape[1]
lons,lats = map.makegrid(nx,ny) # get lat/lons of ny by nx evenly space grid.
x,y = map(lons,lats) # compute map prj coordinates

levels = [-1000,-800,-600,-400,-200,0,200,400,600,800,1000]
cntr = map.contourf(x,y,datain, levels,cmap=cm.RdBu)
cbar = map.colorbar(cntr,location='bottom',pad='15%')

# Add some more info to the map
cstl = map.drawcoastlines(linewidth=.5)
meri = map.drawmeridians(np.arange(0,360,60), linewidth=.2, labels=[1,0,0,1], labelstyle='+/-', color='grey' ) 
para = map.drawparallels(np.arange(-90,90,30), linewidth=.2, labels=[1,0,0,1], labelstyle='+/-', color='grey') 
boun = map.drawmapboundary(linewidth=0.5, color='grey')

Problems starts here: Now I want to use Robinson projection instead of the Cylindrical Equidistant, but this is giving me strange output

In [125]:
map = Basemap(projection='robin',resolution='c',lat_0=0,lon_0=0)

datain = np.flipud( msk_array )
ny = datain.shape[0]; nx = datain.shape[1]
lons,lats = map.makegrid(nx,ny) # get lat/lons of ny by nx evenly space grid.
x,y = map(lons,lats) # compute map prj coordinates

levels = [-1000,-800,-600,-400,-200,0,200,400,600,800,1000]
cntr = map.contourf(x,y,datain, levels,cmap=cm.RdBu)
cbar = map.colorbar(cntr,location='bottom',pad='15%')

# Add some additional info to the map
cstl = map.drawcoastlines(linewidth=.5)
meri = map.drawmeridians(np.arange(0,360,60), linewidth=.2, labels=[1,0,0,1], labelstyle='+/-', color='grey' ) 
para = map.drawparallels(np.arange(-90,90,30), linewidth=.2, labels=[1,0,0,1], labelstyle='+/-', color='grey')
boun = map.drawmapboundary(linewidth=0.5, color='grey')
WARNING: x coordinate not montonically increasing - contour plot
may not be what you expect.  If it looks odd, your can either
adjust the map projection region to be consistent with your data, or
(if your data is on a global lat/lon grid) use the shiftgrid
function to adjust the data to be consistent with the map projection
region (see examples/contour_demo.py).

Map has vertical lines in the northern longitudes and data is offset compare to coastlines.

Now trying the same using the basemap.interp function

In [123]:
map = Basemap(projection='robin',resolution='c',lat_0=0,lon_0=0)

datain = np.flipud( msk_array )

nx = raster.RasterXSize
ny = raster.RasterYSize

xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of x points on the grid
yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of y points on the grid

lons = np.arange(-180,180,0.25) #from raster.GetGeoTransform()
lats  = np.arange(-90,90,0.25) 

lons, lats = np.meshgrid(lons,lats) 
xout,yout = map(lons, lats)
dataout = mpl_toolkits.basemap.interp(datain, xin, yin, xout, yout, order=1)

levels = [-1000,-800,-600,-400,-200,0,200,400,600,800,1000]
cntr = map.contourf(xout,yout,dataout, levels,cmap=cm.RdBu)
cbar = map.colorbar(cntr,location='bottom',pad='15%')

# Add some more info to the map
cstl = map.drawcoastlines(linewidth=.5)
meri = map.drawmeridians(np.arange(0,360,60), linewidth=.2, labels=[1,0,0,1], labelstyle='+/-', color='grey' ) 
para = map.drawparallels(np.arange(-90,90,30), linewidth=.2, labels=[1,0,0,1], labelstyle='+/-', color='grey')
boun = map.drawmapboundary(linewidth=0.5, color='grey')

There is no error message and map has no vertical lines anymore, but data is still offset compare to coastlines

Next step: clueless...