why find_optimal_celestial_wcs() always return PC coefficients with 1.0, 0.0, 1.0, 0.0

[Eririf]

I tried several images but find_optimal_celestial_wcs() keeps returning PC1_1, PC2_2 = 1.0, PC1_2, PC2_1 = 0. I want to write it to the final coadd image.

files = glob.glob(path + '*.FITS')
hdus = [fits.open(file)[0] for file in files[:]]
from reproject.mosaicking import find_optimal_celestial_wcs
wcs_out, shape_out = find_optimal_celestial_wcs(hdus, auto_rotate=True,frame='icrs')

WCSAXES =                    2 / Number of coordinate axes                      
CRPIX1  =      4169.3332638275 / Pixel coordinate of reference point            
CRPIX2  =      4087.4118749305 / Pixel coordinate of reference point            
PC1_1   =     0.99925654189794 / Coordinate transformation matrix element       
PC1_2   =    0.038553384730345 / Coordinate transformation matrix element       
PC2_1   =   -0.038553384730345 / Coordinate transformation matrix element       
PC2_2   =     0.99925654189794 / Coordinate transformation matrix element       
CDELT1  = -0.00046184030975759 / [deg] Coordinate increment at reference point  
CDELT2  =  0.00046184030975759 / [deg] Coordinate increment at reference point  
CUNIT1  = 'deg'                / Units of coordinate increment and value        
CUNIT2  = 'deg'                / Units of coordinate increment and value        
CTYPE1  = 'RA---TAN'           / Right ascension, gnomonic projection           
CTYPE2  = 'DEC--TAN'           / Declination, gnomonic projection               
CRVAL1  =      134.94692177133 / [deg] Coordinate value at reference point      
CRVAL2  =      38.095156010732 / [deg] Coordinate value at reference point      
LONPOLE =                180.0 / [deg] Native longitude of celestial pole       
LATPOLE =      38.095156010732 / [deg] Native latitude of celestial pole        
MJDREF  =                  0.0 / [d] MJD of fiducial time                       
RADESYS = 'ICRS'               / Equatorial coordinate system     

[Eririf]

I suggest to change from

        # Set rotation matrix (use PC instead of CROTA2 since PC is the
        # recommended approach)
        pc = np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]])
        wcs_final.wcs.pc = pc

to

       wcs_final.wcs.pc = pc * resolution