fmaussion
commented
4 years ago Thanks a lot for the report! We discussed this together a while ago, and now I had time to think about it more.
If we want the slope to be described in degree, we would need np.tan() ?!
After some more thought, I don't think so. Since tan(alpha) = y/x, alpha = arctan(y/x).
To stay consistent with the other slope definitions, it might be best to calculate here aswell the slope ratio (height difference against length change).
Yes that's true. After a quick search, we do not seem to be consistent everywhere here yet: I found 12 occurrences or arctan in the codebase. I wonder what's the "correct thing to do" here. I think there are good arguments to use the gradient only for the numerics (although I'll have to get back to Cuffey and Patterson to be sure), but then what about the other cases?
Currently, we use the "atan slope definition" in these places:
- in
centerlines.py, we use it to filter the centerlines according to a threshold parameter which the user defines in degrees in the parameter config file. I think this is pretty straightforward and uncontroversial - in
test_numerics.py, where I used the atan definition to choose a proper default value for the parameter above -> if we change the definition above, we should also change and rerun this test - in
gis.py, where we compute geographical attributes diagnostics of 2d glaciers, like slope, aspect... These values are not used in the model, they are here to convey statistics about the glaciers. Here again, I wonder if we shouldn't keep the atan definition because is is more correct - in
workflow.py, to give the slope statistics to users - this is the point raised by @lilianschuster (again, not used by the model)
So, to conclude, we use the gradient definition for all things related to numerics (bed inversion, forward dynamical model), but the atan definition for diagnostics and statistics.
Now the question is - should we stick to the gradient definition everywhere to be consistent?
In _workflow.py, in the method climate_statistics(), the slope is defined as an arctan of the height difference and length change.
@entity_task(log)def climate_statistics(gdir, add_climate_period=1995):slope = np.append(slope, np.arctan(-np.gradient(hgt, dx)))If we want the slope to be described in degree, we would need np.tan() ?!
To stay consistent with the other slope definitions, it might be best to calculate here aswell the slope ratio (height difference against length change). Suggested code change:
slope = np.append(slope, -np.gradient(hgt, dx))???