diff --git a/changelog_unreleased/new-pint.rst b/changelog_unreleased/new-pint.rst
new file mode 100644
index 00000000..f6d93317
--- /dev/null
+++ b/changelog_unreleased/new-pint.rst
@@ -0,0 +1 @@
+* Fixed compatibility with latest `pint` and `pint-xarray` libraries.
diff --git a/docs/usage.ipynb b/docs/usage.ipynb
index 5c9c9795..3b37756e 100644
--- a/docs/usage.ipynb
+++ b/docs/usage.ipynb
@@ -715,7 +715,7 @@
     "import primap2.tests\n",
     "\n",
     "ds = primap2.tests.examples.opulent_ds()\n",
-    "ds[\"CO2\"].pr.loc[{\"product\": \"milk\"}] = np.nan\n",
+    "ds[\"CO2\"].pr.loc[{\"product\": \"milk\"}].pint.magnitude[:] = np.nan\n",
     "\n",
     "ds.pr.coverage(\"product\", \"entity\", \"area\")"
    ]
@@ -966,7 +966,7 @@
     "temp = da.sum(dim=\"area (ISO3)\")\n",
     "temp = temp.expand_dims({\"area (ISO3)\": [\"LATAM\"]})\n",
     "# delete data from the country level for the years 2005-2009 (inclusive)\n",
-    "da.loc[{\"time\": slice(\"2005\", \"2009\")}] = np.nan\n",
+    "da.loc[{\"time\": slice(\"2005\", \"2009\")}].pint.magnitude[:] = np.nan\n",
     "# add regional data to the array\n",
     "da = xr.concat([da, temp], dim=\"area (ISO3)\")\n",
     "da"
@@ -1130,7 +1130,9 @@
    "source": [
     "# delete all data about the years 2005-2009 (inclusive) from the\n",
     "# KYOTOGHG data\n",
-    "ds[\"KYOTOGHG\"].loc[{\"time\": slice(\"2005\", \"2009\")}] = np.nan\n",
+    "ds[\"KYOTOGHG\"].loc[{\"time\": slice(\"2005\", \"2009\")}].pint.magnitude[\n",
+    "    :\n",
+    "] = np.nan\n",
     "ds[\"KYOTOGHG\"]"
    ]
   },
@@ -1176,9 +1178,9 @@
    "source": [
     "# delete all data about the years 2005-2009 from the individual gas data\n",
     "sel = {\"time\": slice(\"2005\", \"2009\")}\n",
-    "ds[\"CO2\"].loc[sel] = np.nan\n",
-    "ds[\"SF6\"].loc[sel] = np.nan\n",
-    "ds[\"CH4\"].loc[sel] = np.nan\n",
+    "ds[\"CO2\"].loc[sel].pint.magnitude[:] = np.nan\n",
+    "ds[\"SF6\"].loc[sel].pint.magnitude[:] = np.nan\n",
+    "ds[\"CH4\"].loc[sel].pint.magnitude[:] = np.nan\n",
     "ds"
    ]
   },
diff --git a/primap2/tests/test_downscale.py b/primap2/tests/test_downscale.py
index f0946f12..087e22eb 100644
--- a/primap2/tests/test_downscale.py
+++ b/primap2/tests/test_downscale.py
@@ -12,7 +12,7 @@
 
 def test_downscale_gas_timeseries(empty_ds):
     for key in empty_ds:
-        empty_ds[key][:] = np.nan
+        empty_ds[key].pint.magnitude[:] = np.nan
     empty_ds["CO2"].loc[{"time": "2002"}] = 1 * ureg("Gg CO2 / year")
     empty_ds["SF6"].loc[{"time": "2002"}] = 1 * ureg("Gg SF6 / year")
     empty_ds["CH4"].loc[{"time": "2002"}] = 1 * ureg("Gg CH4 / year")
@@ -59,7 +59,7 @@ def test_downscale_gas_timeseries(empty_ds):
 
 def test_downscale_timeseries(empty_ds):
     for key in empty_ds:
-        empty_ds[key][:] = np.nan
+        empty_ds[key].pint.magnitude[:] = np.nan
     t = empty_ds.loc[{"area (ISO3)": "BOL"}].copy()
     t["area (ISO3)"] = ["CAMB"]  # here, the sum of COL, ARG, MEX, and BOL
     ds = xr.concat([empty_ds, t], dim="area (ISO3)")
diff --git a/primap2/tests/test_overview.py b/primap2/tests/test_overview.py
index e14c8527..78f243ec 100644
--- a/primap2/tests/test_overview.py
+++ b/primap2/tests/test_overview.py
@@ -58,7 +58,7 @@ def test_array_empty(empty_ds):
 
 def test_array_coverage(empty_ds):
     da = empty_ds["CO2"]
-    da[:] = np.nan
+    da.pint.magnitude[:] = np.nan
     da.name = None
 
     da.pr.loc[{"time": "2001", "area": "COL"}] = 12.0 * ureg("Gg CO2 / year")
@@ -86,7 +86,7 @@ def test_array_coverage(empty_ds):
 def test_array_coverage_multidim(opulent_ds):
     da = opulent_ds["CO2"]
 
-    da.pr.loc[{"product": "milk"}] = np.nan
+    da.pr.loc[{"product": "milk"}].pint.magnitude[:] = np.nan
 
     expected = pd.DataFrame(
         index=da.pr["animal"].values,
@@ -116,7 +116,7 @@ def test_array_coverage_error(opulent_ds):
 
 def test_set_coverage(opulent_ds):
     ds = opulent_ds
-    ds["CO2"].pr.loc[{"product": "milk"}] = np.nan
+    ds["CO2"].pr.loc[{"product": "milk"}].pint.magnitude[:] = np.nan
 
     expected = pd.DataFrame(
         index=ds.pr["product"].values,
@@ -135,7 +135,7 @@ def test_set_coverage(opulent_ds):
 
 def test_set_coverage_entity(opulent_ds):
     ds = opulent_ds
-    ds["CO2"].pr.loc[{"product": "milk"}] = np.nan
+    ds["CO2"].pr.loc[{"product": "milk"}].pint.magnitude[:] = np.nan
 
     expected = pd.DataFrame(
         index=list(ds.keys()),
@@ -165,7 +165,7 @@ def test_set_coverage_boolean(opulent_ds):
 def test_set_coverage_entity_other_dim_not_existing(opulent_ds):
     ds = opulent_ds
 
-    ds["CO2"].pr.loc[{"product": "milk"}] = np.nan
+    ds["CO2"].pr.loc[{"product": "milk"}].pint.magnitude[:] = np.nan
 
     entites_expected = [x for x in ds.keys() if x != "population"]
 
diff --git a/primap2/tests/test_setters.py b/primap2/tests/test_setters.py
index 6d2c8ed5..9e4d9f2c 100644
--- a/primap2/tests/test_setters.py
+++ b/primap2/tests/test_setters.py
@@ -79,7 +79,7 @@ def test_exists_default_error(
     def test_exists_empty_default(
         self, da: xr.DataArray, ts: np.ndarray, co2: pint.Unit, new
     ):
-        da.loc[{"area (ISO3)": "COL"}] = np.nan
+        da.loc[{"area (ISO3)": "COL"}].pint.magnitude[:] = np.nan
         actual = da.pr.set("area", "COL", ts * co2, **new)
         expected = da
         expected.loc[{"area (ISO3)": "COL"}] = ts[..., np.newaxis] * co2
@@ -88,7 +88,7 @@ def test_exists_empty_default(
     def test_exists_somena_default(
         self, da: xr.DataArray, ts: np.ndarray, co2: pint.Unit, new
     ):
-        da.loc[{"area (ISO3)": "COL"}] = np.nan
+        da.loc[{"area (ISO3)": "COL"}].pint.magnitude[:] = np.nan
         da.loc[{"area (ISO3)": "COL", "time": "2001"}] = 2 * co2
         with pytest.raises(
             ValueError,
@@ -125,7 +125,7 @@ def test_exists_fillna(self, da: xr.DataArray, ts: np.ndarray, co2: pint.Unit, n
         assert_aligned_equal(actual, expected)
 
     def test_mixed_default(self, da: xr.DataArray, ts: np.ndarray, co2: pint.Unit):
-        da.loc[{"area (ISO3)": "COL"}] = np.nan
+        da.loc[{"area (ISO3)": "COL"}].pint.magnitude[:] = np.nan
         actual = da.pr.set(
             "area",
             ["COL", "CUB"],
@@ -210,7 +210,7 @@ def test_mixed_from_data_array_overwrite(
             "area", ["CUB", "COL"], 2 * da.pr.loc[{"area": "COL"}], existing="overwrite"
         )
         expected = da.reindex({"area (ISO3)": [*da["area (ISO3)"].values, "CUB"]})
-        expected.loc[{"area (ISO3)": "COL"}] = np.nan
+        expected.loc[{"area (ISO3)": "COL"}].pint.magnitude[:] = np.nan
         expected = expected.fillna(da.loc[{"area (ISO3)": "COL"}] * 2)
         assert_aligned_equal(actual, expected)
 
@@ -428,7 +428,9 @@ def test_existing_overwrite(self, minimal_ds: xr.Dataset, new):
         assert_ds_aligned_equal(actual, expected)
 
     def test_existing_fillna(self, minimal_ds: xr.Dataset, new):
-        minimal_ds["CO2"].pr.loc[{"area": "COL", "time": "2001"}] = np.nan
+        minimal_ds["CO2"].pr.loc[{"area": "COL", "time": "2001"}].pint.magnitude[
+            :
+        ] = np.nan
         actual = minimal_ds.pr.set(
             "area", "COL", minimal_ds.pr.loc[{"area": "MEX"}], existing="fillna", **new
         )
diff --git a/setup.cfg b/setup.cfg
index 086462bd..5789b4df 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -33,8 +33,8 @@ setup_requires =
     setuptools_scm
 install_requires =
     xarray
-    pint>=0.20.1
-    pint_xarray>=0.2
+    pint>=0.23
+    pint_xarray>=0.3
     numpy
     pandas
     openscm_units>=0.5.1