EarthObservation / RVT_py

Relief Visualization Toolbox in Python
Apache License 2.0
25 stars 6 forks source link

Padding for SVF and SLOPE #12

Closed NejcCoz closed 1 year ago

NejcCoz commented 1 year ago

The output image of some visualizations is cropped compared to the input image. I have noticed this with SVF (openness, negative openness, probably also ASVF) and SLOPE visualizations (see image), I'm not sure if it happens with any others (SLRM and MSTP are ok). With SVF, the output is missing the same number of pixels as search radius, and 1 pixel edge with SLOPE. In the background is MSTP image, which has the same size as input.

To get those results I have used parameter values from defaults:

vis_out = {
   vis_type: default_1.get_slope(
   dict_arrays["array"],
   resolution_x=dict_arrays["resolution"][0],
   resolution_y=dict_arrays["resolution"][1]
)

vis_out = default_1.get_sky_view_factor(
   dict_arrays["array"],
   dict_arrays["resolution"][0],
   compute_svf=False,
   compute_opns=True
)
vis_out["neg_opns"] = default_1.get_neg_opns(
    dict_arrays["array"],
    dict_arrays["resolution"][0]
)

image

zm8597 commented 1 year ago

Hey :) Can you send me a minimal example script and small test DEM to reproduce. I have tried to reproduce this issue with my test cases but I couldn't.

NejcCoz commented 1 year ago

An example script and sample are in the ZIP.

I noticed something when I was preparing the example script... I am reading the array with rasterio and changing nodata values to np.nan before I pass them to RVT, could there be a problem?

tii.zip

NejcCoz commented 1 year ago

Hey Žiga. Today I had some time to look into the slope function. So, 1px padding is applied to the edges of the raster, and then removed at the end right before output. The problem describes above happens on the nodata edges of the image (not edges of an array). Could some sort of padding (dilation?) on the nodata edge solve this problem?

I assume the sfv function (and maybe some others) have the same problem?

zm8597 commented 1 year ago

Hey sorry for my horrible response time. I was very busy and forgot to check your example 🤦‍♂️. I just checked the example and I know exactly what the problem is. If we take a look at slope, it uses a 3x3 kernel so if one of the pixels inside this kernel is np.nan then the result will also be np.nan. We are using edge padding to remove edge effect (at edges kernel needs additional values), but this only works when the image is square (not rotated). The image in the example is slightly rotated (surrounded with nodata pixels, because array/raster is always square) so when we are trying to calculate out value of edge pixel (edge pixel of the image not the raster) one of the pixel of kernel is np.nan and causing also out pixel to be np.nan. The padding only adds pixels to the raster edge and not the image edge. To solve this issue we would need to interpolate neighbouring np.nan pixels of the image edge. I will take a look if any function can easily do that or if there is any easy algorithm to implement that.