GDAL with HDF5 and VRT

import numpy as np
from osgeo import gdal

gdal.UseExceptions()

def load_data():
    src = "vrt://earthaccess_data/20020601090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc?a_ullr=-179.995,89.995,180.005,-89.995&sd_name=analysed_sst"
    ds = gdal.Open(src)
    band = ds.GetRasterBand(1)
    arr = band.ReadAsArray().astype("float32", casting="unsafe")
    mask = band.GetMaskBand().ReadAsArray()
    ma = np.ma.masked_array(arr, np.logical_not(mask), fill_value=np.nan)
    np.multiply(ma, band.GetScale(), out=ma, casting="unsafe")
    np.add(ma, band.GetOffset(), out=ma, casting="unsafe")
    return ma.filled(fill_value=np.nan)

if __name__ == "__main__":
    da = load_data()