Skew method limitations and areas of improvement

[brentwilder]

Issue

There are tradeoffs between skew methods. Differences impact the horizon and therefore the net radiation via the terrain shadows. This minimally impacts sky view factor as well, although less so as n-angles increases.

Skew overview

The method skew() is called to rotate the image in order to compute horizons along a line for numpy array.

ARS Code

In the original ARS code, they defined it in this way,

def skew(arr, angle, fwd=True, fill_min=True):

    if angle == 0:
        return arr

    if angle > 45 or angle < -45:
        raise ValueError('skew angle must be between -45 and 45 degrees')

    nlines, nsamps = arr.shape

    if angle >= 0.0:
        negflag = False
    else:
        negflag = True
        angle = -angle

    slope = np.tan(angle * np.pi / 180.0)
    max_skew = int((nlines - 1) * slope + 0.5)

    o_nsamps = nsamps
    if fwd:
        o_nsamps += max_skew
    else:
        o_nsamps -= max_skew

    b = np.zeros((nlines, o_nsamps))
    if fill_min:
        b += np.min(arr)


    # NOTE , the critical piece for this issue is here:
    for line in range(nlines):
        o = line if negflag else nlines - line - 1
        offset = int(o * slope + 0.5)

        if fwd:
            b[line, offset:offset+nsamps] = arr[line, :]
        else:
            b[line, :] = arr[line, offset:offset+o_nsamps]

    return b

The offset is snapped via int, and so the result is that you end up with many streaking/lines in the data. However, the benefit of this approach, for there is no streaking is the values are exact ( no interpolation).

Interpolation approach

In the more recent implementation created in, #12 and #14 , the last part of the method was changed to perform linear interpolation.

   # same code as above ....
    line_indices = np.arange(nlines)
    o_indices = line_indices if negflag else nlines - line_indices - 1
    offsets = o_indices * slope

    output_grid = np.arange(o_nsamps)
    source_grid = np.arange(nsamps)

    if fwd:
        for i in range(nlines):
            offset = offsets[i]
            source = output_grid - offset
            mask = (source >= 0) & (source <= nsamps - 1)
            y_k = np.interp(source[mask], source_grid, arr[i, :])
            b[i, mask] = y_k

    else:
        for i in range(nlines):
            offset = offsets[i]
            output = source_grid - offset
            mask = (output >= 0) & (output <= o_nsamps - 1)
            y_k = np.interp(output_grid, output[mask], arr[i, source_grid[mask]])
            b[i, :] = y_k

    return b

However, this interpolation error manifests into a small bias for mostly shaded areas in early morning and late evening. This is because the horizon angles are smoothed out slightly. Which creates a small positive bias in net radiation.

Proposed solution?

• Perhaps better interpolation strategy (cubic, or other?)

• scipy.ndimage also has some tools for this but would require a new dependency. But could be worth testing.

[jomey]

Capturing another interpolation idea that @brentwilder and I briefly chatted about yesterday.

We could keep the current linear interpolation, but only apply the interpolated values to pixels where we pass a threshold difference.
Makeshift example for one line of pixels:

orig_array = [11,20,33,5,20]
interp_array = [11,20,33,18,20]
# Determine where values have a large difference
interp_mask = [0,0,0,1,0]
# Apply only to the pixels with threshold difference
new_array = [11,20,33,18,20]

[brentwilder]

@jomey So I did try this suggestion above, don't have any figures to show, but I was unable to get this to work. I was finding that my implementation of this, the thresholding was confusing ridgelines for striping. Was not able to overcome this in my iterations.

Although, I realized that we actually already install scipy within our iSnobal .yaml and so this would enable us to test these scipy.ndimage smoothing techniques for free. I was able to get median_filter to take care of the striping with just a little bit of tweaking. I will want to think about ways to make it more independent for any scene, although I think it is already an improvement.

[jomey]

Nice work finding that filter to address the striping and get back to the original angles for unaffected areas. I rather use a 'stock' method over some custom logic that we come up with.