How to extract coordinates from a raster in QGIS

Question

I have the Harmonised Global Soil Database in raster format (tif). I want to extract the coordinates (cell centres in lon/lat) along with the field value in order to link them with the soil data containing in the mdb file.

Is there any way to extract coordinates from a raster file?

Answer 1

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!"