Hi, I'm attempting to run a simple canopy height model (CHM) using a surface and bare-earth terrain model. I figured it would be quicker to do this in numpy, so I am using the RasterToNumPyArray methods in arcpy:
import arcpy
import os
arcpy.CheckOutExtensions('Spatial')
arcpy.env.overwriteOuptut = True
datadir = 'D:/Data'
dsm = os.path.join(datadir, 'dsm.tif')
dtm = os.path.join(datadir, 'dtm.tif')
arcpy.env.outputCoordinateSystem = dsm
arcpy.env.cellSize = dsm
dsm_arr = arcpy.RasterToNumPyArray(dsm, nodata_to_value=0)
dtm_arr = arcpy.RasterToNumPyArray(dtm, nodata_to_value=0)
diff_arr = dsm_arr - dtm_arr
# Set things below 2m difference to 0
diff_arr[diff_arr < 2] = 0
diff = arcpy.NumPyArrayToRaster(diff_arr, value_to_nodata=0)
diff.save(os.path.join(datadir, 'diff.tif'))
arcpy.CheckInExtension('Spatial')
The result is that the diff.tif output is way off. The cell sizes are right, the SRS remains the same, and relatively speaking the values are correct (it's dropping negligible values I don't care about), but the georeferencing is off by thousands of miles. What am I doing wrong? Do I need to explicitly feed in a corner or extent? If so, how would I do that?