diff --git a/changes/1585.mosaic_pipeline.rst b/changes/1585.mosaic_pipeline.rst
new file mode 100644
index 000000000..977888ff6
--- /dev/null
+++ b/changes/1585.mosaic_pipeline.rst
@@ -0,0 +1 @@
+Roundtrip L3 wcsinfo especially when skycell specifications are used
diff --git a/romancal/pipeline/mosaic_pipeline.py b/romancal/pipeline/mosaic_pipeline.py
index 7294ef3dd..fb5f2315f 100644
--- a/romancal/pipeline/mosaic_pipeline.py
+++ b/romancal/pipeline/mosaic_pipeline.py
@@ -13,6 +13,7 @@
 from astropy import units as u
 from astropy.modeling import models
 from gwcs import WCS, coordinate_frames
+from stcal.alignment import util as wcs_util
 
 import romancal.datamodels.filetype as filetype
 from romancal.datamodels import ModelLibrary
@@ -114,13 +115,13 @@ def process(self, input):
 
                     # check to see if there exists a skycell on disk if not create it
                     if not isfile(skycell_file_name):
-                        # extract the wcs info from the record for generate_tan_wcs
+                        # extract the wcs info from the record for skycell_to_wcs
                         log.info(
                             "Creating skycell image at ra: %f  dec %f",
                             float(skycell_record["ra_center"]),
                             float(skycell_record["dec_center"]),
                         )
-                        skycell_wcs = generate_tan_wcs(skycell_record)
+                        skycell_wcs = skycell_to_wcs(skycell_record)
                         # skycell_wcs.bounding_box = bounding_box
 
                         # For resample to use an external grid we need to pass it the skycell gwcs object
@@ -128,6 +129,7 @@ def process(self, input):
                         wcs_tree = {"wcs": skycell_wcs}
                         wcs_file = asdf.AsdfFile(wcs_tree)
                         wcs_file.write_to("skycell_wcs.asdf")
+
                         self.resample.output_wcs = "skycell_wcs.asdf"
                         self.resample.output_shape = (
                             int(skycell_record["nx"]),
@@ -162,41 +164,96 @@ def process(self, input):
         return result
 
 
-def generate_tan_wcs(skycell_record):
-    """extract the wcs info from the record for generate_tan_wcs
-    we need the scale, ra, dec, bounding_box"""
+def skycell_to_wcs(skycell_record):
+    """From a skycell record, generate a GWCS
+
+    Parameters
+    ----------
+    skycell_record : dict
+        A skycell record, or row, from the skycell patches table.
+
+    Returns
+    -------
+    wcsobj : wcs.GWCS
+        The GWCS object from the skycell record.
+    """
+    wcsinfo = dict()
+
+    # The scale is given in arcseconds per pixel. Convert to degrees.
+    wcsinfo["pixel_scale"] = float(skycell_record["pixel_scale"]) / 3600.0
+
+    # Remaining components of the wcsinfo block
+    wcsinfo["ra_ref"] = float(skycell_record["ra_projection_center"])
+    wcsinfo["dec_ref"] = float(skycell_record["dec_projection_center"])
+    wcsinfo["x_ref"] = float(skycell_record["x0_projection"])
+    wcsinfo["y_ref"] = float(skycell_record["y0_projection"])
+    wcsinfo["orientat"] = float(skycell_record["orientat_projection_center"])
+    wcsinfo["rotation_matrix"] = None
 
-    scale = float(skycell_record["pixel_scale"])
-    ra_center = float(skycell_record["ra_projection_center"])
-    dec_center = float(skycell_record["dec_projection_center"])
-    shiftx = float(skycell_record["x0_projection"])
-    shifty = float(skycell_record["y0_projection"])
+    # Bounding box of the skycell. Note that the center of the pixels are at (0.5, 0.5)
     bounding_box = (
         (-0.5, -0.5 + skycell_record["nx"]),
         (-0.5, -0.5 + skycell_record["ny"]),
     )
 
-    # components of the model
-    # shift = models.Shift(shiftx) & models.Shift(shifty)
+    wcsobj = wcsinfo_to_wcs(wcsinfo, bounding_box=bounding_box)
+    return wcsobj
+
+
+def wcsinfo_to_wcs(wcsinfo, bounding_box=None, name="wcsinfo"):
+    """Create a GWCS from the L3 wcsinfo meta
+
+    Parameters
+    ----------
+    wcsinfo : dict or MosaicModel.meta.wcsinfo
+        The L3 wcsinfo to create a GWCS from.
 
-    # select a scale for the skycell image, this will come from INS and may
-    # be optimized for the different survey programs
-    scale_x = scale
-    scale_y = scale
-    # set the pixelsscale to 0.1 arcsec/pixel
-    pixelscale = models.Scale(scale_x / 3600.0) & models.Scale(scale_y / 3600.0)
+    bounding_box : None or 4-tuple
+        The bounding box in detector/pixel space. Form of input is:
+        (x_left, x_right, y_bottom, y_top)
 
-    pixelshift = models.Shift(-1.0 * shiftx) & models.Shift(-1.0 * shifty)
+    name : str
+        Value of the `name` attribute of the GWCS object.
+
+    Returns
+    -------
+    wcs : wcs.GWCS
+        The GWCS object created.
+    """
+    pixelshift = models.Shift(-wcsinfo["x_ref"], name="crpix1") & models.Shift(
+        -wcsinfo["y_ref"], name="crpix2"
+    )
+    pixelscale = models.Scale(wcsinfo["pixel_scale"], name="cdelt1") & models.Scale(
+        wcsinfo["pixel_scale"], name="cdelt2"
+    )
     tangent_projection = models.Pix2Sky_TAN()
-    celestial_rotation = models.RotateNative2Celestial(ra_center, dec_center, 180.0)
-    det2sky = pixelshift | pixelscale | tangent_projection | celestial_rotation
+    celestial_rotation = models.RotateNative2Celestial(
+        wcsinfo["ra_ref"], wcsinfo["dec_ref"], 180.0
+    )
+
+    matrix = wcsinfo.get("rotation_matrix", None)
+    if matrix:
+        matrix = np.array(matrix)
+    else:
+        orientat = wcsinfo.get("orientat", 0.0)
+        matrix = wcs_util.calc_rotation_matrix(
+            np.deg2rad(orientat), v3i_yangle=0.0, vparity=1
+        )
+        matrix = np.reshape(matrix, (2, 2))
+    rotation = models.AffineTransformation2D(matrix, name="pc_rotation_matrix")
+    det2sky = (
+        pixelshift | rotation | pixelscale | tangent_projection | celestial_rotation
+    )
+
     detector_frame = coordinate_frames.Frame2D(
         name="detector", axes_names=("x", "y"), unit=(u.pix, u.pix)
     )
     sky_frame = coordinate_frames.CelestialFrame(
         reference_frame=coordinates.ICRS(), name="icrs", unit=(u.deg, u.deg)
     )
-    wcsobj = WCS([(detector_frame, det2sky), (sky_frame, None)])
-    wcsobj.bounding_box = bounding_box
+    wcsobj = WCS([(detector_frame, det2sky), (sky_frame, None)], name=name)
+
+    if bounding_box:
+        wcsobj.bounding_box = bounding_box
 
     return wcsobj
diff --git a/romancal/pipeline/tests/test_mosaic_pipeline.py b/romancal/pipeline/tests/test_mosaic_pipeline.py
new file mode 100644
index 000000000..79bcd2d9b
--- /dev/null
+++ b/romancal/pipeline/tests/test_mosaic_pipeline.py
@@ -0,0 +1,133 @@
+"""Unit tests for the mosaic pipeline"""
+
+import numpy as np
+
+import romancal.pipeline.mosaic_pipeline as mp
+
+
+def test_skycell_to_wcs():
+    """Test integrity of skycell_to_wcs"""
+
+    skycell = np.void(
+        (
+            "r274dp63x31y81",
+            269.7783307416819,
+            66.04965143695566,
+            1781.5,
+            1781.5,
+            355.9788,
+            3564,
+            3564,
+            67715.5,
+            -110484.5,
+            269.6657957545588,
+            65.9968687812357,
+            269.6483032937494,
+            66.09523979539262,
+            269.89132874168854,
+            66.10234971630734,
+            269.9079118635897,
+            66.00394719483091,
+            0.1,
+            274.2857142857143,
+            63.0,
+            0.0,
+            463181,
+        ),
+        dtype=[
+            ("name", "<U20"),
+            ("ra_center", "<f8"),
+            ("dec_center", "<f8"),
+            ("x_center", "<f4"),
+            ("y_center", "<f4"),
+            ("orientat", "<f4"),
+            ("nx", "<i4"),
+            ("ny", "<i4"),
+            ("x0_projection", "<f4"),
+            ("y0_projection", "<f4"),
+            ("ra_corn1", "<f8"),
+            ("dec_corn1", "<f8"),
+            ("ra_corn2", "<f8"),
+            ("dec_corn2", "<f8"),
+            ("ra_corn3", "<f8"),
+            ("dec_corn3", "<f8"),
+            ("ra_corn4", "<f8"),
+            ("dec_corn4", "<f8"),
+            ("pixel_scale", "<f4"),
+            ("ra_projection_center", "<f8"),
+            ("dec_projection_center", "<f8"),
+            ("orientat_projection_center", "<f4"),
+            ("index", "<i8"),
+        ],
+    )
+
+    wcs = mp.skycell_to_wcs(skycell)
+
+    assert np.allclose(
+        wcs(
+            skycell["x0_projection"], skycell["y0_projection"], with_bounding_box=False
+        ),
+        (skycell["ra_projection_center"], skycell["dec_projection_center"]),
+    )
+    assert np.allclose(
+        wcs(skycell["x_center"], skycell["y_center"]),
+        (skycell["ra_center"], skycell["dec_center"]),
+    )
+    assert np.allclose(wcs(0.0, 0.0), (skycell["ra_corn1"], skycell["dec_corn1"]))
+    assert np.allclose(
+        wcs(0.0, skycell["ny"] - 1), (skycell["ra_corn2"], skycell["dec_corn2"])
+    )
+    assert np.allclose(
+        wcs(skycell["nx"] - 1, skycell["ny"] - 1),
+        (skycell["ra_corn3"], skycell["dec_corn3"]),
+    )
+    assert np.allclose(
+        wcs(skycell["nx"] - 1, 0.0), (skycell["ra_corn4"], skycell["dec_corn4"])
+    )
+
+
+def test_wcsinfo_to_wcs():
+    """Test integrity of wcsinfo_to_wcs"""
+    wcsinfo = {
+        "ra_ref": 269.83219987378925,
+        "dec_ref": 66.04081466149024,
+        "x_ref": 2069.0914958388985,
+        "y_ref": 2194.658767532754,
+        "rotation_matrix": [
+            [-0.9999964196507396, -0.00267594575838714],
+            [-0.00267594575838714, 0.9999964196507396],
+        ],
+        "pixel_scale": 3.036307317109957e-05,
+        "pixel_shape": [4389, 4138],
+        "ra_center": 269.82284964811464,
+        "dec_center": 66.0369888162117,
+        "ra_corn1": 269.98694025887136,
+        "dec_corn1": 65.97426875366378,
+        "ra_corn2": 269.98687579251805,
+        "dec_corn2": 66.09988065827382,
+        "ra_corn3": 269.6579498847431,
+        "dec_corn3": 66.099533603104,
+        "ra_corn4": 269.6596332616879,
+        "dec_corn4": 65.97389321243348,
+        "orientat": 359.8466793994546,
+    }
+
+    wcs = mp.wcsinfo_to_wcs(wcsinfo)
+
+    assert np.allclose(
+        wcs(wcsinfo["x_ref"], wcsinfo["y_ref"]), (wcsinfo["ra_ref"], wcsinfo["dec_ref"])
+    )
+    assert np.allclose(
+        wcs(4389 / 2.0, 4138 / 2.0), (wcsinfo["ra_center"], wcsinfo["dec_center"])
+    )
+    assert np.allclose(wcs(0.0, 0.0), (wcsinfo["ra_corn1"], wcsinfo["dec_corn1"]))
+    assert np.allclose(
+        wcs(0.0, wcsinfo["pixel_shape"][1]), (wcsinfo["ra_corn2"], wcsinfo["dec_corn2"])
+    )
+    assert np.allclose(
+        wcs(wcsinfo["pixel_shape"][0], wcsinfo["pixel_shape"][1]),
+        (wcsinfo["ra_corn3"], wcsinfo["dec_corn3"]),
+    )
+    assert np.allclose(
+        wcs(wcsinfo["pixel_shape"][0], 0.0), (wcsinfo["ra_corn4"], wcsinfo["dec_corn4"])
+    )
diff --git a/romancal/regtest/test_mos_pipeline.py b/romancal/regtest/test_mos_pipeline.py
index 8fba060af..dfa2247b7 100644
--- a/romancal/regtest/test_mos_pipeline.py
+++ b/romancal/regtest/test_mos_pipeline.py
@@ -5,8 +5,10 @@
 import pytest
 import roman_datamodels as rdm
 
+from romancal.pipeline import mosaic_pipeline
 from romancal.pipeline.mosaic_pipeline import MosaicPipeline
 
+from . import util
 from .regtestdata import compare_asdf
 
 # mark all tests in this module
@@ -116,3 +118,11 @@ def test_added_background(output_model):
 def test_added_background_level(output_model):
     # DMS400
     assert any(output_model.meta.individual_image_meta.background["level"] != 0)
+
+
+def test_wcsinfo_wcs_roundtrip(output_model):
+    """Test that the contents of wcsinfo reproduces the wcs"""
+    wcs_from_wcsinfo = mosaic_pipeline.wcsinfo_to_wcs(output_model.meta.wcsinfo)
+
+    ra_mad, dec_mad = util.comp_wcs_grids_arcs(output_model.meta.wcs, wcs_from_wcsinfo)
+    assert (ra_mad + dec_mad) / 2.0 < 1.0e-5
diff --git a/romancal/regtest/test_mos_skycell_pipeline.py b/romancal/regtest/test_mos_skycell_pipeline.py
index ce5953d96..eddbc3a31 100644
--- a/romancal/regtest/test_mos_skycell_pipeline.py
+++ b/romancal/regtest/test_mos_skycell_pipeline.py
@@ -1,8 +1,10 @@
 import pytest
 import roman_datamodels as rdm
 
+from romancal.pipeline import mosaic_pipeline
 from romancal.pipeline.mosaic_pipeline import MosaicPipeline
 
+from . import util
 from .regtestdata import compare_asdf
 
 # mark all tests in this module
@@ -56,3 +58,11 @@ def test_resample_ran(output_model):
 def test_location_name(output_model):
     # test that the location_name matches the skycell selected
     assert output_model.meta.basic.location_name == "r274dp63x31y81"
+
+
+def test_wcsinfo_wcs_roundtrip(output_model):
+    """Test that the contents of wcsinfo reproduces the wcs"""
+    wcs_from_wcsinfo = mosaic_pipeline.wcsinfo_to_wcs(output_model.meta.wcsinfo)
+
+    ra_mad, dec_mad = util.comp_wcs_grids_arcs(output_model.meta.wcs, wcs_from_wcsinfo)
+    assert (ra_mad + dec_mad) / 2.0 < 1.0e-5
diff --git a/romancal/regtest/util.py b/romancal/regtest/util.py
new file mode 100644
index 000000000..31bf2605f
--- /dev/null
+++ b/romancal/regtest/util.py
@@ -0,0 +1,33 @@
+"""Test utilities"""
+
+import numpy as np
+from astropy.stats import mad_std
+
+
+def comp_wcs_grids_arcs(wcs_a, wcs_b, npix=4088, interval=10):
+    """Compare world grids produced by the two wcs
+
+    Parameters
+    ----------
+    wcs_a, wcs_b : gwcs.WCS
+        The wcs object to compare.
+
+    npix : int
+        The size of the grid to produce.
+
+    interval : int
+        The interval to check over.
+
+    Returns
+    -------
+    mad_std : float
+        The numpy MAD_STD in arcseconds
+    """
+    xx, yy = np.meshgrid(np.linspace(0, npix, interval), np.linspace(0, npix, interval))
+    ra_a, dec_a = wcs_a(xx, yy, with_bounding_box=False)
+    ra_b, dec_b = wcs_b(xx, yy, with_bounding_box=False)
+
+    ra_mad = mad_std(ra_a - ra_b, ignore_nan=True) * 60.0 * 60.0 * 1000.0
+    dec_mad = mad_std(dec_a - dec_b, ignore_nan=True) * 60.0 * 60.0 * 1000.0
+
+    return ra_mad, dec_mad
diff --git a/romancal/resample/resample.py b/romancal/resample/resample.py
index e320c6104..c359f1593 100644
--- a/romancal/resample/resample.py
+++ b/romancal/resample/resample.py
@@ -787,13 +787,13 @@ def gwcs_into_l3(model, wcs):
     l3_wcsinfo.projection = "TAN"
     l3_wcsinfo.pixel_shape = model.shape
 
+    # Fill out image-local information
     pixel_center = [(v - 1) / 2.0 for v in model.shape[::-1]]
     world_center = wcs(*pixel_center)
     l3_wcsinfo.ra_center = world_center[0]
     l3_wcsinfo.dec_center = world_center[1]
     l3_wcsinfo.pixel_scale_local = compute_scale(wcs, world_center)
     l3_wcsinfo.orientat_local = calc_pa(wcs, *world_center)
-
     try:
         footprint = utils.create_footprint(wcs, model.shape)
     except Exception as excp:
@@ -809,6 +809,7 @@ def gwcs_into_l3(model, wcs):
         l3_wcsinfo.dec_corn4 = footprint[3][1]
         l3_wcsinfo.s_region = compute_s_region_keyword(footprint)
 
+    # Fill out wcs-general information
     try:
         l3_wcsinfo.x_ref = -transform["crpix1"].offset.value
         l3_wcsinfo.y_ref = -transform["crpix2"].offset.value
@@ -818,10 +819,18 @@ def gwcs_into_l3(model, wcs):
         )
         l3_wcsinfo.x_ref = pixel_center[0]
         l3_wcsinfo.y_ref = pixel_center[1]
-    world_ref = wcs(l3_wcsinfo.x_ref, l3_wcsinfo.y_ref)
+
+    world_ref = wcs(l3_wcsinfo.x_ref, l3_wcsinfo.y_ref, with_bounding_box=False)
     l3_wcsinfo.ra_ref = world_ref[0]
     l3_wcsinfo.dec_ref = world_ref[1]
-    l3_wcsinfo.pixel_scale = compute_scale(wcs, world_ref)
+
+    try:
+        cdelt1 = transform["cdelt1"].factor.value
+        cdelt2 = transform["cdelt2"].factor.value
+        l3_wcsinfo.pixel_scale = (cdelt1 + cdelt2) / 2.0
+    except IndexError:
+        l3_wcsinfo.pixel_scale = compute_scale(wcs, world_ref)
+
     l3_wcsinfo.orientat = calc_pa(wcs, *world_ref)
 
     try:
@@ -853,9 +862,11 @@ def calc_pa(wcs, ra, dec):
         The position angle in degrees.
 
     """
-    delta_pix = [v for v in wcs.world_to_pixel_values(ra, dec)]
+    delta_pix = [v for v in wcs.invert(ra, dec, with_bounding_box=False)]
     delta_pix[1] += 1
-    delta_coord = wcs.pixel_to_world(*delta_pix)
+    delta_coord = SkyCoord(
+        *wcs(*delta_pix, with_bounding_box=False), frame="icrs", unit="deg"
+    )
     coord = SkyCoord(ra, dec, frame="icrs", unit="deg")
 
     return coord.position_angle(delta_coord).degree
diff --git a/romancal/resample/tests/test_resample.py b/romancal/resample/tests/test_resample.py
index cd6b603b8..bd782c607 100644
--- a/romancal/resample/tests/test_resample.py
+++ b/romancal/resample/tests/test_resample.py
@@ -1076,8 +1076,8 @@ def test_l3_wcsinfo(multiple_exposures):
             "dec_corn3": 0.005043059486913547,
             "ra_corn4": 9.998944914237953,
             "dec_corn4": 0.00038355958555111314,
-            "orientat_local": 9.826983262839223,
-            "orientat": 9.826978421513601,
+            "orientat_local": 20.999999978134802,
+            "orientat": 20.999999978134802,
             "projection": "TAN",
             "s_region": (
                 "POLYGON ICRS 10.005109345783163 -0.001982743978690467 10.006897960220385 "