healpix_analyse.dino ==================== .. py:module:: healpix_analyse.dino .. autoapi-nested-parse:: dino.py ======= DINOv3 (SAT-493M) embeddings of NESTED HEALPix data for the ``healpix_analyse`` package. The idea -------- A ViT is not a "grid" architecture: after the patch embedding it only sees a set of tokens. In the NESTED HEALPix scheme every cell at ``parent_level`` contains exactly ``4**(level - parent_level)`` cells of ``level``, stored contiguously and ordered by a Morton (Z-order) code of their local ``(x, y)`` face coordinates. A NESTED block can therefore be reshaped *exactly* into a square image of side ``2**(level - parent_level)`` pixels, fed to an unmodified DINOv3 network, and every 16x16 patch token maps back to one HEALPix cell at ``level - 4``. This module provides - :func:`nested_to_tiles` NESTED block -> square tiles ``[M, C, S, S]`` - :func:`tiles_to_nested` the exact inverse (for debugging / round trips) - :func:`load_dinov3_sat` load a DINOv3 SAT-493M backbone (hub or HF) - :func:`GetDINOV3SAT` the user-facing function Orientation ----------- Inside a HEALPix base face the local axes ``(x, y)`` are rotated by 45 degrees with respect to the local East/North directions: the "north" corner of a face is at ``(x, y) = (max, max)``. Tiles are built with ``row = S - 1 - y`` and ``col = x`` so that North points to the top-right corner of the image. The orientation is consistent for all tiles of the same base face and rotated for the polar faces, which is a mild nuisance for a satellite backbone (no gravity direction) and should be handled by rotation augmentations when fine-tuning. Dependencies: numpy, torch, ``healpix-geo`` (the GRID4EARTH HEALPix library -- ``healpy`` is deliberately not used anywhere in this package), and either the ``dinov3`` torch-hub repo (``facebookresearch/dinov3``) or ``transformers`` (HF weights). Attributes ---------- .. autoapisummary:: healpix_analyse.dino.SAT493M_MEAN healpix_analyse.dino.SAT493M_STD Classes ------- .. autoapisummary:: healpix_analyse.dino.DINOEmbedding Functions --------- .. autoapisummary:: healpix_analyse.dino.set_ellipsoid healpix_analyse.dino.cell_centres_lonlat healpix_analyse.dino.cells_at_lonlat healpix_analyse.dino.cell_vectors healpix_analyse.dino.cell_neighbours healpix_analyse.dino.nested_to_xy healpix_analyse.dino.xy_to_nested healpix_analyse.dino.nested_to_tiles healpix_analyse.dino.tiles_to_nested healpix_analyse.dino.tile_grid_cell_ids healpix_analyse.dino.healpix_gsd_m healpix_analyse.dino.tangent_grid_lonlat healpix_analyse.dino.offsets_to_lonlat healpix_analyse.dino.lonlat_to_offsets healpix_analyse.dino.parent_cover_px healpix_analyse.dino.sample_healpix healpix_analyse.dino.pow2_cover_px healpix_analyse.dino.tangent_tiles healpix_analyse.dino.load_dinov3_sat healpix_analyse.dino.GetDINOV3SAT Module Contents --------------- .. py:data:: SAT493M_MEAN :type: Tuple[float, float, float] :value: (0.43, 0.411, 0.296) .. py:data:: SAT493M_STD :type: Tuple[float, float, float] :value: (0.213, 0.156, 0.143) .. py:function:: set_ellipsoid(name: str) -> str Choose the reference ellipsoid ("sphere", "wgs84", ...) and return the old one. Call it once, with whatever the store's ``dggs.ellipsoid.name`` says, before any embedding is computed. .. py:function:: cell_centres_lonlat(level: int, cell_id: numpy.ndarray) -> Tuple[numpy.ndarray, numpy.ndarray] Centres of NESTED cells at ``level``, as (lon, lat) in degrees. .. py:function:: cells_at_lonlat(level: int, lon: numpy.ndarray, lat: numpy.ndarray) -> numpy.ndarray NESTED cell ids at ``level`` containing the given (lon, lat), in degrees. .. py:function:: cell_vectors(level: int, cell_id: numpy.ndarray) -> numpy.ndarray Unit vectors of the cell centres, shape ``[3, N]``. .. py:function:: cell_neighbours(level: int, cell_id: numpy.ndarray) -> numpy.ndarray The eight NESTED neighbours of each cell, shape ``[8, N]``; -1 where absent. .. py:function:: nested_to_xy(rel_id: numpy.ndarray) -> Tuple[numpy.ndarray, numpy.ndarray] Decode a *relative* NESTED index into local ``(x, y)`` coordinates. ``rel_id`` is the index of a cell inside its ancestor block (``cell_id - parent_id * 4**k``); the result satisfies ``0 <= x, y < 2**k``. This is the HEALPix ``nest2xyf`` convention: ``x`` occupies the even bits, ``y`` the odd bits. .. py:function:: xy_to_nested(x: numpy.ndarray, y: numpy.ndarray) -> numpy.ndarray Inverse of :func:`nested_to_xy`. .. py:function:: nested_to_tiles(data: ArrayLike, cell_id: ArrayLike, level: int, parent_level: int, *, fill: str = 'mean', duplicates: str = 'mean') -> Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray] Reorder a NESTED HEALPix dataset into square tiles, one per ``parent_level`` cell present in ``cell_id``. :Parameters: * **data** (:py:class:`array [N`, :py:class:`C]` or :py:class:`[N]`) -- Pixel values at HEALPix ``level``. NaN marks a missing pixel. * **cell_id** (:py:class:`int array [N]`) -- NESTED cell ids at ``level``. Duplicates are not allowed. * **level** (:py:class:`int`) -- Resolution of ``data`` (``nside = 2**level``). * **parent_level** (:py:class:`int`) -- Resolution of the tiles. ``k = level - parent_level >= 1``; the tiles have side ``S = 2**k``. * **fill** (``{"mean", "zero", "nan"}``) -- Value given to the pixels of a tile that are absent from ``cell_id`` (or NaN): the per-tile, per-band mean of the available pixels, zero, or NaN. * **duplicates** (``{"mean", "first", "error"}``) -- What to do when the same cell id appears several times, which happens when data projected from another grid puts two source pixels in the same HEALPix cell: average them (NaN-aware), keep the first one, or raise. :returns: * **tiles** (:py:class:`float32 array [M`, :py:class:`C`, :py:class:`S`, :py:class:`S]`) * **parent_ids** (:py:class:`int64 array [M]`) -- NESTED ids of the tiles at ``parent_level`` (sorted). * **valid** (:py:class:`bool array [M`, :py:class:`S`, :py:class:`S]`) -- Availability mask of every pixel. * **coverage** (:py:class:`float array [M]`) -- Fraction of available pixels per tile. .. py:function:: tiles_to_nested(tiles: numpy.ndarray, parent_ids: numpy.ndarray, parent_level: int, level: int) -> Tuple[numpy.ndarray, numpy.ndarray] Inverse of :func:`nested_to_tiles` (full tiles, no mask). :Parameters: * **tiles** (:py:class:`array [M`, :py:class:`C`, :py:class:`S`, :py:class:`S]`) * **parent_ids** (:py:class:`int array [M]`) * **parent_level, level** (:py:class:`int`) :returns: * **data** (:py:class:`array [M * S * S`, :py:class:`C]`) * **cell_id** (:py:class:`int64 array [M * S * S]`) -- NESTED ids at ``level``, in tile-major order. .. py:function:: tile_grid_cell_ids(parent_ids: numpy.ndarray, parent_level: int, sub_level: int) -> numpy.ndarray NESTED ids (at ``sub_level``) of the ``[S, S]`` image grid of every tile. Returns an int64 array ``[M, S, S]`` indexed as ``[tile, row, col]`` with the same orientation convention as :func:`nested_to_tiles`. .. py:function:: healpix_gsd_m(level: int, radius: float = EARTH_RADIUS_M) -> float Square root of the HEALPix cell area at ``level``, in metres. .. py:function:: tangent_grid_lonlat(centre_lon: numpy.ndarray, centre_lat: numpy.ndarray, size: int, gsd_rad: float, *, shift: Tuple[float, float] = (0.0, 0.0), radius: float = EARTH_RADIUS_M) -> Tuple[numpy.ndarray, numpy.ndarray] Longitude / latitude of north-up gnomonic grids centred on given points. :Parameters: * **centre_lon, centre_lat** (:py:class:`array [M]`) -- Tangent points, in degrees. * **size** (:py:class:`int` or :py:class:`(int`, :py:class:`int)`) -- Side of the grid in pixels; a pair gives (rows, columns). * **gsd_rad** (:py:class:`float`) -- Pixel spacing in tangent units (ground sampling / Earth radius). * **shift** (:py:class:`(float`, :py:class:`float)`) -- Extra offset of the grid origin, in pixels (row, column). Used by the sliding-window oversampling. :returns: **lon, lat** (:py:class:`array [M`, :py:class:`rows`, :py:class:`cols]`) -- Row 0 is the northernmost, column 0 the westernmost: the images are north-up and east-right, whatever the HEALPix face. .. py:function:: offsets_to_lonlat(centre_lon: numpy.ndarray, centre_lat: numpy.ndarray, xi: numpy.ndarray, eta: numpy.ndarray) -> Tuple[numpy.ndarray, numpy.ndarray] Gnomonic offsets (east ``xi``, north ``eta``, in tangent units) to lon/lat. ``xi`` and ``eta`` are broadcast against the ``[M]`` tangent points, giving ``[M, *xi.shape]`` outputs. .. py:function:: lonlat_to_offsets(centre_lon: numpy.ndarray, centre_lat: numpy.ndarray, lon: numpy.ndarray, lat: numpy.ndarray) -> Tuple[numpy.ndarray, numpy.ndarray] Inverse of :func:`offsets_to_lonlat`: lon/lat to gnomonic offsets. ``centre_*`` has shape ``[M]`` and ``lon`` / ``lat`` shape ``[M, ...]``; the returned ``xi`` (east) and ``eta`` (north) have the shape of ``lon``. .. py:function:: parent_cover_px(centre_ids: numpy.ndarray, parent_level: int, gsd_m: float, *, margin_px: int = 8, radius: float = EARTH_RADIUS_M) -> Tuple[int, int] Size (rows, columns) of the smallest tangent image covering every parent cell. A HEALPix cell is a parallelogram, not a square, so a north-up image that must contain it entirely is larger than the cell -- and rectangular. Both sides are rounded up to a multiple of the DINOv3 patch. .. py:function:: sample_healpix(data: numpy.ndarray, cell_id: numpy.ndarray, level: int, lon: numpy.ndarray, lat: numpy.ndarray, *, interpolation: str = 'bilinear') -> Tuple[numpy.ndarray, numpy.ndarray] Sample a partial-sky NESTED HEALPix map at arbitrary lon/lat. Cells absent from ``cell_id`` and non-finite values are treated as missing: the interpolation weights are renormalised over the available neighbours, and the returned mask is False where nothing was available. :returns: * **values** (:py:class:`float32 array [...`, :py:class:`C]`) * **valid** (:py:class:`bool array [...]`) .. py:function:: pow2_cover_px(centre_ids: numpy.ndarray, parent_level: int, gsd_m: float, *, margin_px: int = 8, radius: float = EARTH_RADIUS_M) -> int Side of the smallest **square power-of-two** image containing every parent cell. A HEALPix cell and a square of the same area cannot contain one another: a square of side ``2**k`` px has exactly the cell's area, so roughly a fifth of the cell always falls outside it. Those cells then have no token of their own. Doubling the side fixes it and keeps the shape a ViT expects. .. py:function:: tangent_tiles(data: ArrayLike, cell_id: ArrayLike, level: int, parent_level: int, *, tile_px: Optional[Union[int, str, Tuple[int, int]]] = None, gsd_m: Optional[float] = None, shift: Tuple[float, float] = (0.0, 0.0), interpolation: str = 'bilinear', fill: str = 'mean', duplicates: str = 'mean', radius: float = EARTH_RADIUS_M) -> Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray] Resample the data onto north-up tangent-plane images, one per tile centre. Tile centres are the HEALPix cells of ``parent_level`` present in the data; each image is a square of ``tile_px`` pixels at a constant ground sampling, so its shape is right whatever the local HEALPix distortion. Unlike :func:`nested_to_tiles` the images may overlap and are not restricted to the content of one NESTED block. :returns: * **tiles** (:py:class:`float32 [M`, :py:class:`C`, :py:class:`S`, :py:class:`S]`) * **centre_ids** (:py:class:`int64 [M] NESTED ids` of :py:class:`the tile centres`) * **valid** (:py:class:`bool [M`, :py:class:`S`, :py:class:`S]`) * **coverage** (:py:class:`float [M]`) * **lon, lat** (:py:class:`float64 [M`, :py:class:`S`, :py:class:`S] geographic position` of :py:class:`every pixel`) .. py:function:: load_dinov3_sat(model_name: str = 'dinov3_vitl16', weights: Optional[str] = None, *, source: str = 'auto', repo: str = 'facebookresearch/dinov3', device: Optional[Union[str, torch.device]] = None) -> torch.nn.Module Load a DINOv3 SAT-493M backbone in evaluation mode. :Parameters: * **model_name** (``{"dinov3_vitl16", "dinov3_vit7b16"}``) -- SAT-493M is released for ViT-L/16 and ViT-7B/16 only. * **weights** (:py:class:`str`, *optional*) -- Path to the ``.pth`` checkpoint (torch-hub route; DINOv3 weights are gated and must be downloaded after accepting the licence) or a Hugging Face model id / local directory (HF route). * **source** (``{"auto", "hub", "hf"}``) -- "hub" uses ``torch.hub.load(repo, model_name, weights=weights)``; "hf" uses ``transformers.AutoModel``; "auto" tries hub then hf. * **repo** (:py:class:`str`) -- Torch-hub repo or a local clone of ``facebookresearch/dinov3``. * **device** (:py:class:`str` or :py:class:`torch.device`, *optional*) .. py:class:: DINOEmbedding Result of :func:`GetDINOV3SAT`. :ivar embedding: One vector per tile of ``parent_level`` (CLS token by default). :vartype embedding: :py:class:`array [M`, :py:class:`Ndino]` :ivar cell_id: NESTED ids of the tiles at ``parent_level``. :vartype cell_id: :py:class:`int64 array [M]` :ivar parent_level: :vartype parent_level: :py:class:`int` :ivar coverage: Fraction of ``level`` pixels available in each tile. :vartype coverage: :py:class:`float array [M]` :ivar patch_embedding: Patch tokens (``return_patches=True``); one per HEALPix cell at ``patch_level = level - 4``. :vartype patch_embedding: :py:class:`array [M * P`, :py:class:`Ndino]` or :py:obj:`None` :ivar patch_cell_id: :vartype patch_cell_id: :py:class:`int64 array [M * P]` or :py:obj:`None` :ivar patch_level: :vartype patch_level: :py:class:`int` or :py:obj:`None` .. py:attribute:: embedding :type: numpy.ndarray .. py:attribute:: cell_id :type: numpy.ndarray .. py:attribute:: parent_level :type: int .. py:attribute:: coverage :type: numpy.ndarray .. py:attribute:: patch_embedding :type: Optional[numpy.ndarray] :value: None .. py:attribute:: patch_cell_id :type: Optional[numpy.ndarray] :value: None .. py:attribute:: patch_level :type: Optional[int] :value: None .. py:attribute:: patch_lon :type: Optional[numpy.ndarray] :value: None .. py:attribute:: patch_lat :type: Optional[numpy.ndarray] :value: None .. py:attribute:: projection :type: str :value: 'tangent' .. py:attribute:: over_sample :type: int :value: 1 .. py:attribute:: gsd_m :type: Optional[float] :value: None .. py:attribute:: tile_px :type: Optional[int] :value: None .. py:function:: GetDINOV3SAT(data: ArrayLike, cell_id: ArrayLike, level: int, parent_level: int, *, projection: str = 'tangent', over_sample: int = 1, context_px: Optional[int] = None, out_cells: Optional[ArrayLike] = None, gsd_m: Optional[float] = None, tile_px: Optional[int] = None, interpolation: str = 'bilinear', model: Optional[torch.nn.Module] = None, weights: Optional[str] = None, model_name: str = 'dinov3_vitl16', bands: Sequence[int] = (0, 1, 2), mean: Sequence[float] = SAT493M_MEAN, std: Sequence[float] = SAT493M_STD, pooling: str = 'cls', return_patches: bool = False, min_coverage: float = 0.0, fill: str = 'mean', duplicates: str = 'mean', batch_size: int = 16, device: Optional[Union[str, torch.device]] = None, autocast: bool = True) -> DINOEmbedding DINOv3 SAT-493M embeddings of NESTED HEALPix data. Every HEALPix cell of ``parent_level`` present in ``cell_id`` gives one square image, which is passed to an unmodified DINOv3 backbone. How that image is built is the important choice: ``projection="tangent"`` (default) The image is resampled onto a north-up gnomonic plane tangent at the cell centre, with a constant ground sampling. Shapes are then correct: HEALPix cells have equal area but not equal shape, and at 52 deg N the two NESTED axes measure 10.6 m and 16.7 m with a 60.5 deg angle between them, so the "square" of the ``nested`` mode is really a parallelogram stretched by a factor 2. A network trained on map-projected imagery cannot be expected to see through that. The image is sized to contain the whole (parallelogram-shaped) cell, so it is rectangular and a little larger than the cell. The tokens sit on a square metric grid while the cells are parallelograms, so the output is *not* built by assigning each token to the cell it falls in -- that would leave some cells with two tokens and others with none, a moire of holes on a map. Instead the cells of the token level inside each tile are enumerated and the token field is read at their positions: every cell carries exactly one embedding, and the output tiles the sphere without holes. ``patch_lon`` / ``patch_lat`` are then the cell centres. The price is a resampling, twice. ``projection="percell"`` One tangent plane **per output cell**: for every cell of ``parent_level`` a north-up image of ``context_px`` pixels is built, tangent at that cell's centre, and the embedding is the patch token holding the centre (or ``pooling="cls"`` for the whole window). The window is centred exactly on the cell, in the cell's own frame, so nothing is interleaved and nothing is resampled back -- this is the exact sliding window, and ``level - parent_level >= 4`` no longer applies since the cell may be finer than a patch. It costs one forward pass per cell, but each window is small: with the default ``context_px`` -- the cell's own footprint -- the total number of tokens equals the number of output cells, exactly what the tiled modes produce, so this is affordable. What is expensive is context: a window ``m`` patches wide costs ``m**2`` tokens per cell. ``projection="nested"`` The historical mode: the NESTED block is reshaped into an image, which is exact in index space and free of resampling, and every token maps onto exactly one cell of ``level - 4``. Geometrically wrong away from the equator; kept for comparison. :Parameters: * **data** (:py:class:`array [N`, :py:class:`C]`) -- Reflectances at ``level``, expected in ``[0, 1]`` (Sentinel-2 ``DN / 10000`` for instance). NaN marks a missing pixel. * **cell_id** (:py:class:`int array [N]`) -- NESTED cell ids at ``level``. * **level** (:py:class:`int`) -- Resolution of ``data``. * **parent_level** (:py:class:`int`) -- Resolution of the tile centres; ``level - parent_level >= 4`` is required (one DINO patch is 16 px = 4 HEALPix levels). ``level - 8`` (256 px tiles, 16x16 patch tokens) is the natural choice. * **projection** (``{"tangent", "nested"}``) -- See above. * **context_px** (:py:class:`int`, *optional*) -- Side of the window in ``percell`` mode, a multiple of 16. The default is the cell's own footprint, so the total number of tokens equals the number of output cells -- the same count as the tiled modes. A larger window gives the network context, at a cost growing as ``(context_px / cell)**2``: for cells of 16 px, 48 px costs 9 times more and 224 px 196 times. With a centred pooling an even number of patches per axis leaves no token on the centre, so the value is raised by one patch. * **out_cells** (:py:class:`int array`, *optional*) -- Restrict ``percell`` mode to these cells instead of every cell of ``parent_level`` present in the data. The way to try it on a small area before paying for the whole scene. * **over_sample** (:py:class:`int`) -- Sliding-window factor, a power of two. ``1`` (default) runs the patches side by side, exactly as the network does. ``n`` runs the network ``n**2`` times on windows shifted by ``16 / n`` pixels and interleaves the tokens, giving a token field ``n`` times denser in each direction, at ``level - 4 + log2(n)``. Cost grows as ``n**2``. * **gsd_m** (:py:class:`float`, *optional*) -- Ground sampling of the tangent grid, in metres. Defaults to the HEALPix native value ``sqrt(cell area)`` at ``level``. * **tile_px** (:py:class:`int` or :py:class:`(int`, :py:class:`int)`, *optional*) -- Size of the images, ``(rows, columns)``. In tangent mode the default is the smallest multiple of 16 containing every parent cell (see :func:`parent_cover_px`); in nested mode it is ``2**(level - parent_level)``. * **interpolation** (``{"bilinear", "nearest"}``) -- Resampling onto the tangent grid. * **model** (:py:class:`nn.Module`, *optional*) -- A backbone already loaded with :func:`load_dinov3_sat`. When ``None`` the model is loaded from ``weights`` / ``model_name``. * **weights, model_name** -- Passed to :func:`load_dinov3_sat` when ``model`` is ``None``. * **bands** (:py:class:`sequence` of :py:class:`3 int`) -- Indices of the (R, G, B) bands in ``data`` -- Sentinel-2 B04, B03, B02. * **mean, std** (:py:class:`sequence` of :py:class:`3 float`) -- Normalisation constants (SAT-493M defaults). * **pooling** (``{"cls", "mean", "cls+mean", "centre"}``) -- Tile vector: CLS token, mean of the patch tokens, or their concatenation (``2 * Ndino``). In ``percell`` mode the default is the patch token holding the cell centre (``"centre"``), which is what describes the cell itself rather than its whole window. * **return_patches** (:py:class:`bool`) -- Also return the patch tokens. * **min_coverage** (:py:class:`float`) -- Drop tiles whose fraction of usable pixels is below this value; 1.0 keeps only the complete ones. * **fill** (``{"mean", "zero"}``) -- How missing pixels are filled before normalisation. * **duplicates** (``{"mean", "first", "error"}``) -- How repeated cell ids are aggregated (see :func:`nested_to_tiles`). * **batch_size** (:py:class:`int`) * **device** (:py:class:`str` or :py:class:`torch.device`, *optional*) * **autocast** (:py:class:`bool`) -- Use bfloat16 autocast on CUDA. :returns: :py:class:`DINOEmbedding`