f_mosaic = 't1.nc'

meta = {'width': dat_f.shape[1],
        'height': dat_f.shape[2],
        'crs': rasterio.crs.CRS(init='epsg:'+fi['CPER']['Reflectance']['Metadata']['Coordinate_System']['EPSG Code'].value.decode("utf-8")),
        'transform': aff_final,
        'count': dat_f.shape[0]}

with netCDF4.Dataset(f_mosaic, mode='w', format="NETCDF4") as t1:
    # Create spatial dimensions
    y = t1.createDimension('y', meta['width'])
    x = t1.createDimension('x', meta['height'])
    wl_dim = t1.createDimension('wl',meta['count'])
    reflectance = t1.createVariable("reflectance","int16",("wl","y","x",),fill_value=null_val,zlib=True)
    reflectance.setncattr('grid_mapping', 'crs')
    crs = t1.createVariable('crs', 'c')
    crs.spatial_ref = meta['crs'].wkt
    crs.epsg_code = meta['crs'].to_string()
    crs.GeoTransform = " ".join(str(x) for x in meta['transform'].to_gdal())

dat_f.to_netcdf(path=f_mosaic,mode='a',format='NETCDF4',encoding={'reflectance':{'zlib':True}})