The objective is to turn a raster surface (double value type) into smooth, visually appealing polygons representing iso-levels. I cannot stress enough that these must be visually attractive. Jaggedy lines that are raster remnants are not acceptable. Here are the parameters of my problem:
- Input: Raster surface of atmospheric forecasts (not a DEM). Z values can fluctuate greatly over a short horizontal distance, but the surface does tend to have sane looking clusters of pixels with peaks and valleys.
- Desired output: Smooth polygons representing iso-levels at irregular intervals. Class breaks are 0-5, 5-20, 20-75, 75-150, and >150.
- I really have artistic license to make these look nice. I have been requested to remove all remnants of jagged raster edges. Obviously, I want to be as true as possible to the underlying data, but the marching orders are something that is simple and easy for the novice user to understand without a lot of details distracting the eye.
- Needs to be automated to run on regular intervals. I can code in several languages, but prefer python or command line utilities.
- I have access to ArcGIS or can use FOSS tools, but prefer to run in Linux, which probably favors the latter.
I'll talk a bit about what I've tried below, but everything I try runs into problem after problem, to the point where it's just getting too complicated. My inner mental dialog is saying "It's 2017, somebody must have really solved this problem by now!" So... am I missing the obvious solution?
First, I tried reclassifying the raster into categories, and vectorizing that to polys. I got similar results from Esri "reclass by ASCII file" and gdal_calc. The big problem here is that it produced lots of areas where there adjacent polygons skipped a category or meet in diagonal corners (I've flagged a few in the diagram), and I'd really like concentric bands. I might be able to reduce the cell size and smooth this enough to solve the problem, but I'm not sure how far down this path I want to go... Even with the polygon simplification in ArcGIS, I still had to smooth the lines to make them look halfway decent, but this produced some strange topology issues with holes in my coverage. I'm sure I could solve this, too, but I feel like I'm just going down multiple rabbit trails at this point.
The other thing I tried was to create contours. Esri's "contour list" tool produced slightly nicer results than the gdal_contour command, but I was able to produce contour lines using both. Ultimately, those were just lines, not polygons. I found a few references to converting to polys, but this feels like I'm re-inventing the wheel, so I wanted to come up for air and see if this is on the right track. (Esri toolbox produced error when run, blog post not automated yet)
There is also the Esri "GA Layer to Contour" tool, but I was having trouble figuring out the inputs to that. It was unclear to me how to create the input layer needed for that tool and the documentation was less than clear and lacking useful examples. It looks like I need to convert my raster to points and interpolate it again? Any hints on this one?
I feel like I have a number of really complicated options, but this seems like a problem that should have been solved a long time ago. I also feel like I have a lot of useful tools, but I'm not putting them together in the right order...