import mpl_toolkits
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
from osgeo import gdal
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)
plt.imshow(msk_array)
<matplotlib.image.AxesImage at 0x99401f0>
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')
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')
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')
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.
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