Large raster at print size

A synthetic 6000 × 6000 terrain model (36 million pixels, 10 m cells in UTM zone 43N) reduced by fit_raster to the 2000 × 2000 pixels a 90 mm map needs at 600 dpi, drawn with the library's elevation colours (batlow). Synthetic data.
import numpy as np
import rioxarray # noqa: F401 (adds the .rio accessor)
import xarray as xr
import eofigs
n = 6000
y = (np.arange(n, dtype="float32") / n)[:, None]
x = (np.arange(n, dtype="float32") / n)[None, :]
dem = 1200 + 2600 * np.exp(-((x - 0.62) ** 2 + (y - 0.35) ** 2) / 0.05) + 400 * np.sin(9 * x) * y
dem = xr.DataArray(
dem.astype("float32"),
dims=("y", "x"),
coords={"y": 3_900_000 - (np.arange(n) + 0.5) * 10, "x": 750_000 + (np.arange(n) + 0.5) * 10},
name="elevation",
).rio.write_crs("EPSG:32643")
elevation = eofigs.variable("elevation")
eofigs.set_style("elsevier")
fig, ax = eofigs.new_figure("single", height_mm=75)
small = eofigs.fit_raster(dem, width_mm=90, height_mm=75, dpi=600) # 6000 x 6000 -> 2000 x 2000
image = small.plot.imshow(ax=ax, cmap=elevation.cmap, add_colorbar=False)
bar = fig.colorbar(image, ax=ax, label=elevation.label, shrink=0.85)
ax.set_xlabel("Easting (km)")
ax.set_ylabel("Northing (km)")
ax.set_title("")
ax.xaxis.set_major_formatter(lambda v, _: f"{v / 1000:g}")
ax.yaxis.set_major_formatter(lambda v, _: f"{v / 1000:g}")
eofigs.tidy_axes(ax, bar, max_ticks=4)
if __name__ == "__main__":
print(dem.shape, "->", small.shape, small.attrs["eofigs_reduction"])
print(eofigs.save_figure(fig, "large_raster"))
Figure check
Figure check (elsevier): passed
PASS size: 90.0 x 75.0 mm (single width).
PASS resolution: 600 dpi (at least 500 for figures with text and images).
PASS font size: Smallest text 7 pt (minimum 7 pt).
PASS fonts: Fonts in use: Arial.
PASS colour-blind safety: Fewer than two categorical colours.
PASS colour maps: Colour maps: cmc.batlow (sequential).