GDAL with NetCDF

import numpy as np
from osgeo import gdal

gdal.UseExceptions()

def load_data():
    src = "NETCDF:earthaccess_data/20020601090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc: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()