ds = gdal.Open('albers_5km.tif')
data = ds.ReadAsArray()
gt = ds.GetGeoTransform()
proj = ds.GetProjection()
xres = gt[1]
yres = gt[5]
xmin = gt[0] + xres * 0.5
xmax = gt[0] + xres * ds.RasterXSize - xres * 0.5
ymin = gt[3] + yres * ds.RasterYSize + yres * 0.5
ymax = gt[3] - yres * 0.5
ds = None
xy_source = np.mgrid[xmin:xmax + xres:xres, ymax + yres:ymin:yres]
Plot GDAL raster using matplotlib Basemap