Rotating and Translating a Vector Field

Question

I have an explicit expression for a vector field with components ($u',v'$) in some local coordinate frame ($x',y'$). The picture on the left below shows an example, where the vector field is constant in some direction inside a rectangle, and is zero everywhere outside the rectangle. Now I would like to calculate what the field would be in some laboratory frame of reference ($x,y$), in which I have positioned the object somewhere else (a translation by ($x_0,y_0$) and rotation by angle $\theta$ counterclockwise), as shown on the right:

It seems from trying to work it out on paper that I need to do the following steps:

1. First take my lab coordinates ($x,y$) and see how they are represented in the object's local frame ($x',y'$). This is done by applying a rotation matrix \begin{align} \begin{pmatrix} x' \\ y' \end{pmatrix} &= \begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} \end{align}

(note that confusingly to me this actually seems to be the transpose (or inverse?) of the rotation matrix for a counterclockwise rotation by angle $\theta$.)

2. Then use these coordinates to calculate my field components using my analytic expression to get ($u'(x',y'),v'(x',y')$).

3. Finally, it seems I need to transform the resulting vectors by applying a rotation matrix again (the inverse of the original one this time): \begin{align} \begin{pmatrix} u \\ v \end{pmatrix} &= \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix} \begin{pmatrix} u' \\ v' \end{pmatrix} \end{align}

This is somehow "projecting" the first components back into the lab frame?

My questions are basically: Is this line of thought correct? Is it expected that I should need to apply rotation matrix twice when transforming a vector field - one being the inverse of the other?

---

Numerical example

Here is some Matlab code which shows that the steps above seem to work:

x0 = 3;     % Translation offsets
y0 = 4;
theta = 40; % Rotation counterclockwise around z axis (in degrees)

x = linspace(-10,10,1000); y = x;[X,Y] = meshgrid(x,y); % Laboratory frame grid

[U,V] = f(X,Y,x0,y0,theta);  % Calculate vector function

imagesc(x,y,sqrt(U.^2+V.^2)) % Plot magnitude of vector function
quiver(X,Y,U,V)              % Plot arrows of vector function

function [U,V] = f(X,Y,x0,y0,theta)
    [Xr,Yr] = transform_coordinates(X,Y,x0,y0,theta);
    [Ur,Vr] = calculate_components(Xr,Yr);
    [U,V] = project_components_back(Ur,Vr,theta);
end

function [Xr,Yr] = transform_coordinates(X,Y,x0,y0,theta)
    X = X-x0;
    Y = Y-y0;
    R = [+cosd(theta) +sind(theta); 
         -sind(theta) +cosd(theta)]; 
    rotXY = [X(:) Y(:)]*R.';
    Xr = reshape(rotXY(:,1), size(X,1), []);
    Yr = reshape(rotXY(:,2), size(Y,1), []);
end

function [U,V] = project_components_back(Ur,Vr,theta)
    R = [+cosd(theta) -sind(theta); 
         +sind(theta) +cosd(theta)]; 
    rotUV = [Ur(:) Vr(:)]*R.';
    U = reshape(rotUV(:,1), size(Ur,1), []);
    V = reshape(rotUV(:,2), size(Vr,1), []);
end

function [Ur,Vr] = calculate_components(Xr,Yr)
    % Some operation which represents my function (in this case, create a rectangle)
    lx = 3; ly = 5; % Rectangle side lengths
    Ur = zeros(size(Xr));
    Vr = Ur;
    Ur(abs(Xr)<lx/2 & abs(Yr)<ly/2) = 0;
    Vr(abs(Xr)<lx/2 & abs(Yr)<ly/2) = 1; % Vector field points along rectangle axis
end