Raster processing in ArcPy seems to "simplify"/"round" my resulting values

Question

I'm working off a DEM and doing some raster algebra with it (real simple, add this decimal value, then divide by this value). The DEM in question is super high resolution LiDAR, that's both spatially and in terms of vertical resolution. So we're working with even small changes in height being significant for this analysis (like centimeters, or hundredth of a meter considering the raster values are in meters).

The problem is the original DEM has very smooth value changes in this gently up-sloping terrain. The resulting raster though is jarring changes, almost like it's been reclassed. Just a big flat area of 0.3, then it jumps up to 1.5 (for example, this is not meters cause it's after the math). There's no gradient between these jumps in value.

I can't figure out what would be causing the issue. Any code that touches the DEM is below:

dem_extract = ExtractByMask(dem, extent_shp)
output_raster = (dem_extract - a) / b
output_raster.save(output_path)

b and a are float numbers, ~0-2.0 range (set by user). The output is Imagine img. Because of smaller a/b values, the output can range from like -1 to 100+. However if the large range is the problem (doesn't make sense as it's a double-value raster, sure the stretch value will change, but that's symbology dependent), I've also tried putting a raster con in addition to the extract to mask any very high elevations too.

Another possible issue is that the DEM is a HUGE mosaic raster (hence the extract first to work with less data and speed up the script). But I don't know why that would make Arc "round" the results and how I would fix that.

---

The problem is definitely in clipping/masking the raster. In just ArcMap I tried the Data Management version of Raster Clip, it came with the same "steps"/"bins" as the ExtractByMask result with breaks in the same locations, however, the pixel values were now just 1, 2, 3..

So the Extract By Mask tool, despite outputting decimal results, is using some trick (like int(value*100) then converting back with float(new_value)/100.0) but in the process is losing precision in a raster that has such high precision to begin with. Tried looking for a way to just force the mask/clip tools to work in floats but not seeing anything in the environment parameters.

My current stopgap solution is to convert the extent shapefile to a raster (with value just being any >0 field (careful with FID cause it starts at 0). Then use Raster Con to extract where the extent raster is > 0.

Answer 1

seems like a couple of issues going on.
python2 messed up with floats and ints are mixed. And ExtractByMask may be returning a raster with ints instead of floats.

perhaps this may help...the caveat is that I'm not an arcPy guy so you may need to fill in some blanks...and...this is all dependent on you having numpy accessible.

#get a numpy version of your dem
dem_np = arcpy.RasterToNumPyArray(dem, nodata_to_value=-9999)

#do your shape thing
dem_extract = ExtractByMask(dem, extent_shp)

#make numpy version
dem_extract_np = arcpy.RasterToNumPyArray(dem_extract, nodata_to_value=-9999)

#comment this out
#output_raster = (dem_extract - a) / b

#instead of relying on the values of dem_extract_np, which you think are suspect, lets use some numpy masking
#dem_extract_np != -9999  is a mask where there are values...so...

#make 'blank' raster
output_raster_np = np.ones(dem_extract.shape) * -9999

#assign values based on the mask
output_raster_np[dem_extract_np != -9999] = dem[dem_extract_np != -9999]

#do your math...phew! numpy is reliable!
output_raster_np = (output_raster_np - float(a)) / float(b)

#and one more operation...hmmmm...I wonder if this was also a problem in your process(unless arc magically does not operate on NODATAs? I have no idea.)
#reset some nodatas
output_raster_np[dem_extract_np == -9999] = -9999

#make arc raster
#you have values to fill in here your_value1 and your_value2
output_raster = arcpy.NumPyArrayToRaster(output_raster_np, lower_left_corner=your_value1, x_cell_size=your_value2, value_to_nodata=-9999)

output_raster.save(output_path)

EDIT: looking it over, I am assuming that ExtractByMask, when used with a non-rectangle polygon, will fill the areas within the output raster extent but not in the polygon with NODATA. It's a reasonable assumption but I could be wrong...lemme know!