Rotating a 2d sub-array using numpy without aliasing effects

Question

I would like to rotate only the positive value pixels in my 2d array some degree about the center point. The data represents aerosol concentrations from a plume dispersion model, and the chimney position is the origin of rotation.

I would like to rotate this dispersion pattern given the wind direction.

The concentrations are first calculated for a wind direction along the x-axis and then translated to their rotated position using a 2d linear rotation about the center point of my array (the chimney position) for all points whose concentration is > 0. The input X,Y to the rotation formula are pixel indexes.

My problem is that the output is aliased since integers become floats. In order to obtain integers, I rounded up or down the output. However, this creates null cells which become increasingly numerous as the angle increases.

Can anyone help me find a solution to my problem? I would like to fix this problem if possible using numpy, or a minimum of packages...

The part of my script that deals with computing the concentrations and rotating the pixel by 50°N is the following. Thank you for your help.

def linear2D_rotation(xcoord,ycoord,azimuth_degrees):
    radians = (90 - azimuth_degrees) * (np.pi / 180) # in radians
    xcoord_rotated = (xcoord * np.cos(radians)) - (ycoord * np.sin(radians))
    ycoord_rotated = (xcoord * np.sin(radians)) + (ycoord * np.cos(radians))
    return xcoord_rotated,ycoord_rotated

u_orient = 50 # wind orientation in degres from North
kernel = np.zeros((NpixelY, NpixelX))  # initialize matrix
Yc = int((NpixelY - 1) / 2)  # position of central pixel
Xc = int((NpixelX - 1) / 2)  # position of central pixel

nk = 0
for Y in list(range(0,NpixelX)):
    for X in list(range(0,NpixelY)):
        # compute concentrations only in positive x-direction
        if (X-Xc)>0:
            # nnumber of pixels to origin point (chimney)
            dx = ((X-Xc)+1)
        dy = ((Y-Yc)+1)
        # distance of point to origin (chimney)
        DX = dx*pixel_size_X
        DY = dy*pixel_size_Y
        # compute diffusivity coefficients
        Sy, Sz = calcul_diffusivity_coeff(DX, stability_class)
        # concentration at ground level below the centerline of the plume
        C = (Q / (2 * np.pi * u * Sy * Sz)) * \
            np.exp(-(DY / (2 * Sy)) ** 2) * \
            (np.exp(-((Z - H) / (2 * Sz)) ** 2) + np.exp(-((Z + H) / (2 * Sz)) ** 2))  # at point away from center line
        C = C * 1e9  # convert MBq to Bq

        # rotate only if concentration value at pixel is positive
        if C > 1e-12:
            X_rot, Y_rot = linear2D_rotation(xcoord=dx, ycoord=dy,azimuth_degrees=u_orient)
            X2 = int(round(Xc+X_rot))
            Y2 = int(round(Yc-Y_rot)) # Y increases downwards
            # pixels that fall out of bounds -> ignore
            if (X2 > (NpixelX - 1)) or (X2 < 0) or (Y2 > (NpixelY - 1)):
                continue
            else:
                # replace new pixel position in kernel array
                kernel[Y2, X2] = C

The original array to be rotated

The rotated array by 40°N showing the data loss

Answer 1

Your problem description is not 100% clear, but here are a few recommendations:

1.) Don't reinvent the wheel. There are standard solutions for things like rotating pixels. Use them! In this case

• scipy.ndimage.affine_transform for performing the rotation
• a homogeneous coordinate matrix for specifying the rotation
• nearest neighbor interpolation (parameter order=0 in code below).

2.) Don't loop where not necessary. The speed you gain by not processing non-positive pixels is nothing against the speed you lose by looping. Compiled functions can ferry around a lot of redundant zeros before hand-written python code catches up with them.

3.) Don't expect a solution that maps pixels one-to-one because it is a fact that there will be points that are no ones nearest neighbor and points that are nearest neighbor to multiple other points. With that in mind, you may want to consider a higher order, smoother interpolation.

Comparing your solution to the standard tools solution we find that the latter

gives a comparable result much faster and without those hole artifacts.

Code (without plotting). Please note that I had to transpose and flipud to align the results :

import numpy as np
from scipy import ndimage as sim
from scipy import stats

def mock_data(n, Theta=50, put_neg=True):
    y, x = np.ogrid[-20:20:1j*n, -9:3:1j*n, ]
    raster = stats.norm.pdf(y)*stats.norm.pdf(x)
    if put_neg:
        y, x = np.ogrid[-5:5:1j*n, -3:9:1j*n, ]
        raster -= stats.norm.pdf(y)*stats.norm.pdf(x)
        raster -= (stats.norm.pdf(y)*stats.norm.pdf(x)).T
    return {'C': raster * 1e-9, 'Theta': Theta}

def rotmat(Theta, offset=None):
    theta = np.radians(Theta)
    c, s = np.cos(theta), np.sin(theta)
    if offset is None:
        return np.array([[c, -s] [s, c]])
    R = np.array([[c, -s, 0], [s, c,0], [0,0,1]])
    to, fro = np.identity(3), np.identity(3)
    offset = np.asanyarray(offset)
    to[:2, 2] = offset
    fro[:2, 2] = -offset
    return to @ R @ fro

def f_pp(C, Theta):
    m, n = C.shape
    clipped = np.maximum(0, 1e9 * data['C'])
    clipped[:, :n//2] = 0
    M = rotmat(Theta, ((m-1)/2, (n-1)/2))
    return sim.affine_transform(clipped, M, order = 0)

def linear2D_rotation(xcoord,ycoord,azimuth_degrees):
    radians = (90 - azimuth_degrees) * (np.pi / 180) # in radians
    xcoord_rotated = (xcoord * np.cos(radians)) - (ycoord * np.sin(radians))
    ycoord_rotated = (xcoord * np.sin(radians)) + (ycoord * np.cos(radians))
    return xcoord_rotated,ycoord_rotated

def f_OP(C, Theta):
    kernel = np.zeros_like(C)
    m, n = C.shape
    for Y in range(m):
        for X in range(n):
            if X > n//2:
                c = C[Y, X] * 1e9
                if c > 1e-12:
                    dx = X - n//2 + 1
                    dy = Y - m//2 + 1
                    X_rot, Y_rot = linear2D_rotation(xcoord=dx, ycoord=dy,azimuth_degrees=Theta)
                    X2 = int(round(n//2+X_rot))
                    Y2 = int(round(m//2-Y_rot)) # Y increases downwards
                    # pixels that fall out of bounds -> ignore
                    if (X2 > (n - 1)) or (X2 < 0) or (Y2 > (m - 1)):
                        continue
                    else:
                        # replace new pixel position in kernel array
                        kernel[Y2, X2] = c
    return kernel

n = 100
data = mock_data(n, 70)