nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

Extract x,y,z from Raster - script

JohnLonsdale

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

Finde mehr Beiträge mit dem Tag

Accepted answers

All comments

FrancesBiles1

Thanks for the clarification Chris.

ChrisSnyder

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

FrancesBiles1

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

Luke_Pinner

..........

Luke_Pinner

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

ChrisSnyder

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

JohnLonsdale

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

Luke_Pinner

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

ChrisSnyder

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.

Nimm an der Unterhaltung teil

Angemeldete Mitglieder können Beiträge verfassen, Updates folgen und mehr. Neu hier? Registriere ein kostenloses Konto.

Schnellzugriffe

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.