Convert raster to a CSV values and lat long based on raster value in QGIS

Question

I am looking to convert a raster into a csv file (or tab delimited) based on values in the raster (Value > 0).

The csv would have three values, latitude of center of cell, longitude of center of cell, and the raster value where the raster value is greater than 0.

The data in unprojected in geographic coordinate system WGS 1984. It is a GTiff.

I know I can go to point file and then add the latitude and longitude and then go to CSV but I wonder if a direct method exists as I have 130 million points.

Working on going Ras > Point (file Geodatabase) > Query Points > CSV in ArcGIS but it is slow.

I feel I need to build in the query at the Ras > Point section to speed things up as the entire raster has 130 billion cells but only 130 million have values. Rest are 0 or NoData.

Looking for QGIS with processing (Grass, SAGA, etc.). Have OSGEO4W installed but can bring up a Linux box as well.

UPDATE As per @Luke answer below I am running this GDAL code in OSGEO4w.

E:\temp>gdal_translate -q -of xyz -co ADD_HEADER_LINE=YES -co COLUMN_SEPARATOR="
," mwf2012.tif /vsistdout | findstr /V /C:"65536" > output.csv

I am not getting the usal GDAL Progress bar (is this normal?), everything else looks good so Ill Let you know the final result!

Answer 1

As GDAL supports writing to X,Y,Z (CSV) ascii, you could use gdal_translate:

gdal_translate -of xyz -co ADD_HEADER_LINE=YES -co COLUMN_SEPARATOR="," input_raster output.csv

To avoid writing NoData values to your output you can write the output to stdout then pipe to grep/findstr to filter it before writing to your csv:

gdal_translate -q -of xyz -co ADD_HEADER_LINE=YES -co COLUMN_SEPARATOR="," input_raster /vsistdout | findstr /V /C:"your_nodata_value" > output.csv 

Or in a *Nix OS:

gdal_translate -q -of xyz -co ADD_HEADER_LINE=YES -co COLUMN_SEPARATOR="," input_raster /vsistdout | grep -v "your_nodata_value" > output.csv 

Running the above command (inc. findstr/grep) on a 60 million pixel (7700x7800 Landsat scene) took 2 minutes on my Windows PC and 1min 15 on my slightly more powerful Linux workstation.

Notes:

1. my nodata_value is negative so I had to escape the negative character on Linux: gdal_translate etc...|grep -v "\-999" > output.csv
2. if your nodata_value is something really common, i.e "0", you'll need to do something a bit more complex in the grep/findstr expression, like ...|grep -v ",0$" or ...|findstr /V /C:",0$"

Answer 2

Try this way with GDAL:

import os, sys

from osgeo import gdal
from osgeo import gdalconst

# get the arguments
InRaster = sys.argv[1]
OutCSV   = sys.argv[2]

# open the raster and get some properties
ds       = gdal.OpenShared(InRaster,gdalconst.GA_ReadOnly)
GeoTrans = ds.GetGeoTransform()
ColRange = range(ds.RasterXSize)
RowRange = range(ds.RasterYSize)
rBand    = ds.GetRasterBand(1) # first band
nData    = rBand.GetNoDataValue()
if nData == None:
    nData = -9999 # set it to something if not set
else:
    print("NoData value is {0}".format(nData))

# specify the centre offset
HalfX    = GeoTrans[1] / 2
HalfY    = GeoTrans[5] / 2

with open(OutCSV,'w') as CSVwrite:
    for ThisRow in RowRange:
        RowData = rBand.ReadAsArray(0,ThisRow,ds.RasterXSize,1)[0]
        for ThisCol in ColRange:
            if RowData[ThisCol] != nData:
                if RowData[ThisCol] > 0:
                    X = GeoTrans[0] + ( ThisCol * GeoTrans[1] )
                    Y = GeoTrans[3] + ( ThisRow * GeoTrans[5] ) # Y is negative so it's a minus
                    # this gives the upper left of the cell, offset by half a cell to get centre
                    X += HalfX
                    Y += HalfY
                    CSVwrite.write('{0},{1},{2}\n'.format(X,Y,RowData[ThisCol])) # I think newline might be different on linux

This reads up each line (row) of the raster one at a time then when it finds the value in the row writes the coords to a CSV file.

This will work just as well in DD as metres, feet, inches etc.. as it's using the defined GeoTransform object which are numbers, smaller in geographic but still numbers just the same.

If you use gdal_translate -of XYZ, as you've discovered, it will export the NoData records as points... I've seen this before and this is the script I use to fix the XYZ file(s); it also converts PRN (space delimited) to CSV.

import os, sys

InF  = sys.argv[1]
OutF = sys.argv[2]

with open(InF,'r') as InFile:
    with open(OutF,'w') as OutFile:
        for ThisLine in InFile:

            lSp = ThisLine.split(" ")
            if lSp[2] != "-9999\n":
                OutFile.write(ThisLine.replace(' ',','))