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.