Specified z-factor being overwritten

[jonschwenk]

I have a DEM in EPSG:4326 with elevation units in feet. I am trying to use the FeaturePreservingSmoothing tool, but with a specified z-factor that accounts for my CRS and units. I am using the Python API.

I compute my z-factor from the central latitude of the DEM as

z = 1/(365130*np.cos(midlat*np.pi/180))

which in my case gives z = 0.00000739

When I specify this z factor in the feature_preserving_smoothing() call:

wbt.feature_preserving_smoothing(dem_path, denoised_path, zfactor=z)

I get the following message:
It appears that the DEM is in geographic coordinates. The z-factor has been updated to 0.000008833933.

I am not sure how WBT is computing z-factor, so I would prefer it to use the z-factor I specified rather than computing its own.

[jblindsay]

This is how WhiteboxTools calculates z-factor. I don't recall now where I came upon that particular equation, but I suspect it was an Esri blog at some point. It's also true that this tool, and all others that accept a z-factor parameter, will override the user-specified value in favour of this calculated one if the CRS is in geographic coordinates. This is somewhat opaque and I should probably explore alternatives to this.

[jonschwenk]

Ok, thanks for that.

WBT z factor is thus computed as z_factor = (1.0 / (113200.0 * mid_lat.cos())). A couple of things here:

1. I might be missing something but I believe the 1113200 is the number of meters per degree at the equator. Except this number should be 111320 (off by a factor of 10). [edit: I mentally appended an extra "1" to the z-factor--my mistake, this is not an issue.]

2. The z-factor is thus in units of m^-1. However, some DEMs don't use meters for the vertical unit (mine uses feet). Thus the constant in the z-factor equation should be the number of feet per degree at the equator (365130). I know I can just convert my DEM to meters, but this could be a gotcha.

From my use case and general perspective, I would prefer WBT to never overwrite a user-input value unless it leads to an error or instability. It probably doesn't matter too much for the denoising, but it absolutely does for slope, aspect etc. However, looking through the slope code, it appears that the same thing is happening, except without a warning to the user.

I'll just use gdal for now to compute these things instead. WBT has been very useful to me though, so thanks for developing it!

[jblindsay]

Although I don't think that this is the same source that I would have originally derived the equation from (it would have been many years ago now), you can find the same equation in Esri's help here. You'll see that they're using 113200 as well. You're right, that this would assume that the units are in meters rather than feet. Obviously, this is only an issue in the US, since outside of the US no one will have a DEM with units of feet...have you folks considered joining the rest of the world in that regard? ;-) Anyhow, I think you're right, WBT should not override a user-specified parameter. I will aim to modify all of the tools that have a z-conversion factor input parameter to favour a user-specified parameter over the calculated one. This will take a bit of time but it is a worthwhile change for certain. Then you will be able to use your feet-unit DEM.

[jonschwenk]

Ah, yes, your z-factor constant is fine. My mind was adding an additional "1" to it, making it off by a factor of 10.

I didn't make the DEM, I'm just using it!

[jblindsay]

Actually, I'm not convinced that my constant is correct. The number of meters in one degree at the equator is 111.32 km or 111320 m. That link that I gave previously used the value 113200 and I'm sure that I've seen that number used in other places on the Esri site, but I am beginning to question whether it's correct. Once we've gotten to the bottom of this, I'll make a change of all affected tools and make a point-release of WhiteboxTools within the next week. This is certainly an issue worthy of an update for. By the way, I'm just kidding with you about the feet thing. We love too joke about such things!

Edit:

Looks like this isn't the first time that I've been suspicious of the Esri number. See here.

[jonschwenk]

Oh, now that makes sense why I was trying to prepend another "1" to your constant. I usually use the 111km value, and the (meters) value I was using for the constant is 111320. There's some discussion here about what a difference the constant makes (it doesn't seem to be too much).

Haha, no worries. I am frustrated by our silly units system all the time.

[jblindsay]

Okay, I've just committed a large number of changes, all of which are related to this issue. The constant value has been updated to 111320 throughout the codebase (which I believe is the correct number of meters in a degree at the equator), and tools that take a z_conversion factor input parameter no longer override the user's value. This should resolve this issue. Please let me know if you continue to experience problems. I will be releasing version 1.4 shortly which will include these changes.