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"
c:\temp\test_zxy.txt = sample(somegrid)
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...
# 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()
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
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
Thanks for the clarification Chris.
gdal_translate -of XYZ inraster outxyz.txt