topic Re: Extract x,y,z from Raster - script in Geoprocessing Questions

Extract x,y,z from Raster - script

JohnLonsdale — Mon, 11 Oct 2010 12:16:16 GMT

I have raster bathymetry data that i would like to extract the x,y,z values from, whilst also ignoring No Data cells (which on narrow survey passes at a 45° angle to North can be over half the data points). I am an experienced ArcMap 9.3.1 end user but a little rough around the gills in writing Python 2.5 scripts. However I have come up with the following:

import arcgisscripting
gp = arcgisscripting.create()

InRaster = "C:\JPL_Jobs\XYZ_TEST\AreaTest2_Clip.img"

wf = open('C:\JPL_Jobs\XYZ_TEST\Create_Test_xyz.txt','w')

Cellsize = gp.GetRasterProperties(InRaster, "CELLSIZEX")
Xul = gp.GetRasterProperties(InRaster, "LEFT")
Yul = gp.GetRasterProperties(InRaster, "TOP")
Ht = gp.GetRasterProperties(InRaster, "ROWCOUNT")
Wd = gp.GetRasterProperties(InRaster, "COLUMNCOUNT")

wf.write("x,y,z\n")
try:
    for j in range(int(Ht)):
        for i in range(int(Wd)):
            XValue = Xul + ((Cellsize * i) + (Cellsize / 2))
            YValue = Yul - ((Cellsize * j) + (Cellsize / 2))
            XValueStr = str(XValue)
            YValueStr = str(YValue)
            d = XValueStr + " " + YValueStr
            cellVal = gp.GetCellValue_management(InRaster, d)
            if cellVal <> "NoData" :
                cellValStr = str(cellVal)
                lineval = XValueStr + ", " + YValueStr + ", " + cellValStr + "\n"
                wf.write(lineval)
    
except:
    print gp.GetMessages()

wf.close()
print "Finished"

This works but my question is, on a very small 6x6 grid it took 11 secs. I am currently looking at a 9400 x 9300 bathymetry raster (a typical size) and wonder whether using this method just simply won't be up to the task, taking far too long to process? Or am I going about this task in the wrong way? I have used the Raster2xyz vba script in the past which has worked fine but i believe ArcGIS 10 is not supporting VB. Our team is switching soon but I am not keen (just yet!).

As an aside, typing 9.3 into arcgisscripting.create() caused the script to fail. Any thoughts on why that might be?

Thank you in advance
John

Re: Extract x,y,z from Raster - script

ChrisSnyder — Mon, 11 Oct 2010 15:44:19 GMT

I would just use the RasterToPoint tool. NoData values are not exported... Specify a .shp output as that is the vector format that SpatialAnalyst writes to directly (FGDB is a post process!). Took my machine 3 min 14sec for a 7364 x 4454 cell raster. Be mindful of keeping the .dbf file under 2.1 GB. Once in .shp (.dbf) format, you quite readily have your x,y,z database.

I should note that RasterToAscii only took 16 seconds overall to run on the same raster, but you would need to include the added overhead of writing/processing a function to give you the cell center x/y pairs.

Re: Extract x,y,z from Raster - script

Luke_Pinner — Sun, 12 Dec 2021 07:13:22 GMT

Use the Spatial Analyst "sample" tool/function - http://webhelp.esri.com/arcgisdesktop/9.3/index.cfm?id=6230&pid=6215&topicname=Sample

E.g in the raster calculator (as at 9.3)

c:\temp\test_zxy.txt = sample(somegrid)

The output is tab-delimited ZXY format.

somegrid    x    y
0.2306625    0.5    9.5
0.2572387    1.5    9.5
0.3041072    2.5    9.5
0.01269815    3.5    9.5
0.9818037    4.5    9.5
etc...

Re: Extract x,y,z from Raster - script

JohnLonsdale — Tue, 12 Oct 2010 07:46:37 GMT

Wow, that is excellent. Thanks for the replies Bill, Chris and Luke. I needed to know if i was going down a blind alley and it appears I was, but it was an interesting exercise for me anyway.

Bill - that was one of my reasons for posting, to get some expert knowledge on programming intricacies built into ArcGIS and issues with the programming process, so I greatly appreciate that. An interesting thing i discovered was that processing a 76x50 grid took 19 min 31 secs from PythonWin yet when attached to an ArcToolbox script it took 8 min 57 sec. Still a ridiculously long time but something of note. We get data from our survey team as an ArcGIS ascii file but sometimes it is adjusted and clients prefer the x,y,z format. So as i get more confident with Python I can play around with your suggestion.

Chris - don't you find that then having to add the x and y fields to the point .shp and populating them can take a long time. I've tried this before and it seemed a bit cumbersome, unless I was doing it wrong - or should I do this out in a programming sequence or modelbuilder? And you are right we often touch the 2GB limits with data.

Luke - thanks for that. We have a demo version of Spatial Analyst so i will have a look at that. The need for SA within our small team seems to grow. Being a little bit more expensive we have held off and sought out alternatives, like Hawths tools or whatever it is called now (GME?).

Again, thank you all for your time and help.
John

Re: Extract x,y,z from Raster - script

ChrisSnyder — Sun, 12 Dec 2021 07:13:24 GMT

Here's some Python code that will take an ascii file (output of the ArcGIS RastertoAscii tool, for which no SA license is needed), and reformat it to a .csv file with x, y, z values (x and y are cell centers). Note that depending on the particulars of your DEM, the .csv file will generally be much larger than the input ascii file. This script could easily be altered so as to write the xyz stuff to a binary format like a .dbf, .shp, FGDB, etc.

# Description: Format an ascii text file (created using the RasterToAscii in ArcGIS) as a comma delimited .csv file containing x, y, z values
# Written By: Chris Snyder, 20101012
import os
inputAsciiDemFile = r"D:\csny490\dems\dem90m\dem90.txt" #ascii file created using the RasterToAscii in ArcGIS 
outputCsvFile = r"D:\csny490\temp\test.csv" #output comma delimited file ascii file (x,y,z fields)
if os.path.exists(outputCsvFile):
    os.remove(outputCsvFile)
inputLines = open(inputAsciiDemFile, 'r')
nColumns = int(inputLines.next().split(" ")[-1])
nRows = int(inputLines.next().split(" ")[-1]) 
xCorner = float(inputLines.next().split(" ")[-1]) 
yCorner = float(inputLines.next().split(" ")[-1]) 
cellSize = float(inputLines.next().split(" ")[-1]) 
nodataValue = int(inputLines.next().split(" ")[-1])
print >> open(outputCsvFile, 'a'), "X,Y,Z"
outputLines = open(outputCsvFile, 'a')
lineCounter = 0
for inputLine in inputLines:
    if divmod(lineCounter,100)[1] == 0:
           print "Processing row " + str(lineCounter) + " of " + str(nRows) + "..."
    lineCounter = lineCounter + 1
    valueList = [float(value) for value in inputLine.split(" ") if value != '\n']
    columnCounter = 0
    for rasterValue in valueList:
        columnCounter = columnCounter + 1
        if rasterValue != nodataValue:
            xCoord = xCorner + (columnCounter - 1) * cellSize + cellSize / 2
            yCoord = yCorner + (nRows - lineCounter) * cellSize + cellSize / 2
            outputLines.write(str(xCoord) + "," + str(yCoord) + "," + str(rasterValue) + "\n")
outputLines.close()
inputLines.close()

Re: Extract x,y,z from Raster - script

Luke_Pinner — Tue, 12 Oct 2010 22:41:02 GMT

Another option, once GDAL 1.8.0 is released (http://www.gdal.org) is to use the GDAL command line tool "gdal_translate" (or the very useful python bindings, but that's a whole nuther ballgame) to convert to XYZ - http://www.gdal.org/frmt_xyz.html

If you want to try this out now, prebuilt binaries of the 1.8.0 development version are available at http://vbkto.dyndns.org/sdk/

The command line syntax would be something like:

gdal_translate -of XYZ inraster outxyz.txt

Re: Extract x,y,z from Raster - script

Luke_Pinner — Wed, 13 Oct 2010 21:43:18 GMT

..........

Re: Extract x,y,z from Raster - script

FrancesBiles1 — Fri, 25 Oct 2019 05:29:38 GMT

Hi Chris -

I was about to try your script, but my coordinates are longitude, latitude, elevation. Is your script able to compute the center of the cells correctly when dealing with degrees, or should this script only be used with coordinates in linear units?

Thanks,

Frances

Re: Extract x,y,z from Raster - script

ChrisSnyder — Mon, 28 Oct 2019 16:33:49 GMT

Francis -

In it's current state it will only work with planar coordinate systems (State Plane, UTM, etc.).

Another way to do what you need is to use the "Raster To Point" tool and just add/calc a Lat and Long field to the table. Since geometry is produced (in addition to tabular records) this method would take longer, but should yield you your necessary outcome.

Chris

Re: Extract x,y,z from Raster - script

FrancesBiles1 — Mon, 28 Oct 2019 23:49:21 GMT

Thanks for the clarification Chris.