2

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?

Improve this question
  • If you do not need very many values you could try gdal.org/gdallocationinfo.html – user30184 Oct 29 '15 at 21:38
  • Complete map has a huge number of cells (933,120,000). How many cells (or rows and columns) are in the clipped raster of your study area? – xunilk Oct 29 '15 at 21:52
  • Thanks for your answer. Unfortunately I have too many cells in my study area. I have selected a specific area of the global map which has 2105x1388 cells. Is it possible to extract a shape file from the raster (tif) file so as to have for each value of the raster the corresponding point (lat/lon)? – nat Oct 30 '15 at 8:12
  • I edited my code to add the raster value in the text file. I hope that it helps. – xunilk Oct 30 '15 at 16:09
2

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).

enter image description here

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.

enter image description here

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!"
Improve this answer
  • Thank you very much for the useful answer! I will try that. In fact I wanted the txt file you created so to handle it in another model coded in FORTRAN, so it won't be a problem. However, I am a very new user of QGIS and not very familiar with python. Could you please attach the whole procedure of creating the txt file? That is to say, do I have to declare first some values or characters in python? What about the points.txt file? Where is it stored? Thank you again! – nat Oct 30 '15 at 11:47
  • The text file is in my default directory. You can specify your own path with: 'your_path/points.txt'. To run this code in the Python Console it is very easy. Copy/paste and save as, eg, coor_pyqgis.py in a well known directory. Open the Python Console/Editor and load the script. Select the raster as active Layer an run the script (last blue icon). By the way, you can accept the answer if it was useful. – xunilk Oct 30 '15 at 12:02
  • Thank you, it almost worked! It stops in line 6240 without printing the rest cells. Also, ID is the increasing number for each row/cell. I need to have also the value of the soil map (in raster file there is 1 band named DN i think). So I need this value together with the points. Is this possible? – nat Oct 30 '15 at 12:43
  • sorry, it was my fault, the output file is correct with all the points printed. However, what about the value "DN"? In fact, in the soil map is the value "MU_GLOBAL" which has to be linked with the soil data. – nat Oct 30 '15 at 12:55
  • Yes, it is possible but you also need to open the raster (with GDAL) to read its values (the easier way it is with 'ReadAsArray'). I will modify my answer later to include this. Now, I'm working with the long/lat raster. – xunilk Oct 30 '15 at 12:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.