astropy: astropy handles spectral wcs units with logs incorrectly

EDIT: Solution given at https://github.com/astropy/astropy/issues/11215#issuecomment-755906300

It looks either wcslib or astropy’s wrapping of it do not handle the units on the wcs object correctly, and instead override the spectral units to metres, when "WAVE-LOG" is used (note the change to metres in both the output of pixel_to_world and on the w.wcs object):

from astropy.wcs import WCS
import numpy as np
from astropy.units import Unit
w0 = 3.57880000000000E+00
w1 = 1.00000000000000E-04
crval = 10 ** w0
cdelt = crval * w1 * np.log(10)
cunit = Unit('Angstrom')
ctype = "WAVE-LOG"
w = WCS(naxis=1)
w.wcs.crval[0] = crval
w.wcs.cdelt[0] = cdelt
w.wcs.ctype[0] = ctype
w.wcs.cunit[0] = cunit
print(w)
# OUT: WCS Keywords
# OUT: Number of WCS axes: 1
# OUT: CTYPE : 'WAVE-LOG'  
# OUT: CRVAL : 3791.4034418330666  
# OUT: CRPIX : 0.0  
# OUT: PC1_1  : 1.0  
# OUT: CDELT : 0.8730029046691138  
# OUT: NAXIS : 0  0
print(w.wcs.cunit)
# OUT: ['Angstrom']
print(w.pixel_to_world(np.arange(10)))
# OUT: [3792.27644474 3793.14944764 3794.02245055 3794.89545345 3795.76845636
# OUT:  3796.64145926 3797.51446217 3798.38746507 3799.26046798 3800.13347088] m
print(w.wcs.cunit)
# OUT: ['m']
### 

System Details

>>> import platform; print(platform.platform())
Linux-5.10.0-1-amd64-x86_64-with-glibc2.31
>>> import sys; print("Python", sys.version)
Python 3.9.1 (default, Dec  8 2020, 07:51:42)
[GCC 10.2.0]
>>> import numpy; print("Numpy", numpy.__version__)
Numpy 1.19.5
>>> import astropy; print("astropy", astropy.__version__)
astropy 4.2