diff --git a/.claude/skills/newop/SKILL.md b/.claude/skills/newop/SKILL.md
new file mode 100644
index 00000000..0f1ae2e5
--- /dev/null
+++ b/.claude/skills/newop/SKILL.md
@@ -0,0 +1,124 @@
+---
+name: newop
+description: Create a new PyLops linear operator following docs/source/adding.rst - class file, docstring, tests, docs entry and example. Use when the user asks to add/implement/port a new operator into PyLops, including porting an existing non-PyLops forward/adjoint implementation from a URL or a local file (e.g. "add a Foo operator", "turn this script into a PyLops operator", "port the operator at <link>").
+---
+
+Goal: add a new, PyLops-compliant `LinearOperator` to the library, complete with
+docstring, registration, tests, docs entry and a gallery example, following
+`docs/source/adding.rst` (the authoritative guide - read it if unsure).
+
+The operator may be written from scratch (from a mathematical description) or
+**ported** from an existing non-PyLops implementation supplied as a **web link**
+or a **local file**. Ask for the operator name and the source only if neither is
+inferable from the invocation.
+
+## 0. Get the source material
+
+- **Web link**: fetch it with `WebFetch` (or the browser tools if the page needs
+  JS). Extract the actual forward/adjoint code, not the prose.
+- **Local file**: read it in full.
+- **Neither**: work from the user's mathematical description, and state the
+  assumed definition of the operator before writing code.
+
+Then write down explicitly, before touching `pylops/`:
+
+- what the forward map does, and its input/output shapes;
+- whether the source's "adjoint" is a true adjoint (\(\mathbf{A}^H\)) or merely an
+  inverse/transpose/approximation - **this is the most common porting bug**;
+- which source parameters become `__init__` arguments, which become derived
+  members, and which are irrelevant (e.g. plotting, I/O, CLI args);
+- whether the operator is real- or complex-linear, and whether it is `explicit`.
+
+If the source adjoint is not the true adjoint, say so and implement the correct
+adjoint - the dot-test in step 4 will fail otherwise. Never relax the dot-test
+tolerance to make a wrong adjoint pass.
+
+## 1. Place the file
+
+- One class per file; file named after the class but **lowercase**
+  (`pylops/basicoperators/diagonal.py` holds `Diagonal`). Choose the subpackage by
+  theme: `basicoperators`, `signalprocessing`, `waveeqprocessing`, `optimization`,
+  etc. Create a new subpackage only if nothing fits.
+- If the operator is just a composition of existing operators, write a **function**
+  returning the composed operator instead of a class (see `pylops.Laplacian`).
+- Start the file with `__all__ = ["<Operator>"]`.
+- Register it: add the import/`__all__` entry in the subpackage `__init__.py`
+  (and its module-level summary table), plus the top-level `pylops/__init__.py`
+  if the operator is meant to be user-facing as `pylops.<Operator>`.
+
+## 2. Write the class
+
+Use `reference/operator_template.py` as the skeleton. Key rules:
+
+- Inherit from `pylops.LinearOperator` and initialize via
+  `super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name)`.
+  Prefer `dims`/`dimsd` over setting `shape` directly; `shape` is derived.
+  Set `explicit=True` only when the operator also exposes a dense matrix `A`.
+- Decorate `_matvec`/`_rmatvec` with `@reshaped` when the operator is
+  n-dimensional, so `x` arrives shaped as `dims` (`dimsd` for `_rmatvec`) and the
+  return value is flattened for you.
+- Use the backend helpers rather than raw NumPy so CuPy/JAX work:
+  `pylops.utils.backend.get_array_module`, `to_cupy_conditional`, and friends.
+  Do not `import numpy` for array creation inside `_matvec`/`_rmatvec`.
+- Type-annotate with `pylops.utils.typing` (`NDArray`, `DTypeLike`,
+  `InputDimsLike`).
+- Keep a `name` argument (default a short string) for `pylops.utils.describe`.
+- Write the `numpydoc` docstring with, at minimum: one-line summary, expanded
+  description, `Parameters`, `Attributes` (when non-obvious), `Raises` (when the
+  `__init__` validates inputs), and a `Notes` section giving the maths of forward
+  and adjoint in `.. math::` blocks. Match the level of detail of neighbouring
+  operators.
+
+## 3. Add tests
+
+Add to the existing `pytests/test_*.py` matching the subpackage, or create a new
+one following the same header (the `TEST_CUPY_PYLOPS` / `backend` guard block).
+Follow `reference/test_template.py`:
+
+- module-level `par*` dicts, parametrized with `@pytest.mark.parametrize("par", [...])`
+  covering real/complex and, where relevant, square/over-/under-determined;
+- an `assert dottest(Op, nr, nc, rtol=..., complexflag=0 if par["imag"] == 0 else 3, backend=backend)`
+  in every test of a new configuration;
+- a forward check against an independently computed expected result
+  (e.g. `Op.todense() @ x`, or the original source implementation's output);
+- an inversion round-trip with `lsqr` / `Op / y` and `assert_array_almost_equal`
+  when the operator is invertible;
+- error-path tests for anything the `__init__` raises.
+
+## 4. Run
+
+Always use `uv`:
+
+```bash
+uv run pytest pytests/test_<file>.py -k <Operator> -q
+make lint_uv
+```
+
+Iterate until the dot-test and all assertions pass cleanly.
+
+## 5. Document
+
+- Add the operator name to the right `autosummary` block in
+  `docs/source/api/index.rst`.
+- Add a gallery example `examples/plot_<operator>.py` (or a tutorial in
+  `tutorials/` for a heavier workflow), following the sphinx-gallery format of
+  `examples/plot_diagonal.py`: `r"""` title/underline/description `"""` header,
+  then `###...` comment blocks separating narrative from code, and matplotlib
+  figures showing forward and adjoint (and inversion, if relevant).
+
+## 6. Final checklist (from `docs/source/adding.rst`)
+
+Report back confirming each item:
+
+- [ ] single class (or function) in its own file, in a suitable `pylops` subpackage
+- [ ] `__init__`, `_matvec`, `_rmatvec` implemented (plus `todense`/`matrix` if cheap)
+- [ ] operator exported from the subpackage and top-level `__init__.py`
+- [ ] numpydoc docstring with `Parameters` and a mathematical `Notes` section
+- [ ] test added, `dottest` passes, forward/inverse checked
+- [ ] listed in `docs/source/api/index.rst`
+- [ ] used in at least one `examples/` script or `tutorials/` script
+- [ ] `make lint_uv` clean
+
+When porting, close with a short note on what differed between the source
+implementation and the PyLops version (adjoint correction, shape/flattening
+conventions, dtype handling, removed I/O).
diff --git a/.claude/skills/newop/reference/operator_template.py b/.claude/skills/newop/reference/operator_template.py
new file mode 100644
index 00000000..316e901a
--- /dev/null
+++ b/.claude/skills/newop/reference/operator_template.py
@@ -0,0 +1,97 @@
+"""Skeleton for a new PyLops operator.
+
+Copy into ``pylops/<subpackage>/<operatorname>.py``, rename, and fill in.
+Delete anything that does not apply.
+"""
+
+__all__ = ["MyOperator"]
+
+import numpy as np
+
+from pylops import LinearOperator
+from pylops.utils._internal import _value_or_sized_to_tuple
+from pylops.utils.backend import get_array_module, to_cupy_conditional
+from pylops.utils.decorators import reshaped
+from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray
+
+
+class MyOperator(LinearOperator):
+    r"""One-line summary of the operator.
+
+    Longer description of what the operator applies in forward mode and what
+    its adjoint does.
+
+    .. versionadded:: X.Y.Z
+
+    Parameters
+    ----------
+    param : :obj:`numpy.ndarray`
+        Description of the main parameter.
+    dims : :obj:`list` or :obj:`int`, optional
+        Number of samples for each dimension of the model.
+    axis : :obj:`int`, optional
+        Axis along which the operator is applied.
+    dtype : :obj:`str`, optional
+        Type of elements in input array.
+    name : :obj:`str`, optional
+        Name of operator (to be used by :func:`pylops.utils.describe.describe`)
+
+    Attributes
+    ----------
+    shape : :obj:`tuple`
+        Operator shape.
+    explicit : :obj:`bool`
+        Operator contains a matrix that can be solved explicitly (``True``) or
+        not (``False``).
+
+    Raises
+    ------
+    ValueError
+        If ``param`` has incompatible size with ``dims``.
+
+    Notes
+    -----
+    In forward mode the operator applies
+
+    .. math::
+        y_i = \ldots \quad \forall i=1,2,\ldots,N
+
+    and in adjoint mode
+
+    .. math::
+        x_i = \ldots \quad \forall i=1,2,\ldots,M
+
+    """
+
+    def __init__(
+        self,
+        param: NDArray,
+        dims: int | InputDimsLike | None = None,
+        axis: int = -1,
+        dtype: DTypeLike = "float64",
+        name: str = "M",
+    ) -> None:
+        self.param = param
+        self.axis = axis
+        dims = param.shape if dims is None else _value_or_sized_to_tuple(dims)
+        # dimsd is the shape of the data (output of the forward)
+        dimsd = dims
+        super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name)
+
+    @reshaped
+    def _matvec(self, x: NDArray) -> NDArray:
+        ncp = get_array_module(x)
+        if type(self.param) is not type(x):
+            self.param = to_cupy_conditional(x, self.param)
+        y = ncp.zeros(self.dimsd, dtype=self.dtype)
+        # ... forward implementation, y = A x
+        return y
+
+    @reshaped
+    def _rmatvec(self, y: NDArray) -> NDArray:
+        ncp = get_array_module(y)
+        if type(self.param) is not type(y):
+            self.param = to_cupy_conditional(y, self.param)
+        x = ncp.zeros(self.dims, dtype=self.dtype)
+        # ... adjoint implementation, x = A^H y (conjugate for complex params!)
+        return x
diff --git a/.claude/skills/newop/reference/test_template.py b/.claude/skills/newop/reference/test_template.py
new file mode 100644
index 00000000..d125f5b1
--- /dev/null
+++ b/.claude/skills/newop/reference/test_template.py
@@ -0,0 +1,57 @@
+"""Skeleton test for a new PyLops operator.
+
+Merge into the ``pytests/test_<subpackage>.py`` that matches the operator.
+Keep the CuPy guard header identical to the one already in that file.
+"""
+
+import os
+
+if int(os.environ.get("TEST_CUPY_PYLOPS", 0)):
+    import cupy as np
+    from cupy.testing import assert_array_almost_equal
+
+    backend = "cupy"
+else:
+    import numpy as np
+    from numpy.testing import assert_array_almost_equal
+
+    backend = "numpy"
+import pytest
+
+from pylops.basicoperators import MyOperator  # noqa: F401  (adjust import)
+from pylops.optimization.basic import lsqr
+from pylops.utils import dottest
+
+par1 = {"ny": 11, "nx": 11, "imag": 0, "dtype": "float64"}  # square real
+par2 = {"ny": 21, "nx": 11, "imag": 0, "dtype": "float64"}  # overdetermined real
+par1j = {"ny": 11, "nx": 11, "imag": 1j, "dtype": "complex128"}  # square complex
+par2j = {"ny": 21, "nx": 11, "imag": 1j, "dtype": "complex128"}  # overdet. complex
+
+
+@pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j)])
+def test_MyOperator(par):
+    """Dot-test, forward and inversion for MyOperator"""
+    param = np.arange(par["nx"]) + 1.0 + par["imag"] * (np.arange(par["nx"]) + 1.0)
+
+    Op = MyOperator(param, dtype=par["dtype"])
+    assert dottest(
+        Op,
+        par["ny"],
+        par["nx"],
+        rtol=1e-6 if par["dtype"] in ("float64", "complex128") else 1e-3,
+        complexflag=0 if par["imag"] == 0 else 3,
+        backend=backend,
+    )
+
+    x = np.ones(par["nx"]) + par["imag"] * np.ones(par["nx"])
+    y = Op * x
+    assert_array_almost_equal(y, Op.todense() @ x, decimal=6)
+
+    xinv = lsqr(Op, y, x0=np.zeros_like(x), niter=300, show=0)[0]
+    assert_array_almost_equal(x, xinv, decimal=4)
+
+
+def test_MyOperator_raises():
+    """Check input validation of MyOperator"""
+    with pytest.raises(ValueError):
+        MyOperator(np.ones(5), dims=(4,))
diff --git a/AIPOLICY.md b/AIPOLICY.md
index 04e00845..ab8f6420 100644
--- a/AIPOLICY.md
+++ b/AIPOLICY.md
@@ -36,5 +36,6 @@ More specifically, we currently provide:
 
 - ``.pi/prompts/optest.md`` / ``.claude/skills/optest``: a skill to increase the test coverage of an operator;
 
+- ``.claude/skills/newop``: a skill to create a new operator from a mathematical description or a plain implementation of forward and adjoint from file or URL;
 
 🤖🤖 **This Policy was written by humans and polished by AI** 🤖🤖
diff --git a/docs/source/api/index.rst b/docs/source/api/index.rst
index 8084d114..246ec9cb 100755
--- a/docs/source/api/index.rst
+++ b/docs/source/api/index.rst
@@ -103,6 +103,7 @@ Signal processing
     Interp
     InterpCubicSpline
     Bilinear
+    Downsample2D
     FFT
     FFT2D
     FFTND
diff --git a/examples/plot_downsample.py b/examples/plot_downsample.py
new file mode 100644
index 00000000..872f8cc1
--- /dev/null
+++ b/examples/plot_downsample.py
@@ -0,0 +1,140 @@
+r"""
+Downsampling
+============
+This example shows how to use the
+:py:class:`pylops.signalprocessing.Downsample2D` operator to reduce the size
+of a 2-dimensional array along both of its directions.
+
+Downsampling is performed in two steps: an anti-aliasing Gaussian filter is
+first applied to the input array, and the smoothed array is subsequently
+subsampled by the required decimation factors. Whilst a naive subsampling of
+the input array would fold any energy above the Nyquist wavenumber of the
+coarse grid back onto the retained wavenumbers (i.e., aliasing), the Gaussian
+filter removes such energy prior to decimation.
+
+As the operator is linear, its adjoint (and, more interestingly, its inverse)
+can also be used to move back from the coarse to the fine grid; the latter
+represents a very simple form of *super-resolution*.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy import datasets
+
+import pylops
+
+plt.close("all")
+np.random.seed(0)
+
+###############################################################################
+# Let's start by creating a 2-dimensional input vector containing an image
+# from the ``scipy.datasets`` family and downsample it by a factor of 4 in
+# both directions.
+x = datasets.face()[::2, ::2, 0].astype(np.float64)
+nz, nx = x.shape
+
+Dop = pylops.signalprocessing.Downsample2D((nz, nx), factors=4)
+y = Dop @ x
+
+print(Dop)
+print(f"Model size: {Dop.dims}, Data size: {Dop.dimsd}")
+
+fig, axs = plt.subplots(1, 2, figsize=(10, 4))
+axs[0].imshow(x, cmap="gray")
+axs[0].set_title(f"Original {Dop.dims}")
+axs[0].axis("tight")
+axs[1].imshow(y, cmap="gray")
+axs[1].set_title(f"Downsampled {Dop.dimsd}")
+axs[1].axis("tight")
+plt.tight_layout()
+
+###############################################################################
+# The role of the anti-aliasing Gaussian filter becomes evident if we take
+# the input and output in the frequency domain. Note how a simple resampling
+# of the input array (i.e., picking one every four samples in each direction)
+# would lead to aliasing of the high wavenumbers, which is instead not present
+# in the downsampled array.
+Fop = pylops.signalprocessing.FFT2D((nz, nx), fftshift_after=True)
+F1op = pylops.signalprocessing.FFT2D((nz // 4, nx // 4), fftshift_after=True)
+
+xf = Fop @ x
+yf = F1op @ y
+yf1 = F1op @ x[::4, ::4]
+
+fig, axs = plt.subplots(1, 3, figsize=(10, 4))
+axs[0].imshow(
+    np.abs(xf)[
+        nz // 2 - nz // 8 : nz // 2 + nz // 8, nx // 2 - nx // 8 : nx // 2 + nx // 8
+    ],
+    cmap="jet",
+    vmin=0,
+    vmax=0.005 * np.abs(xf).max(),
+)
+axs[0].set_title("Original (centered)")
+axs[0].axis("tight")
+axs[1].imshow(np.abs(yf), cmap="jet", vmin=0, vmax=0.005 * np.abs(yf).max())
+axs[1].set_title("Downsampled")
+axs[1].axis("tight")
+axs[2].imshow(np.abs(yf1), cmap="jet", vmin=0, vmax=0.005 * np.abs(yf1).max())
+axs[2].set_title("Resampled (no filter)")
+axs[2].axis("tight")
+plt.tight_layout()
+
+###############################################################################
+# Similarly, if we take a synthetic image containing a rapidly oscillating
+# pattern, where aliasing is easy to spot, we can see the difference between
+# our downsampled image with the one obtained by simply picking one every
+# four samples in each direction (which is equivalent to using ``sigma=0``).
+nz1, nx1 = 201, 201
+iz, ix = np.meshgrid(np.arange(nz1), np.arange(nx1), indexing="ij")
+xosc = np.sin(0.1 * np.sqrt((iz - nz1 // 2) ** 2 + (ix - nx1 // 2) ** 2) ** 2 / 10.0)
+
+Dop = pylops.signalprocessing.Downsample2D((nz1, nx1), factors=4)
+Dop_noaa = pylops.signalprocessing.Downsample2D((nz1, nx1), factors=4, sigma=0.0)
+
+fig, axs = plt.subplots(1, 3, figsize=(12, 4))
+axs[0].imshow(xosc, cmap="gray")
+axs[0].set_title("Original")
+axs[0].axis("tight")
+axs[1].imshow(Dop_noaa @ xosc, cmap="gray")
+axs[1].set_title("Subsampled (aliased)")
+axs[1].axis("tight")
+axs[2].imshow(Dop @ xosc, cmap="gray")
+axs[2].set_title("Downsampled (anti-aliased)")
+axs[2].axis("tight")
+plt.tight_layout()
+
+###############################################################################
+# Finally, we consider the inverse problem: given the downsampled data, can we
+# retrieve the original, finely sampled image? As the operator has many more
+# columns than rows, this problem is heavily underdetermined and we must
+# regularize it. Here we simply ask for a smooth solution by penalizing the
+# Laplacian of the model. We compare the estimated model with the adjoint,
+# which spreads each coarse sample back over the fine grid.
+x = datasets.face()[400:528:2, 400:528:2, 0].astype(np.float64)
+nz, nx = x.shape
+
+Dop = pylops.signalprocessing.Downsample2D((nz, nx), factors=2)
+y = Dop @ x
+
+xadj = Dop.H @ y
+D2op = pylops.Laplacian((nz, nx), weights=(1, 1), dtype="float64")
+xinv = pylops.optimization.leastsquares.regularized_inversion(
+    Dop, y.ravel(), [D2op], epsRs=[np.sqrt(0.1)], **dict(iter_lim=200)
+)[0]
+xinv = xinv.reshape(nz, nx)
+
+fig, axs = plt.subplots(1, 4, figsize=(14, 4))
+axs[0].imshow(x, cmap="gray", vmin=0, vmax=255)
+axs[0].set_title("Original")
+axs[0].axis("tight")
+axs[1].imshow(y, cmap="gray", vmin=0, vmax=255)
+axs[1].set_title("Downsampled")
+axs[1].axis("tight")
+axs[2].imshow(xadj, cmap="gray")
+axs[2].set_title("Adjoint")
+axs[2].axis("tight")
+axs[3].imshow(xinv, cmap="gray", vmin=0, vmax=255)
+axs[3].set_title("Inverse")
+axs[3].axis("tight")
+plt.tight_layout()
diff --git a/pylops/signalprocessing/__init__.py b/pylops/signalprocessing/__init__.py
index c9d58262..de21b40a 100755
--- a/pylops/signalprocessing/__init__.py
+++ b/pylops/signalprocessing/__init__.py
@@ -18,6 +18,7 @@
     Interp                          Interpolation operator.
     InterpCubicSpline               Cubic Spline Interpolation operator.
     Bilinear                        Bilinear interpolation operator.
+    Downsample2D                    2D downsampling operator.
     FFT                             One dimensional Fast-Fourier Transform.
     FFT2D                           Two dimensional Fast-Fourier Transform.
     FFTND                           N-dimensional Fast-Fourier Transform.
@@ -59,6 +60,7 @@
 from .interp import *
 from .interpspline import *
 from .bilinear import *
+from .downsample2d import *
 from .radon2d import *
 from .radon3d import *
 from .fourierradon2d import *
@@ -96,6 +98,7 @@
     "Interp",
     "InterpCubicSpline",
     "Bilinear",
+    "Downsample2D",
     "Radon2D",
     "Radon3D",
     "FourierRadon2D",
diff --git a/pylops/signalprocessing/downsample2d.py b/pylops/signalprocessing/downsample2d.py
new file mode 100644
index 00000000..f811eb1d
--- /dev/null
+++ b/pylops/signalprocessing/downsample2d.py
@@ -0,0 +1,235 @@
+__all__ = ["Downsample2D"]
+
+from typing import Literal
+
+import numpy as np
+
+from pylops import LinearOperator
+from pylops.signalprocessing import Convolve2D
+from pylops.utils._internal import _value_or_sized_to_tuple
+from pylops.utils.backend import get_array_module, get_normalize_axis_index
+from pylops.utils.decorators import reshaped
+from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray, SamplingLike
+
+
+def _gaussian_kernel1d(sigma: float, truncate: float) -> NDArray:
+    """Create a normalized, symmetric 1d Gaussian kernel.
+
+    The kernel is truncated at ``truncate`` standard deviations, leading to a
+    kernel of size :math:`2 \\lfloor \\text{truncate} \\sigma + 0.5 \\rfloor + 1`.
+    A unitary kernel (i.e., ``[1.]``) is returned when ``sigma=0``.
+    """
+    if sigma == 0.0:
+        return np.ones(1)
+    radius = int(truncate * sigma + 0.5)
+    x = np.arange(-radius, radius + 1)
+    h = np.exp(-0.5 * (x / sigma) ** 2)
+    return h / h.sum()
+
+
+class Downsample2D(LinearOperator):
+    r"""2D downsampling operator.
+
+    Downsample a two (or more) dimensional array along a pair of ``axes`` by
+    applying an anti-aliasing Gaussian filter followed by subsampling with
+    a given decimation factor in each of the two directions.
+
+    Parameters
+    ----------
+    dims : :obj:`list` or :obj:`int`
+        Number of samples for each dimension.
+    factors : :obj:`int` or :obj:`tuple`, optional
+        Decimation factors along each of the two ``axes``. If a single value is
+        provided, the same factor is used in both directions.
+    sigma : :obj:`float` or :obj:`tuple`, optional
+        Standard deviations (in number of samples) of the Gaussian filter along
+        each of the two ``axes``. If a single value is provided, the same
+        standard deviation is used in both directions. If ``None``, the
+        standard deviations are set to ``(factor - 1) / 2`` for each direction.
+    truncate : :obj:`float`, optional
+        Number of standard deviations at which the Gaussian filter is
+        truncated. The filter has ``2 * int(truncate * sigma + 0.5) + 1``
+        samples along each direction.
+    axes : :obj:`tuple`, optional
+        Axes along which downsampling is applied.
+    method : :obj:`str`, optional
+        Method used to calculate the Gaussian filtering (``auto``, ``direct``
+        or ``fft``) - see :func:`scipy.signal.convolve` for details.
+    dtype : :obj:`str`, optional
+        Type of elements in input array.
+    name : :obj:`str`, optional
+        Name of operator (to be used by :func:`pylops.utils.describe.describe`)
+
+    Attributes
+    ----------
+    h : :obj:`numpy.ndarray`
+        2d Gaussian filter applied prior to subsampling.
+    Cop : :obj:`pylops.signalprocessing.Convolve2D`
+        Gaussian filtering operator.
+    dims : :obj:`tuple`
+        Shape of the array after the adjoint, but before flattening.
+
+        For example, ``x_reshaped = (Op.H * y.ravel()).reshape(Op.dims)``.
+    dimsd : :obj:`tuple`
+        Shape of the array after the forward, but before flattening.
+
+        For example, ``y_reshaped = (Op * x.ravel()).reshape(Op.dimsd)``.
+    shape : :obj:`tuple`
+        Operator shape.
+    explicit : :obj:`bool`
+        Operator contains a matrix that can be solved explicitly (``True``) or
+        not (``False``).
+
+    Raises
+    ------
+    ValueError
+        If ``dims`` has less than 2 dimensions, if ``axes``, ``factors``, or
+        ``sigma`` do not contain 2 elements, if any element of ``factors`` is
+        smaller than 1 or larger than half the size of the corresponding axis,
+        or if any element of ``sigma`` is negative.
+
+    See Also
+    --------
+    pylops.signalprocessing.Convolve2D : 2D convolution operator
+    pylops.Restriction : Restriction (or sampling) operator
+
+    Notes
+    -----
+    The Downsample2D operator reduces the size of a two-dimensional array
+    :math:`\mathbf{x}` of size :math:`n_0 \times n_1` by a factor
+    :math:`f_0` and :math:`f_1` along the first and second direction,
+    respectively. Direct subsampling of the input array would however lead to
+    aliasing of any energy above the Nyquist wavenumber of the coarse grid;
+    for this reason the array is first smoothed by a separable Gaussian kernel
+
+    .. math::
+        h[p, q] = g_{\sigma_0}[p]\, g_{\sigma_1}[q], \qquad
+        g_\sigma[p] = \frac{e^{-p^2 / (2\sigma^2)}}
+        {\sum_{p'} e^{-p'^2 / (2\sigma^2)}}
+
+    with :math:`|p| \leq r_0`, :math:`|q| \leq r_1`, and
+    :math:`r_i = \lfloor \tau \sigma_i + 0.5 \rfloor` where :math:`\tau` is the
+    ``truncate`` parameter. In forward mode, filtering and subsampling are
+    applied one after the other
+
+    .. math::
+        y[i, j] = \sum_{p=-r_0}^{r_0} \sum_{q=-r_1}^{r_1}
+        h[p, q] \, x[f_0 i - p, f_1 j - q]
+        \quad \forall i=0,\ldots,\lceil n_0 / f_0 \rceil - 1,
+        \; j=0,\ldots,\lceil n_1 / f_1 \rceil - 1
+
+    where the input array is assumed to be zero-padded outside of its
+    boundaries. Since the adjoint of subsampling is zero-interleaving and the
+    adjoint of convolution is correlation, in adjoint mode the data is first
+    spread over the fine grid and then correlated with the same kernel
+
+    .. math::
+        x[k, l] = \sum_{p=-r_0}^{r_0} \sum_{q=-r_1}^{r_1}
+        h[p, q] \, \tilde{y}[k + p, l + q], \qquad
+        \tilde{y}[k, l] =
+        \begin{cases}
+        y[k / f_0, l / f_1] & k \bmod f_0 = 0 \land l \bmod f_1 = 0\\
+        0 & \text{otherwise}
+        \end{cases}
+
+    Note that, as the Gaussian kernel is real and symmetric, the operator
+    is effectively the composition of a self-adjoint smoothing operator and a
+    restriction operator.
+
+    """
+
+    def __init__(
+        self,
+        dims: InputDimsLike,
+        factors: int | InputDimsLike = 2,
+        sigma: float | SamplingLike | None = None,
+        truncate: float = 4.0,
+        axes: InputDimsLike = (-2, -1),
+        method: Literal["auto", "direct", "fft"] | None = "fft",
+        dtype: DTypeLike = "float64",
+        name: str = "D",
+    ) -> None:
+        # check dims
+        dims = _value_or_sized_to_tuple(dims)
+        if len(dims) < 2:
+            msg = "dims must contain at least 2 dimensions"
+            raise ValueError(msg)
+
+        # check axes
+        if len(axes) != 2:
+            msg = "axes must contain 2 elements"
+            raise ValueError(msg)
+        axes = tuple(get_normalize_axis_index()(ax, len(dims)) for ax in axes)
+
+        # check factors
+        factors = _value_or_sized_to_tuple(factors, repeat=2)
+        if len(factors) != 2:
+            msg = "factors must contain 2 elements"
+            raise ValueError(msg)
+
+        for f, ax in zip(factors, axes, strict=True):
+            if f < 1:
+                msg = "factors must be greater or equal to 1"
+                raise ValueError(msg)
+            if f > dims[ax] // 2:
+                msg = (
+                    f"factor={f} is larger than the half of the "
+                    f"number of samples ({dims[ax]}) along axis={ax}"
+                )
+                raise ValueError(msg)
+
+        if sigma is None:
+            sigma = tuple((f - 1) / 2.0 for f in factors)
+        else:
+            sigma = _value_or_sized_to_tuple(sigma, repeat=2)
+        if len(sigma) != 2:
+            msg = "sigma must contain 2 elements"
+            raise ValueError(msg)
+        if any(s < 0 for s in sigma):
+            msg = "sigma must be positive"
+            raise ValueError(msg)
+
+        self.axes = axes
+        self.factors = factors
+        self.sigma = sigma
+        self.truncate = truncate
+
+        # data dimensions after subsampling
+        dimsd = list(dims)
+        for f, ax in zip(factors, axes, strict=True):
+            dimsd[ax] = int(np.ceil(dims[ax] / f))
+
+        super().__init__(
+            dtype=np.dtype(dtype), dims=dims, dimsd=tuple(dimsd), name=name
+        )
+
+        # separable gaussian filter and associated convolution operator
+        h0 = _gaussian_kernel1d(sigma[0], truncate)
+        h1 = _gaussian_kernel1d(sigma[1], truncate)
+        self.h = np.outer(h0, h1).astype(self.dtype)
+        self.Cop = Convolve2D(
+            dims,
+            h=self.h,
+            offset=(h0.size // 2, h1.size // 2),
+            axes=axes,
+            method=method,
+            dtype=dtype,
+        )
+
+        # slices used to subsample the filtered model
+        self.slices = tuple(
+            slice(None, None, factors[axes.index(ax)]) if ax in axes else slice(None)
+            for ax in range(len(dims))
+        )
+
+    @reshaped
+    def _matvec(self, x: NDArray) -> NDArray:
+        y = self.Cop._matvec(x.ravel()).reshape(self.dims)
+        return y[self.slices]
+
+    @reshaped
+    def _rmatvec(self, x: NDArray) -> NDArray:
+        ncp = get_array_module(x)
+        y = ncp.zeros(self.dims, dtype=self.dtype)
+        y[self.slices] = x
+        return self.Cop._rmatvec(y.ravel()).reshape(self.dims)
diff --git a/pylops/signalprocessing/radon2d.py b/pylops/signalprocessing/radon2d.py
index 77bd0313..97c75f7b 100644
--- a/pylops/signalprocessing/radon2d.py
+++ b/pylops/signalprocessing/radon2d.py
@@ -144,7 +144,7 @@ def Radon2D(
     taxis: NDArray,
     haxis: NDArray,
     pxaxis: NDArray,
-    kind: Literal["linear", "parabolic", "hyperbolic"] = "linear",
+    kind: Literal["linear", "parabolic", "hyperbolic"] | Callable = "linear",
     centeredh: bool = True,
     interp: bool = True,
     onthefly: bool = False,
diff --git a/pylops/signalprocessing/radon3d.py b/pylops/signalprocessing/radon3d.py
index 96eabd0f..cbc7a935 100644
--- a/pylops/signalprocessing/radon3d.py
+++ b/pylops/signalprocessing/radon3d.py
@@ -164,7 +164,7 @@ def Radon3D(
     hxaxis: NDArray,
     pyaxis: NDArray,
     pxaxis: NDArray,
-    kind: Literal["linear", "parabolic", "hyperbolic"] = "linear",
+    kind: Literal["linear", "parabolic", "hyperbolic"] | Callable = "linear",
     centeredh: bool = True,
     interp: bool = True,
     onthefly: bool = False,
diff --git a/pytests/test_downsample.py b/pytests/test_downsample.py
new file mode 100644
index 00000000..2bb7e56a
--- /dev/null
+++ b/pytests/test_downsample.py
@@ -0,0 +1,156 @@
+import os
+
+if int(os.environ.get("TEST_CUPY_PYLOPS", 0)):
+    import cupy as np
+    from cupy.testing import assert_array_almost_equal
+
+    backend = "cupy"
+else:
+    import numpy as np
+    from numpy.testing import assert_array_almost_equal
+
+    backend = "numpy"
+
+import numpy as npp
+import pytest
+from scipy.ndimage import gaussian_filter
+
+from pylops.optimization.basic import lsqr
+from pylops.signalprocessing import Downsample2D
+from pylops.utils import dottest
+
+par1 = {
+    "ny": 21,
+    "nx": 15,
+    "factors": 3,
+    "imag": 0,
+    "dtype": "float64",
+}  # same factor, real
+par2 = {
+    "ny": 20,
+    "nx": 16,
+    "factors": (2, 4),
+    "imag": 0,
+    "dtype": "float64",
+}  # different factors, real
+par3 = {
+    "ny": 11,
+    "nx": 13,
+    "factors": 1,
+    "imag": 0,
+    "dtype": "float64",
+}  # unitary factor, real
+par1j = {
+    "ny": 21,
+    "nx": 15,
+    "factors": 3,
+    "imag": 1j,
+    "dtype": "complex128",
+}  # same factor, complex
+par2j = {
+    "ny": 20,
+    "nx": 16,
+    "factors": (2, 4),
+    "imag": 1j,
+    "dtype": "complex128",
+}  # different factors, complex
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        {"dims": (10,)},
+        {"dims": (10, 10), "axes": (0,)},
+        {"dims": (10, 10), "factors": (2, 2, 2)},
+        {"dims": (10, 10), "factors": 0},
+        {"dims": (10, 10), "factors": 11},
+        {"dims": (10, 10), "sigma": (1.0, 1.0, 1.0)},
+        {"dims": (10, 10), "sigma": -1.0},
+    ],
+)
+def test_Downsample2D_raises(kwargs):
+    """Check input validation of Downsample2D"""
+    with pytest.raises(ValueError):
+        Downsample2D(**kwargs)
+
+
+def test_Downsample2D_sigma():
+    """Check that a null sigma leads to pure subsampling"""
+    x = np.random.normal(0.0, 1.0, (12, 9))
+    Dop = Downsample2D((12, 9), factors=(3, 3), sigma=0.0)
+    assert Dop.h.shape == (1, 1)
+    assert_array_almost_equal(Dop @ x, x[::3, ::3], decimal=10)
+
+
+def test_Downsample2D_ndim():
+    """Check that Downsample2D can be applied to a subset of axes of a
+    3-dimensional array
+    """
+    Dop = Downsample2D((7, 9, 5), factors=2, axes=(0, 1))
+    assert Dop.dimsd == (4, 5, 5)
+    assert dottest(Dop, *Dop.shape, rtol=1e-6, backend=backend)
+
+
+@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par1j), (par2j)])
+def test_Downsample2D(par):
+    """Dot-test and shapes for Downsample2D"""
+    Dop = Downsample2D(
+        (par["ny"], par["nx"]), factors=par["factors"], dtype=par["dtype"]
+    )
+    factors = (
+        (par["factors"], par["factors"])
+        if isinstance(par["factors"], int)
+        else par["factors"]
+    )
+    assert Dop.dimsd == (
+        int(npp.ceil(par["ny"] / factors[0])),
+        int(npp.ceil(par["nx"] / factors[1])),
+    )
+    assert dottest(
+        Dop,
+        *Dop.shape,
+        rtol=1e-6,
+        complexflag=0 if par["imag"] == 0 else 3,
+        backend=backend,
+    )
+
+
+@pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j)])
+def test_Downsample2D_scipy(par):
+    """Compare Downsample2D forward with scipy
+    gaussian filtering plus subsampling"""
+    factors = (
+        (par["factors"], par["factors"])
+        if isinstance(par["factors"], int)
+        else par["factors"]
+    )
+    sigma = tuple((f - 1) / 2.0 for f in factors)
+
+    shape = (par["ny"], par["nx"])
+    x = np.random.normal(0.0, 1.0, shape) + par["imag"] * np.random.normal(
+        0.0, 1.0, shape
+    )
+    Dop = Downsample2D(
+        (par["ny"], par["nx"]), factors=par["factors"], dtype=par["dtype"]
+    )
+    y = Dop @ x
+
+    xnp = np.asnumpy(x) if backend == "cupy" else x
+    ynp = gaussian_filter(xnp, sigma=sigma, truncate=4.0, mode="constant")[
+        :: factors[0], :: factors[1]
+    ]
+    assert_array_almost_equal(y, np.asarray(ynp), decimal=10)
+
+
+@pytest.mark.parametrize("par", [(par3)])
+def test_Downsample2D_inverse(par):
+    """Invert Downsample2D when no decimation is applied (factors=1) as in
+    this case the operator is a square, invertible smoothing operator
+    """
+    x = np.random.normal(0.0, 1.0, (par["ny"], par["nx"]))
+    Dop = Downsample2D(
+        (par["ny"], par["nx"]), factors=par["factors"], sigma=0.6, dtype=par["dtype"]
+    )
+    y = Dop @ x
+    xinv = lsqr(Dop, y.ravel(), x0=np.zeros(Dop.shape[1]), niter=500, show=0)[0]
+    assert_array_almost_equal(x.ravel(), xinv, decimal=3)
diff --git a/pytests/test_radon.py b/pytests/test_radon.py
index 787fa635..8f299c73 100644
--- a/pytests/test_radon.py
+++ b/pytests/test_radon.py
@@ -126,6 +126,18 @@ def test_unknown_engine():
         _ = Radon3D(None, None, None, None, None, engine="foo")
 
 
+@pytest.mark.skipif(
+    int(os.environ.get("TEST_CUPY_PYLOPS", 0)) == 1, reason="Not CuPy enabled"
+)
+def test_Radon2D_unknown_kind():
+    """Check error is raised if unknown (and non-callable) kind is passed"""
+    t = np.arange(11, dtype=np.float64) * 0.005
+    h = np.arange(21, dtype=np.float64)
+    px = np.linspace(0, 2e-2, 21, dtype=np.float64)
+    with pytest.raises(NotImplementedError, match="Wrong kind of basis function"):
+        _ = Radon2D(t, h, px, kind="foo")
+
+
 @pytest.mark.skipif(
     int(os.environ.get("TEST_CUPY_PYLOPS", 0)) == 1, reason="Not CuPy enabled"
 )
@@ -189,6 +201,33 @@ def test_Radon2D(par, dtype):
     assert_array_almost_equal(x.ravel(), xinv, decimal=1)
 
 
+@pytest.mark.skipif(
+    int(os.environ.get("TEST_CUPY_PYLOPS", 0)) == 1, reason="Not CuPy enabled"
+)
+def test_Radon2D_callable_kind():
+    """Dot-test for Radon2D operator when kind is a custom callable"""
+    dt, dh = 0.005, 1
+    t = np.arange(par1["nt"], dtype=np.float64) * dt
+    h = np.arange(par1["nhx"], dtype=np.float64) * dh
+    px = np.linspace(0, par1["pxmax"], par1["npx"], dtype=np.float64)
+
+    def _linear(x, t, px):
+        return t + px * x
+
+    Rop = Radon2D(
+        t,
+        h,
+        px,
+        centeredh=par1["centeredh"],
+        interp=par1["interp"],
+        kind=_linear,
+        onthefly=False,
+        engine="numpy",
+        dtype=np.float64,
+    )
+    assert dottest(Rop, par1["nhx"] * par1["nt"], par1["npx"] * par1["nt"])
+
+
 @pytest.mark.skipif(
     int(os.environ.get("TEST_CUPY_PYLOPS", 0)) == 1, reason="Not CuPy enabled"
 )