diff --git a/analysis/power_spectrum.py b/analysis/power_spectrum.py
index 93bb8ee..cb8ab73 100644
--- a/analysis/power_spectrum.py
+++ b/analysis/power_spectrum.py
@@ -5,6 +5,24 @@ kmax = 2e0
 Nk = 50
 AliasingCorr=False
 
+def crop_field(field, Ncrop):
+
+    if Ncrop is None or Ncrop == 0:
+        return
+    
+    elif Ncrop > 0:
+        field.data = field.data[Ncrop:-Ncrop, Ncrop:-Ncrop, Ncrop:-Ncrop]
+        d0 = field.L0/field.N0
+        d1 = field.L1/field.N1
+        d2 = field.L2/field.N2
+        field.N0 -= 2*Ncrop
+        field.N1 -= 2*Ncrop
+        field.N2 -= 2*Ncrop
+        field.L0 = field.N0*d0
+        field.L1 = field.N1*d1
+        field.L2 = field.N2*d2
+
+
 def get_power_spectrum(field, kmin=kmin, kmax=kmax, Nk=Nk, G=None):
     from pysbmy.power import PowerSpectrum
     from pysbmy.fft import FourierGrid
@@ -91,7 +109,8 @@ def plot_power_spectra(filenames,
                        figsize=(8,4),
                        dpi=300,
                        ax=None,
-                       fig=None,):
+                       fig=None,
+                       Ncrop=None,):
     
     import matplotlib.pyplot as plt
     from pysbmy.field import read_field 
@@ -110,6 +129,7 @@ def plot_power_spectra(filenames,
     
     for i, filename in enumerate(filenames):
         field = read_field(filename)
+        crop_field(field, Ncrop)
         _, G, k, _ = add_power_spectrum_to_plot(ax=ax,
                                    field=field,
                                    Pk_ref=Pk_ref,
@@ -163,7 +183,9 @@ def plot_cross_correlations(filenames_A,
                             figsize=(8,4),
                             dpi=300,
                             ax=None,
-                            fig=None,):
+                            fig=None,
+                            Ncrop=None,
+                            ):
     
     import matplotlib.pyplot as plt
     from pysbmy.field import read_field 
@@ -181,9 +203,11 @@ def plot_cross_correlations(filenames_A,
         markers = [None for f in filenames_A]
 
     field_B = read_field(filename_B)
+    crop_field(field_B, Ncrop)
     
     for i, filename_A in enumerate(filenames_A):
         field_A = read_field(filename_A)
+        crop_field(field_A, Ncrop)
         _, G, k, _ = add_cross_correlations_to_plot(ax=ax,
                                        field_A=field_A,
                                        field_B=field_B,
@@ -255,6 +279,7 @@ if __name__ == "__main__":
     parser.add_argument('-t','--title', type=str, default=None, help='Title of the plot.')
     parser.add_argument('-yrp', '--ylim_power', type=float, nargs=2, default=[0.9,1.1], help='Y-axis limits.')
     parser.add_argument('-yrc', '--ylim_corr', type=float, nargs=2, default=[0.99,1.001], help='Y-axis limits.')
+    parser.add_argument('--crop', type=int, default=None, help='Remove the outter N pixels of the fields.')
 
     args = parser.parse_args()
 
@@ -264,7 +289,9 @@ if __name__ == "__main__":
 
     if args.reference is not None:
         from pysbmy.field import read_field
-        G, _, Pk_ref = get_power_spectrum(read_field(args.directory+args.reference), kmin=kmin, kmax=kmax, Nk=Nk)
+        F_ref = read_field(args.directory+args.reference)
+        crop_field(F_ref, args.crop)
+        G, _, Pk_ref = get_power_spectrum(F_ref, kmin=kmin, kmax=kmax, Nk=Nk)
     else:
         Pk_ref = None
         G = None
@@ -294,7 +321,9 @@ if __name__ == "__main__":
                            kmax=kmax,
                            Nk=Nk,
                            ax=axes[0],
-                           fig=fig)
+                           fig=fig,
+                           Ncrop=args.crop,
+                           )
         
         plot_cross_correlations(filenames_A=filenames,
                                 filename_B=args.directory+args.reference,
@@ -311,7 +340,9 @@ if __name__ == "__main__":
                                 kmax=kmax,
                                 Nk=Nk,
                                 ax=axes[1],
-                                fig=fig)
+                                fig=fig,
+                                Ncrop=args.crop,
+                                )
         
         axes[1].legend(loc='lower left')
         axes[0].set_title("Power Spectrum")
@@ -334,7 +365,9 @@ if __name__ == "__main__":
                                      bound2=0.02,
                                      kmin=kmin,
                                      kmax=kmax,
-                                     Nk=Nk)
+                                     Nk=Nk,
+                                    Ncrop=args.crop,
+                                    )
         ax.legend()
         if args.title is not None:
             ax.set_title(args.title)
@@ -353,7 +386,9 @@ if __name__ == "__main__":
                                           bound2=0.002,
                                           kmin=kmin,
                                           kmax=kmax,
-                                          Nk=Nk)
+                                          Nk=Nk,
+                                        Ncrop=args.crop,
+                                        )
         ax.legend(loc='lower left')
         if args.title is not None:
             ax.set_title(args.title)