How to create a 3D point grid?

[David_Brooks]

I need to create a 3D point grid such as the screenshot below, that is defined by a Max/Min of X,Y and Z and  with a defined point spacing.

My first thought was to create a raster with the xy and spacing/pixel size to match the variables above, then convert to point. Or to use the fishnet tool to extract the "label point" output. But the issue is how to duplicate in the vertical, which may need duplicating up to 200x depending on the scenario.

This might be easier to achieve in python, but am unsure of how to go about it. The XYZ min/max would all be defined by a collection of terrain rasters. The pixel/point spacing would just be a manual input.

Unless there's a better way, my thinking is to duplicate rasters with a fixed value from Zmin to Zmax with increments matching the point spacing, then create 3D points from these.

[JohannesLindner]

A very basic implementation:

wkid = 25832  # coordinate system
x_range = [423000, 423100, 10]  # start, stop, step
y_range = [5970000, 5970100, 10]
z_range = [1, 10, 1]

points_3d = arcpy.management.CreateFeatureclass("memory", "points_3d", "POINT", has_z="YES", spatial_reference=wkid)
with arcpy.da.InsertCursor(points_3d, ["SHAPE@"]) as cursor:
    for x in range(*x_range):
        for y in range(*y_range):
            for z in range(*z_range):
                p = arcpy.Point(x, y, z)
                g = arcpy.PointGeometry(p, wkid)
                cursor.insertRow([g])

[David_Brooks]

@JohannesLindner thanks a lot for the script. I gave it a whirl, but it took 1hr 40mins to create a point array at the size I needed. However, I persisted with my previous approach (creating rasters) and it's only taken 10mins for the same area. Here's the interim script:

import arcpy
from arcpy import env
from arcpy.sa import *
inputPolygon = arcpy.GetParameterAsText(0)
outputPoint = arcpy.GetParameterAsText(1)
cellsize = arcpy.GetParameterAsText(2)
maxDepth = arcpy.GetParameterAsText(3)
arcpy.env.workspace = r"INSERT"
#create base raster
rasterBase = "memory\\rasterBase"
arcpy.conversion.PolygonToRaster(inputPolygon, 'OBJECTID', rasterBase, '', '', cellsize, 'TRUE')
#reclassify base raster to match depth values and convert to point2D then point3D
depthRange = list(range(0,int(maxDepth)+1))
outRaster = "memory\outRaster"
for depth in depthRange:    
    MyRemapRange = RemapRange([[0,9999,depth]])
    outRas = Reclassify(rasterBase, 'Value', MyRemapRange, 'NODATA')
    outRas.save(outRaster)
    outPoint2D = "memory\outPoint2D"
    arcpy.conversion.RasterToPoint(outRaster, outPoint2D, 'Value')
    outPoint3D = outputPoint+"_"+str(depth)+"m"
    arcpy.ddd.FeatureTo3DByAttribute(outPoint2D, outPoint3D, 'grid_code')
    arcpy.AddMessage("Saved 3D Point_"+str(depth)+"m")
arcpy.AddMessage("Finished...")

with the added advantage that I can defined the extent of the point array with a polygon, and control the depth and resolution via parameters.

[Robert_LeClair]

In a way this looks reminiscent of a Space Time Cube but the creation of those is from a varierty of data sources.  My guess is you're wanting to create the empty structure perhaps or not so much?

[David_Brooks]

@Robert_LeClair BINGO! I'm working with a geophysical team on an offshore project. We have a ground model in Kingdom software, and they can only provide me with isopach rasters (base of soil unit rasters) so I plan to create an array of points, then loop through each point and assign a soil type based on the existence of an isopach below it. It's the only way  to create a voxel ground model from isopachs as far as I can see. 

Have tried using "empirical bayesian kriging 3D", but you don't want to interpolate data when you have isosurfaces, as the data is absolute.

[David_Brooks]

Next challenge is working out how to get regularly gridded 3D point data into a compatible voxel format without using the Space Time Cube or Bayesian Kriging tools. Currently failing 🤔