KartozaSomeone wrote to me asking if it would be possible to generate an XYZ ASCII file from a single band raster layer in QGIS. No doubt there are more efficient ways (this approach is pretty slow but it works), but I thought it would be fun to show how you can iterate over a raster, writing out the value of each cell into a text file (along with the centroid coordinates for that cell).
To use the script, you should save it to your local machine, then open the python console and load the script in the python editor. Next select a single band raster and then run the script in the editor. If your raster is quite large, it will take some time to run. I have spent zero time trying to optimise the script - if someone has an idea for doing it faster, send me a patch and I will update the example above.
The generated output dataset will look something like this:
Longitude,Latitude,VI
8.31259406548,7.86128343221,10
8.31264849753,7.86128343221,16
8.31270292958,7.86128343221,18
8.31248520138,7.8613378416,15
8.31253963343,7.8613378416,17
8.31259406548,7.8613378416,24
8.31264849753,7.8613378416,46
8.31270292958,7.8613378416,47
...
The resulting script can be used with programs like gdal_grid or loaded back into QGIS as a vector layer using the Delimited Text provider:
Comments
Stefan Keller 9 years, 8 months ago
Hi Tim!
Link | ReplyThat blog post is a fun Python code snippet! I'd like to add following comments and questions:
Why not simply use QGIS menu "Raster->Conversion->Translate"? (which in fact to me should be moved to layer "Save as...")
Given one prefers scripting, I would point explicitly to GDAL gdal_translate.
And given one prefers more a GUI what one could add a QtGui.QFileDialog(...) before line 9 "f = file('/tmp/test.xyz', 'w')"?
Yours, Stefan
Comment awaiting approval 5 years ago
Comment awaiting approval 5 years ago
Comment awaiting approval 4 years, 11 months ago
Comment awaiting approval 4 years, 11 months ago
Tim Sutton 9 years, 8 months ago
Hi Stefan
Link | ReplyYes gdal_translate is muuuuuuch faster so if you are only interested in performance you should definately use it. If you are interested in a knowing how to step wise iterate through a raster and do something with each, cell's value using the QGIS api then this snippet if aimed at you. Note that it would be faster also to use QGIS raster block API - you can see I actually get the block out in line 9:
> b = p.block(0, p.extent(), w, h)
but then don't use it (sloppy code I know :-) ). It would be a nice topic for a future post - how to read and use raster blocks from QGIS in python.
In this specific case, the person wanted all the 0 values excluded from the output file, and specific headers written (both are probably possible using gdal if I dug a little around the options).
As for the GUI file selector, I would probably make a proper plugin of this snippet and have an icon and a file dialog if I wanted to spend more time on it.
Thanks for your comments!
Tim
FFXIV 9 years, 4 months ago
<strong>FFXIV</strong>
Link | ReplyIf he shows up, we stick with him.
New Comment