diff --git a/btk/catalog.py b/btk/catalog.py
index 3a1071bdf..33b50e663 100644
--- a/btk/catalog.py
+++ b/btk/catalog.py
@@ -29,6 +29,7 @@ def __init__(self, raw_catalog, verbose=False):
         """
         self.verbose = verbose
         self.table = self._prepare_table(raw_catalog)
+        self._raw_catalog = raw_catalog
 
         if self.verbose:
             print("Catalog loaded")
diff --git a/btk/draw_blends.py b/btk/draw_blends.py
index 1b698219d..280b16dac 100644
--- a/btk/draw_blends.py
+++ b/btk/draw_blends.py
@@ -76,24 +76,22 @@ def get_catsim_galaxy(entry, filt, survey, no_disk=False, no_bulge=False, no_agn
         raise SourceNotVisible
 
     if disk_flux > 0:
-        beta_radians = np.radians(entry["pa_disk"])
         if bulge_flux > 0:
             assert entry["pa_disk"] == entry["pa_bulge"], "Sersic components have different beta."
         a_d, b_d = entry["a_d"], entry["b_d"]
         disk_hlr_arcsecs = np.sqrt(a_d * b_d)
         disk_q = b_d / a_d
         disk = galsim.Exponential(flux=disk_flux, half_light_radius=disk_hlr_arcsecs).shear(
-            q=disk_q, beta=beta_radians * galsim.radians
+            q=disk_q, beta=entry["pa_disk"] * galsim.degrees
         )
         components.append(disk)
 
     if bulge_flux > 0:
-        beta_radians = np.radians(entry["pa_bulge"])
         a_b, b_b = entry["a_b"], entry["b_b"]
         bulge_hlr_arcsecs = np.sqrt(a_b * b_b)
         bulge_q = b_b / a_b
         bulge = galsim.DeVaucouleurs(flux=bulge_flux, half_light_radius=bulge_hlr_arcsecs).shear(
-            q=bulge_q, beta=beta_radians * galsim.radians
+            q=bulge_q, beta=entry["pa_bulge"] * galsim.degrees
         )
         components.append(bulge)
 
@@ -112,7 +110,6 @@ class DrawBlendsGenerator(ABC):
     each mini-batch analyzed separately. The results are then combined to output a
     dict with results of entire batch. If the number of cpus is greater than one, then each of
     the mini-batches are run in parallel.
-
     """
 
     compatible_catalogs = ("Catalog",)
@@ -167,8 +164,9 @@ def __init__(
             apply_shear (float): Whether to apply the shear specified in catalogs to galaxies.
                                 If set to True, sampling function must add 'g1', 'g2' columns.
             augment_data (float): If set to True, augment data by adding a random rotation to every
-                                galaxy drawn. Rotation added is recorded via `btk_rotation` column
-                                output.
+                                galaxy drawn. Rotation added is proapaged to the `pa_bulge`
+                                and `pa_disk` columns if using the `CatsimGenerator`. It is also
+                                stored in the `btk_rotation` column.
         """
         self.blend_generator = BlendGenerator(
             catalog, sampling_function, batch_size, shifts, indexes, verbose
@@ -428,7 +426,6 @@ def render_blend(self, blend_catalog, psf, filt, survey, seedseq_blend, extra_da
 
         Returns:
             Images of blend and isolated galaxies as `numpy.ndarray`.
-
         """
         sky_level = mean_sky_level(survey, filt).to_value("electron")
         blend_catalog.add_column(
@@ -517,9 +514,10 @@ def render_single(self, entry, filt, psf, survey, extra_data):
 
         pix_stamp_size = int(self.stamp_size / survey.pixel_scale.to_value("arcsec"))
         try:
-            gal = get_catsim_galaxy(entry, filt, survey)
             if self.augment_data:
-                gal.rotate(galsim.Angle(entry["btk_rotation"], unit=galsim.degrees))
+                entry["pa_bulge"] = (entry["pa_bulge"] + entry["btk_rotation"]) % 360
+                entry["pa_disk"] = (entry["pa_disk"] + entry["btk_rotation"]) % 360
+            gal = get_catsim_galaxy(entry, filt, survey)
             if self.apply_shear:
                 if "g1" in entry.keys() and "g2" in entry.keys():
                     gal = gal.shear(g1=entry["g1"], g2=entry["g2"])