Raster to XYZ CSV

[Andrew_Needham]

Hi all,

I am trying to improve my python skills, has anyone got some tips to improve this code.

The purpose of this is to quickly as possible convert a single band raster (Elevation) into a CSV table that has the rasters x,y,z values.

Obviously this can also be done with various ESRI tools, but it can take hours for large areas, whereas this code will do it much quicker.

I'm still learning python, any tips would be appreciated.  Thanks

Also see attached for code

################################

import arcpy
from arcpy import env
import csv,os
import numpy

# Set local variables
inRaster = arcpy.GetParameterAsText(0)
csv_filepath = arcpy.GetParameterAsText(1)
roundValue = int(arcpy.GetParameterAsText(2))
createPointFC = arcpy.GetParameterAsText(3)
out_path = arcpy.GetParameterAsText(4)
out_name = arcpy.GetParameterAsText(5)

'''
# Set Fixed variables for testing
inRaster = r"D:\Input\Raster_clip.tif"
csv_filepath = r"D:\Output\Coordinates.csv"
roundValue = 2
createPointFC = "true" # false is default
out_path = r"D:\Output\Working.gdb"
out_name = "points"
'''

#variables
fld_names = ["X", "Y", "Z"] # For point file and CSV headers
noData = -100000000 #Value for no data to exclude later

# Execute Raster To NumPy Array for Dict of xyz values
arcpy.AddMessage("Creating raster dictionary...")
theRaster = inRaster
ras = arcpy.Raster(theRaster)
ar = arcpy.RasterToNumPyArray(theRaster, nodata_to_value=noData)
# Get Cellsize of raster
cellsize = (ras.meanCellHeight + ras.meanCellWidth) / 2
# Get XY extent
xmin = ras.extent.XMin
ymax = ras.extent.YMax
# Shift XY to center of cell
centerX = xmin + (cellsize/2)
centerY = ymax - (cellsize/2)
centerXj = centerX
# Amount of cells
yheight, xwidth = ar.shape
# Dictionary of (x,y): values
d={}

# Loop through X and Y for all cells and add X,Y,Z data to Dictionary
for j in range(yheight):#74
for i in range(xwidth):#84
value = ar.item(j, i)
if value != noData: # remove No Data values
rValue = round(value, roundValue)
rCenterX = round(centerX,roundValue)
rCenterY = round(centerY,roundValue)
d[rCenterX, rCenterY] = rValue
centerX = centerX + cellsize
centerY = centerY - cellsize
centerX = centerXj

# Create CSV and populate from Dictionary
arcpy.AddMessage("Creating CSV...")
my_dict = d
with open(csv_filepath, 'w') as csv_file:
writer = csv.writer(csv_file)
#writer.writerow(fld_names)
for key in my_dict.keys():
newKey = str(key).replace('(', '').replace(')', '')
csv_file.write("{}, {}\n".format(newKey,my_dict[key]))
csv_file.close()

# Create Point FC from Dictionary
if createPointFC == "true":
arcpy.AddMessage("Creating Point Featureclass...")
# Execute CreateFeatureclass
spatial_reference = arcpy.Describe(inRaster).spatialReference
ptFC = arcpy.CreateFeatureclass_management(out_path, out_name, geometry_type='POINT', has_z='ENABLED', spatial_reference=spatial_reference)
arcpy.management.AddField(ptFC, "X", "Float")
arcpy.management.AddField(ptFC, "Y", "Float")
arcpy.management.AddField(ptFC, "Z", "Float")
#fc = out_path &"\"& out_name
#fc = r"D:\Output\Working.gdb\Points"
fullPath = os.path.join(out_path, out_name)
with arcpy.da.InsertCursor(fullPath,['SHAPE@X', 'SHAPE@Y', 'SHAPE@Z'] + fld_names ) as c: #+ fld_names
for x,y in my_dict.keys():
z = my_dict[x,y]
shapeXYZ = x,y,z
c.insertRow((x,y,z,x,y,z))
del c

[DanPatterson]

Not sure what parts of the script is doing, so just an example for consideration.

Starting with a 2d numpy array and convert it to 2 forms of xyz arrays.  The last, a structured array with field names, makes it easier to convert to a csv.

    a = np.random.randint(0, 10, size=(4, 3))
    a_st = []
    idx = np.indices(a.shape)
    x = idx[0].flatten()
    y = idx[1].flatten()
    vals = a.flatten()
    xy = list(zip(x, y, vals))
    dt = [('X', '<f8'), ('Y', '<f8'), ('Val', '<i4')]
    xy = np.asarray(xy, dtype=dt)
    print("syntax... a\n{!r:}\nxy ...\n{}".format(a, xy))
    return a, xy

# example

a
array([[7, 2, 7],
       [1, 2, 7],
       [3, 4, 3],
       [9, 7, 3]])
# -- xy ...
[(0., 0., 7) (0., 1., 2) (0., 2., 7) (1., 0., 1) (1., 1., 2) (1., 2., 7)
 (2., 0., 3) (2., 1., 4) (2., 2., 3) (3., 0., 9) (3., 1., 7) (3., 2., 3)]

# -- as file for NumPyArrayToTable
(array([[7, 2, 7],
        [1, 2, 7],
        [3, 4, 3],
        [9, 7, 3]]),

XY
 array([(0., 0., 7), (0., 1., 2), (0., 2., 7), (1., 0., 1), (1., 1., 2),
        (1., 2., 7), (2., 0., 3), (2., 1., 4), (2., 2., 3), (3., 0., 9),
        (3., 1., 7), (3., 2., 3)],
       dtype=[('X', '<f8'), ('Y', '<f8'), ('Val', '<i4')]))

And my helper script which I use to bring back arrays as featureclasses (points with x, y and z values) or out to csv.

def save_txt(a, name="arr.txt", sep=", ", dt_hdr=True):
    """Save a NumPy structured/recarray to text.

    Parameters
    ----------
    a : array
        Input array.
    name : filename
        Output filename and path otherwise save to script folder.
    sep : separator
        Column separator, include a space if needed.
    dt_hdr : boolean
        If True, add dtype names to the header of the file.

    """
    a_names = ", ".join(a.dtype.names)
    hdr = ["", a_names][dt_hdr]  # use "" or names from input array
    s = np.array(a.tolist(), dtype=np.unicode_)
    widths = [max([len(i) for i in s[:, j]])
              for j in range(s.shape[1])]
    frmt = sep.join(["%{}s".format(i) for i in widths])
    np.savetxt(name, a, fmt=frmt, header=hdr, comments="")
    print("\nFile saved...")

Result

X, Y, Val
0.0, 0.0, 6
0.0, 1.0, 3
0.0, 2.0, 0
1.0, 0.0, 2
1.0, 1.0, 5
1.0, 2.0, 3
2.0, 0.0, 3
2.0, 1.0, 9
2.0, 2.0, 4
3.0, 0.0, 7
3.0, 1.0, 2
3.0, 2.0, 4

Hope you find a few tidbits therein.

glad to see you using numpy.