Repository: gammapy/gammapy
Pull request #5264 (2024-05-16 15:45:37)

Introduce COPY_IF_NEEDED to allow compatibility with numpy 2.0
<!-- These are hidden commments. You can leave them or delete them. -->

<!-- Don't be shy! Github issues or pull requests are welcome any time. -->
<!-- But if you first want to chat, use the Gammapy Slack. -->
<!-- See https://gammapy.org/contact.html -->

**Description**
<!-- What is the motivation for this pull request? -->
<!-- Briefly: what changes are done here? -->
<!-- If this is to address a Github issue, mention its number and use the specific keywords to create a link -->
<!-- See docs.github.com/en/issues/tracking-your-work-with-issues/linking-a-pull-request-to-an-issue -->
<!-- Example: This PR is to resolve #42 -->

This pull request makes gammapy compatible with numpy 2.0 while keeping backward compatibility. 
The main issue is connected with the `np.array` copy semantics . See https://github.com/numpy/numpy/pull/25168
`Quantity(copy=False)` is now strictly interpreted and will fail if a copy is unavoidable. Because the behavior has been fully ported in astropy. See https://github.com/astropy/astropy/pull/16142. 
This creates a hundred fails in the devdeps CI. 
This is solved in a backward compatible manner by following the same approach as https://github.com/astropy/astropy/issues/16167 introducing a global variable `COPY_IF_NEEDED` set to `False` if numpy<2.0 and `None` otherwise.

An additional issue appeared in `WcsNDMap.to_nd_region_map_histogram` probably because units aren't propagated the same way with `np.apply_along_axis`. The backward compatible solution is to tear the units out before applying it.

**Dear reviewer**
<!-- Let the reviewer and Gammapy team know what you want: -->
<!-- * Is this ready for review? Or is it work in progress and you want some feedback? -->
<!-- * Do you want to go through review here? Or if someone should just finish this up and merge it? -->
<!-- Do you have any specific questions, e.g. about API or implementation? -->
<!-- Do you include a test executing the new code you're adding (to make sure it runs)? -->
<!-- Do you include some documentation? Is it needed? -->



Patch:
diff --git a/gammapy/data/pointing.py b/gammapy/data/pointing.py
index 4470c5c12..eb46e5452 100644
--- a/gammapy/data/pointing.py
+++ b/gammapy/data/pointing.py
@@ -18,6 +18,7 @@ from astropy.io import fits
 from astropy.table import Table
 from astropy.units import Quantity
 from astropy.utils import lazyproperty
+from gammapy.utils.compat import COPY_IF_NEEDED
 from gammapy.utils.deprecation import GammapyDeprecationWarning
 from gammapy.utils.fits import earth_location_from_dict, earth_location_to_dict
 from gammapy.utils.scripts import make_path
@@ -442,8 +443,12 @@ class FixedPointingInfo:
             alt = self.fixed_altaz.alt
             az = self.fixed_altaz.az
             return SkyCoord(
-                alt=u.Quantity(np.full(obstime.shape, alt.deg), u.deg, copy=False),
-                az=u.Quantity(np.full(obstime.shape, az.deg), u.deg, copy=False),
+                alt=u.Quantity(
+                    np.full(obstime.shape, alt.deg), u.deg, copy=COPY_IF_NEEDED
+                ),
+                az=u.Quantity(
+                    np.full(obstime.shape, az.deg), u.deg, copy=COPY_IF_NEEDED
+                ),
                 frame=frame,
             )
 
diff --git a/gammapy/estimators/points/core.py b/gammapy/estimators/points/core.py
index fb3aa0e7f..87dd6856c 100644
--- a/gammapy/estimators/points/core.py
+++ b/gammapy/estimators/points/core.py
@@ -267,7 +267,7 @@ class FluxPoints(FluxMaps):
         table = self.to_table(sed_type=sed_type, format=format)
 
         if ".fits" not in filename.suffixes:
-            table.write(filename, overwrite=overwrite)
+            table.write(filename)
             return
 
         primary_hdu = fits.PrimaryHDU()
diff --git a/gammapy/irf/background.py b/gammapy/irf/background.py
index 4adf3997c..b2b62c7f9 100644
--- a/gammapy/irf/background.py
+++ b/gammapy/irf/background.py
@@ -8,6 +8,7 @@ import matplotlib.pyplot as plt
 from matplotlib.colors import LogNorm
 from gammapy.maps import MapAxes, MapAxis
 from gammapy.maps.axes import UNIT_STRING_FORMAT
+from gammapy.utils.compat import COPY_IF_NEEDED
 from gammapy.visualization.utils import add_colorbar
 from .core import IRF
 from .io import gadf_is_pointlike
@@ -73,7 +74,7 @@ class BackgroundIRF(IRF):
         data = table[bkg_name].quantity[0].T
 
         if data.unit == "" or isinstance(data.unit, u.UnrecognizedUnit):
-            data = u.Quantity(data.value, "s-1 MeV-1 sr-1", copy=False)
+            data = u.Quantity(data.value, "s-1 MeV-1 sr-1", copy=COPY_IF_NEEDED)
             log.warning(
                 "Invalid unit found in background table! Assuming (s-1 MeV-1 sr-1)"
             )
diff --git a/gammapy/irf/core.py b/gammapy/irf/core.py
index fd3124240..a0fdabb77 100644
--- a/gammapy/irf/core.py
+++ b/gammapy/irf/core.py
@@ -10,6 +10,7 @@ from astropy.io import fits
 from astropy.table import Table
 from astropy.utils import lazyproperty
 from gammapy.maps import Map, MapAxes, MapAxis, RegionGeom
+from gammapy.utils.compat import COPY_IF_NEEDED
 from gammapy.utils.integrate import trapz_loglog
 from gammapy.utils.interpolation import (
     ScaledRegularGridInterpolator,
@@ -205,7 +206,7 @@ class IRF(metaclass=abc.ABCMeta):
     @property
     def quantity(self):
         """Quantity as a `~astropy.units.Quantity` object."""
-        return u.Quantity(self.data, unit=self.unit, copy=False)
+        return u.Quantity(self.data, unit=self.unit, copy=COPY_IF_NEEDED)
 
     @quantity.setter
     def quantity(self, val):
@@ -216,7 +217,7 @@ class IRF(metaclass=abc.ABCMeta):
         value : `~astropy.units.Quantity`
            Quantity.
         """
-        val = u.Quantity(val, copy=False)
+        val = u.Quantity(val, copy=COPY_IF_NEEDED)
         self.data = val.value
         self._unit = val.unit
 
@@ -288,7 +289,7 @@ class IRF(metaclass=abc.ABCMeta):
         for key, value in kwargs.items():
             coord = kwargs.get(key, value)
             if coord is not None:
-                coords_default[key] = u.Quantity(coord, copy=False)
+                coords_default[key] = u.Quantity(coord, copy=COPY_IF_NEEDED)
 
         data = self._interpolate(coords_default.values(), method=method)
 
diff --git a/gammapy/maps/axes.py b/gammapy/maps/axes.py
index 91df33cd1..c12ef6d72 100644
--- a/gammapy/maps/axes.py
+++ b/gammapy/maps/axes.py
@@ -13,6 +13,7 @@ from astropy.table import Column, Table, hstack
 from astropy.time import Time
 from astropy.utils import lazyproperty
 import matplotlib.pyplot as plt
+from gammapy.utils.compat import COPY_IF_NEEDED
 from gammapy.utils.interpolation import interpolation_scale
 from gammapy.utils.time import time_ref_from_dict, time_ref_to_dict
 from .utils import INVALID_INDEX, INVALID_VALUE, edges_from_lo_hi
@@ -281,7 +282,7 @@ class MapAxis:
     def edges(self):
         """Return an array of bin edges."""
         pix = np.arange(self.nbin + 1, dtype=float) - 0.5
-        return u.Quantity(self.pix_to_coord(pix), self._unit, copy=False)
+        return u.Quantity(self.pix_to_coord(pix), self._unit, copy=COPY_IF_NEEDED)
 
     @property
     def edges_min(self):
@@ -453,7 +454,7 @@ class MapAxis:
     def center(self):
         """Return an array of bin centers."""
         pix = np.arange(self.nbin, dtype=float)
-        return u.Quantity(self.pix_to_coord(pix), self._unit, copy=False)
+        return u.Quantity(self.pix_to_coord(pix), self._unit, copy=COPY_IF_NEEDED)
 
     @lazyproperty
     def bin_width(self):
@@ -802,7 +803,7 @@ class MapAxis:
         """
         pix = pix - self._pix_offset
         values = self._transform.pix_to_coord(pix=pix)
-        return u.Quantity(values, unit=self.unit, copy=False)
+        return u.Quantity(values, unit=self.unit, copy=COPY_IF_NEEDED)
 
     def wrap_coord(self, coord):
         """Wrap coords between axis edges for a periodic boundary condition
@@ -862,7 +863,7 @@ class MapAxis:
         """
         if self._boundary_type == BoundaryEnum.periodic:
             coord = self.wrap_coord(coord)
-        coord = u.Quantity(coord, self.unit, copy=False).value
+        coord = u.Quantity(coord, self.unit, copy=COPY_IF_NEEDED).value
         pix = self._transform.coord_to_pix(coord=coord)
         return np.array(pix + self._pix_offset, ndmin=1)
 
@@ -885,7 +886,7 @@ class MapAxis:
         """
         if self._boundary_type == BoundaryEnum.periodic:
             coord = self.wrap_coord(coord)
-        coord = u.Quantity(coord, self.unit, copy=False, ndmin=1).value
+        coord = u.Quantity(coord, self.unit, copy=COPY_IF_NEEDED, ndmin=1).value
         edges = self.edges.value
         idx = np.digitize(coord, edges) - 1
 
diff --git a/gammapy/maps/coord.py b/gammapy/maps/coord.py
index 8420f9146..3129b4709 100644
--- a/gammapy/maps/coord.py
+++ b/gammapy/maps/coord.py
@@ -4,6 +4,7 @@ import html
 import numpy as np
 from astropy import units as u
 from astropy.coordinates import SkyCoord
+from gammapy.utils.compat import COPY_IF_NEEDED
 
 __all__ = ["MapCoord"]
 
@@ -105,13 +106,13 @@ class MapCoord:
     @property
     def theta(self):
         """Theta co-latitude angle in radians."""
-        theta = u.Quantity(self.lat, unit="deg", copy=False).to_value("rad")
+        theta = u.Quantity(self.lat, unit="deg", copy=COPY_IF_NEEDED).to_value("rad")
         return np.pi / 2.0 - theta
 
     @property
     def phi(self):
         """Phi longitude angle in radians."""
-        phi = u.Quantity(self.lon, unit="deg", copy=False).to_value("rad")
+        phi = u.Quantity(self.lon, unit="deg", copy=COPY_IF_NEEDED).to_value("rad")
         return phi
 
     @property
@@ -257,10 +258,10 @@ class MapCoord:
             lon, lat, frame = skycoord_to_lonlat(self.skycoord, frame=frame)
             data = copy.deepcopy(self._data)
             if isinstance(self.lon, u.Quantity):
-                lon = u.Quantity(lon, unit="deg", copy=False)
+                lon = u.Quantity(lon, unit="deg", copy=COPY_IF_NEEDED)
 
             if isinstance(self.lon, u.Quantity):
-                lat = u.Quantity(lat, unit="deg", copy=False)
+                lat = u.Quantity(lat, unit="deg", copy=COPY_IF_NEEDED)
 
             data["lon"] = lon
             data["lat"] = lat
diff --git a/gammapy/maps/core.py b/gammapy/maps/core.py
index 415b50160..a52fcb011 100644
--- a/gammapy/maps/core.py
+++ b/gammapy/maps/core.py
@@ -11,6 +11,7 @@ from astropy import units as u
 from astropy.io import fits
 import matplotlib.pyplot as plt
 import gammapy.utils.parallel as parallel
+from gammapy.utils.compat import COPY_IF_NEEDED
 from gammapy.utils.random import InverseCDFSampler, get_random_state
 from gammapy.utils.scripts import make_path
 from gammapy.utils.types import JsonQuantityDecoder
@@ -128,7 +129,7 @@ class Map(abc.ABC):
     @property
     def quantity(self):
         """Map data as a `~astropy.units.Quantity` object."""
-        return u.Quantity(self.data, self.unit, copy=False)
+        return u.Quantity(self.data, self.unit, copy=COPY_IF_NEEDED)
 
     @quantity.setter
     def quantity(self, val):
@@ -139,7 +140,7 @@ class Map(abc.ABC):
         val : `~astropy.units.Quantity`
            Quantity.
         """
-        val = u.Quantity(val, copy=False)
+        val = u.Quantity(val, copy=COPY_IF_NEEDED)
 
         self.data = val.value
         self._unit = val.unit
@@ -1662,7 +1663,7 @@ class Map(abc.ABC):
         cumsum = self.cumsum(axis_name=axis_name)
         cumsum = cumsum.pad(pad_width=1, axis_name=axis_name, mode="edge")
         return u.Quantity(
-            cumsum.interp_by_coord(coords, **kwargs), cumsum.unit, copy=False
+            cumsum.interp_by_coord(coords, **kwargs), cumsum.unit, copy=COPY_IF_NEEDED
         )
 
     def normalize(self, axis_name=None):
@@ -1872,7 +1873,7 @@ class Map(abc.ABC):
             else:
                 raise ValueError("Map Arithmetic: Inconsistent geometries.")
         else:
-            q = u.Quantity(other, copy=False)
+            q = u.Quantity(other, copy=COPY_IF_NEEDED)
 
         out = self.copy() if copy else self
         out.quantity = operator(out.quantity, q)
@@ -1899,25 +1900,25 @@ class Map(abc.ABC):
         return self._arithmetics(np.add, other, copy=True)
 
     def __iadd__(self, other):
-        return self._arithmetics(np.add, other, copy=False)
+        return self._arithmetics(np.add, other, copy=COPY_IF_NEEDED)
 
     def __sub__(self, other):
         return self._arithmetics(np.subtract, other, copy=True)
 
     def __isub__(self, other):
-        return self._arithmetics(np.subtract, other, copy=False)
+        return self._arithmetics(np.subtract, other, copy=COPY_IF_NEEDED)
 
     def __mul__(self, other):
         return self._arithmetics(np.multiply, other, copy=True)
 
     def __imul__(self, other):
-        return self._arithmetics(np.multiply, other, copy=False)
+        return self._arithmetics(np.multiply, other, copy=COPY_IF_NEEDED)
 
     def __truediv__(self, other):
         return self._arithmetics(np.true_divide, other, copy=True)
 
     def __itruediv__(self, other):
-        return self._arithmetics(np.true_divide, other, copy=False)
+        return self._arithmetics(np.true_divide, other, copy=COPY_IF_NEEDED)
 
     def __le__(self, other):
         return self._arithmetics(np.less_equal, other, copy=True)
@@ -1950,13 +1951,13 @@ class Map(abc.ABC):
         return self._boolean_arithmetics(np.logical_xor, other, copy=True)
 
     def __iand__(self, other):
-        return self._boolean_arithmetics(np.logical_and, other, copy=False)
+        return self._boolean_arithmetics(np.logical_and, other, copy=COPY_IF_NEEDED)
 
     def __ior__(self, other):
-        return self._boolean_arithmetics(np.logical_or, other, copy=False)
+        return self._boolean_arithmetics(np.logical_or, other, copy=COPY_IF_NEEDED)
 
     def __ixor__(self, other):
-        return self._boolean_arithmetics(np.logical_xor, other, copy=False)
+        return self._boolean_arithmetics(np.logical_xor, other, copy=COPY_IF_NEEDED)
 
     def __array__(self):
         return self.data
diff --git a/gammapy/maps/hpx/ndmap.py b/gammapy/maps/hpx/ndmap.py
index 341bd706a..0d885d74f 100644
--- a/gammapy/maps/hpx/ndmap.py
+++ b/gammapy/maps/hpx/ndmap.py
@@ -6,6 +6,7 @@ from astropy.coordinates import SkyCoord
 from astropy.io import fits
 from regions import PointSkyRegion
 import matplotlib.pyplot as plt
+from gammapy.utils.compat import COPY_IF_NEEDED
 from gammapy.utils.units import unit_from_fits_image_hdu
 from ..coord import MapCoord
 from ..geom import pix_tuple_to_idx
@@ -755,7 +756,9 @@ class HpxNDMap(HpxMap):
                 idx = np.clip(idx, 0, len(ax.center) - 2)
 
                 w = ax.center[idx + 1] - ax.center[idx]
-                c = u.Quantity(coords[ax.name], ax.center.unit, copy=False).value
+                c = u.Quantity(
+                    coords[ax.name], ax.center.unit, copy=COPY_IF_NEEDED
+                ).value
 
                 if i & (1 << j):
                     wt *= (c - ax.center[idx].value) / w.value
diff --git a/gammapy/maps/wcs/geom.py b/gammapy/maps/wcs/geom.py
index 6995b6030..f5b21796f 100644
--- a/gammapy/maps/wcs/geom.py
+++ b/gammapy/maps/wcs/geom.py
@@ -17,6 +17,7 @@ from astropy.wcs.utils import (
 )
 from regions import RectangleSkyRegion
 from gammapy.utils.array import round_up_to_even, round_up_to_odd
+from gammapy.utils.compat import COPY_IF_NEEDED
 from ..axes import MapAxes
 from ..coord import MapCoord, skycoord_to_lonlat
 from ..geom import Geom, get_shape, pix_tuple_to_idx
@@ -670,8 +671,8 @@ class WcsGeom(Geom):
             coords = self._wcs.wcs_pix2world(pix[0], pix[1], 0)
 
         coords = (
-            u.Quantity(coords[0], unit="deg", copy=False),
-            u.Quantity(coords[1], unit="deg", copy=False),
+            u.Quantity(coords[0], unit="deg", copy=COPY_IF_NEEDED),
+            u.Quantity(coords[1], unit="deg", copy=COPY_IF_NEEDED),
         )
 
         coords += self.axes.pix_to_coord(pix[self._slice_non_spatial_axes])
@@ -847,7 +848,9 @@ class WcsGeom(Geom):
         # TODO: for non-negative cdelt a negative solid angle is returned
         #  find out why and fix properly
 
-        value = np.abs(u.Quantity(area_low_right + area_up_left, "sr", copy=False))
+        value = np.abs(
+            u.Quantity(area_low_right + area_up_left, "sr", copy=COPY_IF_NEEDED)
+        )
         if self.is_regular:
             value = value.reshape(self.data_shape_image)
         return value
diff --git a/gammapy/maps/wcs/ndmap.py b/gammapy/maps/wcs/ndmap.py
index 34e4a3e1e..cf57345bb 100644
--- a/gammapy/maps/wcs/ndmap.py
+++ b/gammapy/maps/wcs/ndmap.py
@@ -726,16 +726,12 @@ class WcsNDMap(WcsMap):
 
         # This is likely not the most efficient way to do this
         data = np.apply_along_axis(
-            lambda a: np.histogram(a, bins=bins_axis.edges, density=density)[0],
+            lambda a: np.histogram(a, bins=bins_axis.edges.value, density=density)[0],
             axis=-1,
-            arr=quantity,
+            arr=quantity.to_value(bins_axis.unit),
         )
 
-        if density:
-            unit = 1.0 / bins_axis.unit
-            data = data.value
-        else:
-            unit = ""
+        unit = 1.0 / bins_axis.unit if density else ""
 
         return RegionNDMap.from_geom(geom=geom_hist, data=data, unit=unit)
 
diff --git a/gammapy/modeling/models/cube.py b/gammapy/modeling/models/cube.py
index 9ffeeb340..3f135ec1c 100644
--- a/gammapy/modeling/models/cube.py
+++ b/gammapy/modeling/models/cube.py
@@ -10,6 +10,7 @@ from astropy.time import Time
 from gammapy.maps import Map, MapAxis, WcsGeom
 from gammapy.modeling import Covariance, Parameters
 from gammapy.modeling.parameter import _get_parameters_str
+from gammapy.utils.compat import COPY_IF_NEEDED
 from gammapy.utils.fits import LazyFitsData
 from gammapy.utils.scripts import make_name, make_path
 from .core import Model, ModelBase, Models
@@ -1189,8 +1190,8 @@ def create_fermi_isotropic_diffuse_model(filename, **kwargs):
         Fermi isotropic diffuse sky model.
     """
     vals = np.loadtxt(make_path(filename))
-    energy = u.Quantity(vals[:, 0], "MeV", copy=False)
-    values = u.Quantity(vals[:, 1], "MeV-1 s-1 cm-2", copy=False)
+    energy = u.Quantity(vals[:, 0], "MeV", copy=COPY_IF_NEEDED)
+    values = u.Quantity(vals[:, 1], "MeV-1 s-1 cm-2", copy=COPY_IF_NEEDED)
 
     kwargs.setdefault("interp_kwargs", {"fill_value": None})
 
diff --git a/gammapy/modeling/models/spatial.py b/gammapy/modeling/models/spatial.py
index 37a7e1f82..fe0f4dd6d 100644
--- a/gammapy/modeling/models/spatial.py
+++ b/gammapy/modeling/models/spatial.py
@@ -21,6 +21,7 @@ from regions import (
 import matplotlib.pyplot as plt
 from gammapy.maps import HpxNDMap, Map, MapCoord, WcsGeom, WcsNDMap
 from gammapy.modeling import Parameter, Parameters
+from gammapy.utils.compat import COPY_IF_NEEDED
 from gammapy.utils.gauss import Gauss2DPDF
 from gammapy.utils.interpolation import interpolation_scale
 from gammapy.utils.regions import region_circle_to_ellipse, region_to_frame
@@ -651,7 +652,7 @@ class GaussianSpatialModel(SpatialModel):
             norm = (1 / (2 * np.pi * sigma * minor_axis)).to_value("sr-1")
 
         exponent = -0.5 * ((1 - np.cos(sep)) / a)
-        return u.Quantity(norm * np.exp(exponent).value, "sr-1", copy=False)
+        return u.Quantity(norm * np.exp(exponent).value, "sr-1", copy=COPY_IF_NEEDED)
 
     def to_region(self, x_sigma=1.5, **kwargs):
         r"""Model outline at a given number of :math:`\sigma`.
@@ -937,7 +938,7 @@ class DiskSpatialModel(SpatialModel):
         in_ellipse = DiskSpatialModel._evaluate_smooth_edge(
             sep - sigma_eff, sigma_eff * edge_width
         )
-        return u.Quantity(norm * in_ellipse, "sr-1", copy=False)
+        return u.Quantity(norm * in_ellipse, "sr-1", copy=COPY_IF_NEEDED)
 
     def to_region(self, **kwargs):
         """Model outline as a `~regions.EllipseSkyRegion`."""
@@ -1465,7 +1466,7 @@ class TemplateSpatialModel(SpatialModel):
 
         val = self.map.interp_by_coord(coord, **self._interp_kwargs)
         val = np.clip(val, 0, a_max=None)
-        return u.Quantity(val, self.map.unit, copy=False)
+        return u.Quantity(val, self.map.unit, copy=COPY_IF_NEEDED)
 
     @property
     def position_lonlat(self):
@@ -1617,7 +1618,7 @@ class TemplateNDSpatialModel(SpatialModel):
         val = self.map.interp_by_coord(coord, **self._interp_kwargs)
         val = np.clip(val, 0, a_max=None)
 
-        return u.Quantity(val, self.map.unit, copy=False)
+        return u.Quantity(val, self.map.unit, copy=COPY_IF_NEEDED)
 
     def write(self, overwrite=False):
         """
diff --git a/gammapy/modeling/models/spectral.py b/gammapy/modeling/models/spectral.py
index 0db2a3523..b238a8648 100644
--- a/gammapy/modeling/models/spectral.py
+++ b/gammapy/modeling/models/spectral.py
@@ -18,6 +18,7 @@ import matplotlib.pyplot as plt
 from gammapy.maps import MapAxis, RegionNDMap
 from gammapy.maps.axes import UNIT_STRING_FORMAT
 from gammapy.modeling import Parameter, Parameters
+from gammapy.utils.compat import COPY_IF_NEEDED
 from gammapy.utils.deprecation import GammapyDeprecationWarning
 from gammapy.utils.integrate import trapz_loglog
 from gammapy.utils.interpolation import (
@@ -646,7 +647,7 @@ class SpectralModel(ModelBase):
 
         def f(x):
             # scale by 1e12 to achieve better precision
-            energy = u.Quantity(x, eunit, copy=False)
+            energy = u.Quantity(x, eunit, copy=COPY_IF_NEEDED)
             y = self(energy).to_value(value.unit)
             return 1e12 * (y - value.value)
 
@@ -1761,7 +1762,7 @@ class TemplateSpectralModel(SpectralModel):
         meta=None,
     ):
         self.energy = energy
-        self.values = u.Quantity(values, copy=False)
+        self.values = u.Quantity(values, copy=COPY_IF_NEEDED)
         self.meta = {} if meta is None else meta
         interp_kwargs = interp_kwargs or {}
         interp_kwargs.setdefault("values_scale", "log")
@@ -1821,7 +1822,9 @@ class TemplateSpectralModel(SpectralModel):
         # Get spectrum values (no interpolation, take closest value for param)
         table_spectra = Table.read(filename, hdu="SPECTRA")
         idx = np.abs(table_spectra["PARAMVAL"] - param).argmin()
-        values = u.Quantity(table_spectra[idx][1], "", copy=False)  # no dimension
+        values = u.Quantity(
+            table_spectra[idx][1], "", copy=COPY_IF_NEEDED
+        )  # no dimension
 
         kwargs.setdefault("interp_kwargs", {"values_scale": "lin"})
         return cls(energy=energy, values=values, **kwargs)
@@ -1923,7 +1926,7 @@ class TemplateNDSpectralModel(SpectralModel):
         ]
 
         val = self.map.interp_by_pix(pixels, **self._interp_kwargs)
-        return u.Quantity(val, self.map.unit, copy=False)
+        return u.Quantity(val, self.map.unit, copy=COPY_IF_NEEDED)
 
     def write(self, overwrite=False):
         """
@@ -2023,7 +2026,7 @@ class EBLAbsorptionNormSpectralModel(SpectralModel):
         # set values log centers
         self.param = param
         self.energy = energy
-        self.data = u.Quantity(data, copy=False)
+        self.data = u.Quantity(data, copy=COPY_IF_NEEDED)
 
         interp_kwargs = interp_kwargs or {}
         interp_kwargs.setdefault("points_scale", ("lin", "log"))
@@ -2117,10 +2120,10 @@ class EBLAbsorptionNormSpectralModel(SpectralModel):
         # Get energy values
         table_energy = Table.read(filename, hdu="ENERGIES")
         energy_lo = u.Quantity(
-            table_energy["ENERG_LO"], "keV", copy=False
+            table_energy["ENERG_LO"], "keV", copy=COPY_IF_NEEDED
         )  # unit not stored in file
         energy_hi = u.Quantity(
-            table_energy["ENERG_HI"], "keV", copy=False
+            table_energy["ENERG_HI"], "keV", copy=COPY_IF_NEEDED
         )  # unit not stored in file
         energy = np.sqrt(energy_lo * energy_hi)
 
diff --git a/gammapy/modeling/models/spectral_crab.py b/gammapy/modeling/models/spectral_crab.py
index 110eac2c6..de2520451 100644
--- a/gammapy/modeling/models/spectral_crab.py
+++ b/gammapy/modeling/models/spectral_crab.py
@@ -2,6 +2,7 @@
 import numpy as np
 from astropy import units as u
 from gammapy.modeling import Parameter
+from gammapy.utils.compat import COPY_IF_NEEDED
 from .spectral import (
     ExpCutoffPowerLawSpectralModel,
     LogParabolaSpectralModel,
@@ -31,7 +32,7 @@ class MeyerCrabSpectralModel(SpectralModel):
         polynomial = np.poly1d(MeyerCrabSpectralModel.coefficients)
         log_energy = np.log10(energy.to_value("TeV"))
         log_flux = polynomial(log_energy)
-        flux = u.Quantity(np.power(10, log_flux), "erg / (cm2 s)", copy=False)
+        flux = u.Quantity(np.power(10, log_flux), "erg / (cm2 s)", copy=COPY_IF_NEEDED)
         return norm * flux / energy**2
 
 
diff --git a/gammapy/modeling/models/temporal.py b/gammapy/modeling/models/temporal.py
index d2a6fa192..45e3b61c3 100644
--- a/gammapy/modeling/models/temporal.py
+++ b/gammapy/modeling/models/temporal.py
@@ -10,6 +10,7 @@ from astropy.time import Time
 from astropy.utils import lazyproperty
 from gammapy.maps import MapAxis, RegionNDMap, TimeMapAxis
 from gammapy.modeling import Parameter
+from gammapy.utils.compat import COPY_IF_NEEDED
 from gammapy.utils.random import InverseCDFSampler, get_random_state
 from gammapy.utils.scripts import make_path
 from gammapy.utils.time import time_ref_from_dict, time_ref_to_dict
@@ -730,7 +731,7 @@ class LightCurveTemplateTemporalModel(TemporalModel):
             coords, method=self.method, values_scale=self.values_scale
         )
         val = np.clip(val, 0, a_max=None)
-        return u.Quantity(val, unit=self.map.unit, copy=False)
+        return u.Quantity(val, unit=self.map.unit, copy=COPY_IF_NEEDED)
 
     def integral(self, t_min, t_max, oversampling_factor=100, **kwargs):
         if self.is_energy_dependent:
diff --git a/gammapy/utils/compat.py b/gammapy/utils/compat.py
new file mode 100644
index 000000000..0363225ad
--- /dev/null
+++ b/gammapy/utils/compat.py
@@ -0,0 +1,12 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+import numpy as np
+from astropy.utils import minversion
+
+__all__ = [
+    "COPY_IF_NEEDED",
+]
+
+# This is copied from astropy.utils.compat
+NUMPY_LT_2_0 = not minversion(np, "2.0.0.dev")
+
+COPY_IF_NEEDED = False if NUMPY_LT_2_0 else None
diff --git a/gammapy/utils/interpolation.py b/gammapy/utils/interpolation.py
index 2f4c1c0ef..3d8968079 100644
--- a/gammapy/utils/interpolation.py
+++ b/gammapy/utils/interpolation.py
@@ -5,6 +5,7 @@ from itertools import compress
 import numpy as np
 import scipy.interpolate
 from astropy import units as u
+from .compat import COPY_IF_NEEDED
 
 __all__ = [
     "interpolate_profile",
@@ -162,7 +163,7 @@ class InterpolationScale:
 
     def __call__(self, values):
         if hasattr(self, "_unit"):
-            values = u.Quantity(values, copy=False).to_value(self._unit)
+            values = u.Quantity(values, copy=COPY_IF_NEEDED).to_value(self._unit)
         else:
             if isinstance(values, u.Quantity):
                 self._unit = values.unit
@@ -178,7 +179,7 @@ class InterpolationScale:
     def inverse(self, values):
         values = self._inverse(values)
         if hasattr(self, "_unit"):
-            return u.Quantity(values, self._unit, copy=False)
+            return u.Quantity(values, self._unit, copy=COPY_IF_NEEDED)
         else:
             return values
 


Test patch:
diff --git a/gammapy/estimators/points/tests/test_core.py b/gammapy/estimators/points/tests/test_core.py
index dc2e023f0..867e35e17 100644
--- a/gammapy/estimators/points/tests/test_core.py
+++ b/gammapy/estimators/points/tests/test_core.py
@@ -211,14 +211,12 @@ class TestFluxPoints:
         assert str(flux_points) == str(actual)
 
     def test_write_ecsv(self, tmp_path, flux_points):
-        filename = tmp_path / "flux_points.ecsv"
-        filename.touch()
         flux_points.write(
-            filename,
+            tmp_path / "flux_points.ecsv",
             sed_type=flux_points.sed_type_init,
             overwrite=True,
         )
-        actual = FluxPoints.read(filename)
+        actual = FluxPoints.read(tmp_path / "flux_points.ecsv")
         actual._data.pop("is_ul", None)
         flux_points._data.pop("is_ul", None)
         assert str(flux_points) == str(actual)
diff --git a/gammapy/irf/tests/test_rad_max.py b/gammapy/irf/tests/test_rad_max.py
index 87dd1f650..c9810a2b1 100644
--- a/gammapy/irf/tests/test_rad_max.py
+++ b/gammapy/irf/tests/test_rad_max.py
@@ -5,6 +5,7 @@ from numpy.testing import assert_allclose
 import astropy.units as u
 from gammapy.irf import EffectiveAreaTable2D, RadMax2D
 from gammapy.maps import MapAxis
+from gammapy.utils.compat import COPY_IF_NEEDED
 from gammapy.utils.testing import mpl_plot_check
 
 
@@ -52,7 +53,7 @@ def test_rad_max_from_irf():
     )
     offset_axis = MapAxis.from_bounds(0 * u.deg, 3 * u.deg, nbin=o_bins, name="offset")
     aeff = EffectiveAreaTable2D(
-        data=u.Quantity(np.ones((e_bins, o_bins)), u.m**2, copy=False),
+        data=u.Quantity(np.ones((e_bins, o_bins)), u.m**2, copy=COPY_IF_NEEDED),
         axes=[energy_axis, offset_axis],
         is_pointlike=True,
     )
diff --git a/gammapy/modeling/models/tests/test_spectral.py b/gammapy/modeling/models/tests/test_spectral.py
index fcea0dc79..d8f3e21a0 100644
--- a/gammapy/modeling/models/tests/test_spectral.py
+++ b/gammapy/modeling/models/tests/test_spectral.py
@@ -32,6 +32,7 @@ from gammapy.modeling.models import (
     TemplateNDSpectralModel,
     TemplateSpectralModel,
 )
+from gammapy.utils.compat import COPY_IF_NEEDED
 from gammapy.utils.scripts import make_path
 from gammapy.utils.testing import (
     assert_quantity_allclose,
@@ -1234,10 +1235,10 @@ def test_template_ND_EBL(tmpdir):
     # Get energy values
     table_energy = Table.read(filename, hdu="ENERGIES")
     energy_lo = u.Quantity(
-        table_energy["ENERG_LO"], "keV", copy=False
+        table_energy["ENERG_LO"], "keV", copy=COPY_IF_NEEDED
     )  # unit not stored in file
     energy_hi = u.Quantity(
-        table_energy["ENERG_HI"], "keV", copy=False
+        table_energy["ENERG_HI"], "keV", copy=COPY_IF_NEEDED
     )  # unit not stored in file
     energy = np.sqrt(energy_lo * energy_hi)
 
diff --git a/gammapy/utils/testing.py b/gammapy/utils/testing.py
index 9fb13e388..81bac784f 100644
--- a/gammapy/utils/testing.py
+++ b/gammapy/utils/testing.py
@@ -9,6 +9,7 @@ from astropy.coordinates import SkyCoord
 from astropy.time import Time
 from astropy.utils.introspection import minversion
 import matplotlib.pyplot as plt
+from .compat import COPY_IF_NEEDED
 
 __all__ = [
     "assert_quantity_allclose",
@@ -174,9 +175,9 @@ def assert_quantity_allclose(actual, desired, rtol=1.0e-7, atol=None, **kwargs):
 
 
 def _unquantify_allclose_arguments(actual, desired, rtol, atol):
-    actual = u.Quantity(actual, subok=True, copy=False)
+    actual = u.Quantity(actual, subok=True, copy=COPY_IF_NEEDED)
 
-    desired = u.Quantity(desired, subok=True, copy=False)
+    desired = u.Quantity(desired, subok=True, copy=COPY_IF_NEEDED)
     try:
         desired = desired.to(actual.unit)
     except u.UnitsError:
@@ -189,7 +190,7 @@ def _unquantify_allclose_arguments(actual, desired, rtol, atol):
         # by default, we assume an absolute tolerance of 0
         atol = u.Quantity(0)
     else:
-        atol = u.Quantity(atol, subok=True, copy=False)
+        atol = u.Quantity(atol, subok=True, copy=COPY_IF_NEEDED)
         try:
             atol = atol.to(actual.unit)
         except u.UnitsError:
@@ -198,7 +199,7 @@ def _unquantify_allclose_arguments(actual, desired, rtol, atol):
                 "are not convertible".format(atol.unit, actual.unit)
             )
 
-    rtol = u.Quantity(rtol, subok=True, copy=False)
+    rtol = u.Quantity(rtol, subok=True, copy=COPY_IF_NEEDED)
     try:
         rtol = rtol.to(u.dimensionless_unscaled)
     except Exception: