bug with to_netcdf when reading dataset via THREDDS

I am trying to extract data for a specific station from data observations accessed via THREDDS. At first I got an error but a colleague found a way to make it work. In the code below it works if APPLY_PATCH=True but it crashes if APPLY_PATCH=False.

I do not know if this behaviour is expected and/or linked to the way the data is 'served' via THREDDS but I would never have found the way around it by myself.

import xarray as xr

CLEF_STATION = "MONTREAL"
NOM_VAR = 'prsn'
APPLY_PATCH = False

url_pavics = "http://pavics.ouranos.ca:8083/twitcher/ows/proxy/thredds/dodsC"
url_donnees_homo = "/".join((url_pavics, "datasets/station_obs/ECCC_AHCCD_gen2_precipitation.ncml"))

ds = xr.open_dataset(url_donnees_homo)

# data for specific station(s)
dsr = ds.sel(station=ds.station_name.str.contains(CLEF_STATION))

# apply patch
if APPLY_PATCH:
    for var in dsr.variables:
        if dsr[var].dtype == 'O':
            dsr[var].encoding = {}

# save to a netcdf file
dsr.to_netcdf("out.nc", format="NETCDF4_CLASSIC")

What happened:

When running this I get the following error

runfile('/Users/sbiner/cloud/OneDrive - Ouranos inc/job/stagiaires_etudiants/2021_tempete_neige/code/issue_xaray.py', wdir='/Users/sbiner/cloud/OneDrive - Ouranos inc/job/stagiaires_etudiants/2021_tempete_neige/code')
Traceback (most recent call last):

  File "src/netCDF4/_netCDF4.pyx", line 4931, in netCDF4._netCDF4.Variable.__setitem__

ValueError: cannot reshape array of size 7 into shape (1,1)


During handling of the above exception, another exception occurred:

Traceback (most recent call last):

  File "/Users/sbiner/cloud/OneDrive - Ouranos inc/job/stagiaires_etudiants/2021_tempete_neige/code/issue_xaray.py", line 14, in <module>
    dsr.to_netcdf("out.nc", format="NETCDF4_CLASSIC")

  File "/Users/sbiner/anaconda3/envs/py38/lib/python3.8/site-packages/xarray/core/dataset.py", line 1689, in to_netcdf
    return to_netcdf(

  File "/Users/sbiner/anaconda3/envs/py38/lib/python3.8/site-packages/xarray/backends/api.py", line 1107, in to_netcdf
    dump_to_store(

  File "/Users/sbiner/anaconda3/envs/py38/lib/python3.8/site-packages/xarray/backends/api.py", line 1154, in dump_to_store
    store.store(variables, attrs, check_encoding, writer, unlimited_dims=unlimited_dims)

  File "/Users/sbiner/anaconda3/envs/py38/lib/python3.8/site-packages/xarray/backends/common.py", line 255, in store
    self.set_variables(

  File "/Users/sbiner/anaconda3/envs/py38/lib/python3.8/site-packages/xarray/backends/common.py", line 297, in set_variables
    writer.add(source, target)

  File "/Users/sbiner/anaconda3/envs/py38/lib/python3.8/site-packages/xarray/backends/common.py", line 146, in add
    target[...] = source

  File "/Users/sbiner/anaconda3/envs/py38/lib/python3.8/site-packages/xarray/backends/netCDF4_.py", line 64, in __setitem__
    data[key] = value

  File "src/netCDF4/_netCDF4.pyx", line 4933, in netCDF4._netCDF4.Variable.__setitem__

  File "<__array_function__ internals>", line 5, in broadcast_to

  File "/Users/sbiner/anaconda3/envs/py38/lib/python3.8/site-packages/numpy/lib/stride_tricks.py", line 411, in broadcast_to
    return _broadcast_to(array, shape, subok=subok, readonly=True)

  File "/Users/sbiner/anaconda3/envs/py38/lib/python3.8/site-packages/numpy/lib/stride_tricks.py", line 348, in _broadcast_to
    it = np.nditer(

ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (1,7)  and requested shape (1,1)

What you expected to happen:

I expect a netcdf file to be written

Environment:

Start with Dataset.to_netcdf and the backend paths shown in the traceback, especially xarray/core/backends/common.py and xarray/backends/netCDF4_.py. Reproduce the THREDDS selection and write operation, then identify why the object-valued variables produce a shape mismatch. Done means the provided dataset can be written without applying the workaround and a regression test covers the behavior.