Yes, it is possible to extract a shapefile but it is very difficult to handle it because you need a lot a memory. It is preferable to "burn" these coordinates in two raster (one for each coordinates values; x and y).
However, I tried to produce the "shapefile" (in fact an "event theme" from a text file) with next approach. First, I selected an arbitrary area (at Brazil) with your exact cell number (2105x1388 cells) and helped with my two plugins (as it can be observed at the next image).
With the next code ran at the Python Console of QGIS, I got a text file to be loaded in QGIS.
layer = iface.activeLayer()
extent = layer.extent()
xmin = extent.xMinimum()
ymax = extent.yMaximum()
rows = layer.height()
columns = layer.width()
xsize = layer.rasterUnitsPerPixelX()
ysize = layer.rasterUnitsPerPixelY()
k = 1
xinit = xmin + xsize/2
yinit = ymax - ysize/2
pfile = open('points.txt', 'w')
pfile.write('id x y \n')
print "Processing..."
for i in range(rows):
for j in range(columns):
x = xinit + j*xsize
y = yinit
pfile.write(str(k) + " " + str(x) + " " + str(y) + "\n")
k +=1
xinit = xmin + xsize/2
yinit -= ysize
pfile.close()
At the next image, you can observe a point layer (cell centres in lon/lat) created as delimited text file but, it is very difficult to handle. It is possible but inconvenient.

Editing Note:
Next code add the raster values to the text file:
from osgeo import gdal
import os
layer = iface.activeLayer()
provider = layer.dataProvider()
my_path = provider.dataSourceUri()
dataset = gdal.Open(my_path)
band = dataset.GetRasterBand(1)
data = band.ReadAsArray(0, 0, band.XSize, band.YSize)
extent = layer.extent()
xmin = extent.xMinimum()
ymax = extent.yMaximum()
rows = layer.height()
columns = layer.width()
xsize = layer.rasterUnitsPerPixelX()
ysize = layer.rasterUnitsPerPixelY()
k = 1
xinit = xmin + xsize/2
yinit = ymax - ysize/2
pfile = open('points.txt', 'w')
pfile.write('id x y \n')
print "Processing..."
for i in range(rows):
for j in range(columns):
x = xinit + j*xsize
y = yinit
pfile.write(str(k) + " " + str(x) + " " + str(y) + " " + str(data[i][j]) + "\n")
k +=1
xinit = xmin + xsize/2
yinit -= ysize
pfile.close()
dataset = None
print "Done!"