3

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!

  • You said you have Arc.. right? If you have Spatial Analyst use Extract by Attributes (Value > 0) so all you have is valid values and NoData then use Raster to Point, which will ignore NoData cells, then Add XY to get the lat/lon for each point. If you don't have Arc AND Spatial Analyst your best bet is with GDAL in python or RasterToNumpyArray resources.arcgis.com/en/help/main/10.2/index.html#/… to speed up iterating cells. – Michael Stimson Jul 7 '16 at 2:56
  • @MichaelMiles-Stimson thnaks, yep that is my Arc process 100% but 123 million points is a pain to create. Breaks Shapefile size rules so needs to be geodatabase but that is fine just dog slow. Not sure how to do it it Numpy (although I am sure it is the way to go) without loading the entire raster into RAM (~600GB). – If you do not know- just GIS Jul 7 '16 at 3:01
  • RasterToNumpyArray lets you specify the origin, rows and cols to load. I'd do it one row (or col) at a time and iterate the col (or row) through the array, when you detect a value you're interested in calculate the lat/lon as the row/col * cell size offset from the origin of the raster. I've not used this one but I have done a similar process in GDAL - just remember the origin in GDAL is the upper left and cell Y is negative; Y = GeoTrans[3] + (row * GeoTrans[5]) and X = GeoTrans[0] + (col * GeoTrans[1]) – Michael Stimson Jul 7 '16 at 3:26
  • @MichaelMiles-Stimson very helpful as usual. Going to try gdal alone first then gdal numpy and i have been meaning to learn how to better use numpy in this iterative manner. – If you do not know- just GIS Jul 7 '16 at 3:28
7

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$"
  • That's handy to know Luke.. I didn't realize there was a way around writing the NoData cells as points; that's certainly a handy hint about needing to escape the negative. +1 from me! – Michael Stimson Jul 7 '16 at 5:24
  • +1 from me too. This is my first time to know that you can avoid exporting no data while exporting using gdal_translate. – ahmadhanb Jul 7 '16 at 5:53
  • @Luke, wow. this is going to be my approach. Going to have to go Windows OSGEO4w route due to not having enough storage attached to linux to run. Ill update question. No progress bar, is this normal? Bethinking else looks ggod. GDAL can appear magic. – If you do not know- just GIS Jul 7 '16 at 12:15
  • @Ifyoudonotknow-justGIS GDAL usually writes the progress bar to stdout. You don't get a progress bar because the "-q" argument stops GDAL writing anything (except the output x,y,z values) to stdout, so only your raster x, y, z values get piped to findstr/grep and then to your output CSV. If you want progress, you could do the conversion in python and use a gdal.TermProgress callback function to print to stderr. – user2856 Jul 7 '16 at 21:47
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(' ',','))
  • +1 for providing a python code to process -9999 – ahmadhanb Jul 7 '16 at 6:09

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.