WCS get_pc() not taking CROTA into account

[astrofrog]

According to the docstring of get_pc, this method should return a matrix that takes into account CROTA if present:

Returns the PC matrix in read-only form. Unlike the
~astropy.wcs.Wcsprm.pc property, this works even when the header
specifies the linear transformation matrix in one of the alternative
CDi_ja or CROTAia forms. This is useful when you want access
to the linear transformation matrix, but don't care how it was
specified in the header.

but this is not the case:

In [21]: wcs = WCS(naxis=2)

In [22]: wcs.wcs.cdelt = [1,2]

In [23]: wcs.wcs.crota = [0, 60]

In [24]: wcs.wcs.get_pc()
Out[24]: 
array([[ 1.,  0.],
       [ 0.,  1.]])

Also, to_header outputs a header without CROTA2.

cc @nden