diff --git a/CHANGES.rst b/CHANGES.rst
index 15ba630c4..c8053471c 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -13,9 +13,13 @@ Bug Fixes
 - ``with_spectral_unit`` has been changed to ``with_spectral_axis_unit`` and actually works
   now. [#1119]
 
+- Fixed numpy error when printing a ``Spectrum1D`` object. [#1123]
+
 Other Changes and Additions
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
+- Updated the format of ``Spectrum1D.__str__`` and ``Spectrum1D.__repr__``. [#1123]
+
 1.12.0 (2023-10-17)
 -------------------
 
diff --git a/docs/arithmetic.rst b/docs/arithmetic.rst
index 663dc3317..aa8b1fbf3 100644
--- a/docs/arithmetic.rst
+++ b/docs/arithmetic.rst
@@ -29,23 +29,11 @@ Arithmetic support includes addition, subtract, multiplication, and division.
     >>> spec2 = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, flux=rng.random(49)*u.Jy)
     >>> spec3 = spec1 + spec2
     >>> spec3  # doctest: +FLOAT_CMP
-    <Spectrum1D(flux=<Quantity [0.85594057, 0.59914105, 0.86545315, 1.29308365, 0.56743587,
-               0.63720232, 1.03919556, 0.33693653, 0.89068491, 1.41613598,
-               0.72461457, 1.20411351, 0.96480272, 0.37496506, 0.70241888,
-               1.36924151, 0.90943254, 0.82210346, 0.98018949, 1.05861761,
-               0.26172516, 1.02205189, 0.51839963, 1.21572449, 0.87754143,
-               1.02493144, 0.95316681, 0.37872965, 0.17723648, 1.21662474,
-               1.39171024, 1.23614795, 1.10636247, 0.97294585, 1.5453743 ,
-               1.01020945, 1.42125961, 1.36636734, 1.11338214, 0.58687869,
-               0.63074156, 0.67475105, 0.54093389, 1.77345469, 1.22990398,
-               0.78162068, 0.63760289, 0.63139356, 0.97112644] Jy>, spectral_axis=<SpectralAxis
+    <Spectrum1D(flux=[0.8559405665668484 ... 0.9711264429515736] Jy (shape=(49,), mean=0.91592 Jy); spectral_axis=<SpectralAxis
        (observer to target:
           radial_velocity=0.0 km / s
           redshift=0.0)
-      [ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12., 13., 14.,
-       15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28.,
-       29., 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., 40., 41., 42.,
-       43., 44., 45., 46., 47., 48., 49.] nm>)>
+      [ 1.  2.  3. ... 47. 48. 49.] nm> (length=49))>
 
 
 Propagation of Uncertainties
diff --git a/docs/manipulation.rst b/docs/manipulation.rst
index 04e45e0e1..f82bac80a 100644
--- a/docs/manipulation.rst
+++ b/docs/manipulation.rst
@@ -59,19 +59,7 @@ along the spectral dimension.
     >>> spec1_gsmooth = gaussian_smooth(spec1, stddev=3)
     >>> spec1_tsmooth = trapezoid_smooth(spec1, width=3)
     >>> gaussian_smooth(spec1, stddev=3) # doctest: +FLOAT_CMP
-    <Spectrum1D(flux=<Quantity [0.25860917, 0.32640733, 0.38562205, 0.43159269, 0.46445898,
-               0.48828727, 0.50822537, 0.52741247, 0.54565787, 0.56035736,
-               0.56880835, 0.57015658, 0.56561802, 0.55700108, 0.54472603,
-               0.52720679, 0.50217622, 0.46914624, 0.43114785, 0.39432907,
-               0.36548049, 0.34885759, 0.34456742, 0.34936895, 0.35949514,
-               0.37376516, 0.39451092, 0.42591854, 0.47068509, 0.5272207 ,
-               0.58920337, 0.64763434, 0.69379067, 0.72157276, 0.7284874 ,
-               0.71560131, 0.68696137, 0.64870578, 0.60757162, 0.56898454,
-               0.53520565, 0.5047466 , 0.47338241, 0.43641424, 0.39090495,
-               0.3368543 , 0.27706956, 0.21605198, 0.15868783] Jy>, spectral_axis=<SpectralAxis [ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12., 13., 14.,
-       15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28.,
-       29., 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., 40., 41., 42.,
-       43., 44., 45., 46., 47., 48., 49.] nm>)>
+    <Spectrum1D(flux=[0.25860917267578276 ... 0.15868783272966752] Jy (shape=(49,), mean=0.48009 Jy); spectral_axis=<SpectralAxis [ 1.  2.  3. ... 47. 48. 49.] nm> (length=49))>
 
 Each of the specific smoothing methods create the appropriate `astropy.convolution.convolve`
 kernel and then call a helper function :func:`~specutils.manipulation.convolution_smooth`
@@ -87,19 +75,7 @@ that takes the spectrum and an astropy 1D kernel.  So, one could also do:
     >>> spec1 = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm,
     ...                    flux=np.random.default_rng(12345).random(49) * u.Jy)
     >>> convolution_smooth(spec1, box1d_kernel) # doctest: +FLOAT_CMP
-    <Spectrum1D(flux=<Quantity [0.18136479, 0.44715327, 0.59679282, 0.62157656, 0.46672605,
-               0.44074408, 0.37261896, 0.48593299, 0.60043103, 0.62093487,
-               0.71297658, 0.62145477, 0.57067218, 0.40165835, 0.47473917,
-               0.6752577 , 0.63680209, 0.58595151, 0.42684533, 0.34521923,
-               0.15387504, 0.19386345, 0.32172965, 0.35723812, 0.51579686,
-               0.42516394, 0.37951329, 0.13814274, 0.27323395, 0.51499156,
-               0.68497705, 0.79611717, 0.75279699, 0.83910701, 0.83822349,
-               0.77869171, 0.80439892, 0.65961564, 0.56881136, 0.40684661,
-               0.46353027, 0.48256952, 0.63312795, 0.4971397 , 0.50011884,
-               0.28525197, 0.31804274, 0.20644074, 0.12015627] Jy>, spectral_axis=<SpectralAxis [ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12., 13., 14.,
-       15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28.,
-       29., 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., 40., 41., 42.,
-       43., 44., 45., 46., 47., 48., 49.] nm>)>
+    <Spectrum1D(flux=[0.1813647873923075 ... 0.1201562712204726] Jy (shape=(49,), mean=0.49378 Jy); spectral_axis=<SpectralAxis [ 1.  2.  3. ... 47. 48. 49.] nm> (length=49))>
 
 In this case, the ``spec1_bsmooth2`` result should be equivalent to the ``spec1_bsmooth`` in
 the section above (assuming the flux data of the input ``spec`` is the same). Note that,
@@ -131,19 +107,7 @@ method applys the median filter across the flux.
     >>> spec1 = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm,
     ...                    flux=np.random.default_rng(12345).random(49) * u.Jy)
     >>> median_smooth(spec1, width=3) # doctest: +FLOAT_CMP
-    <Spectrum1D(flux=<Quantity [0.22733602, 0.31675834, 0.67625467, 0.67625467, 0.39110955,
-               0.39110955, 0.33281393, 0.59830875, 0.67275604, 0.67275604,
-               0.94180287, 0.66723745, 0.66723745, 0.44183967, 0.44183967,
-               0.6974535 , 0.6974535 , 0.6974535 , 0.32647286, 0.22013496,
-               0.1598956 , 0.1598956 , 0.34010018, 0.34010018, 0.46519315,
-               0.26642103, 0.19329439, 0.12946908, 0.12946908, 0.59856801,
-               0.60162124, 0.8547419 , 0.72478136, 0.86055132, 0.86055132,
-               0.86055132, 0.9293378 , 0.54618601, 0.49498794, 0.45177871,
-               0.45177871, 0.45177871, 0.66503892, 0.33089093, 0.33982834,
-               0.2588534 , 0.33982834, 0.2588534 , 0.00502233] Jy>, spectral_axis=<SpectralAxis [ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12., 13., 14.,
-       15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28.,
-       29., 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., 40., 41., 42.,
-       43., 44., 45., 46., 47., 48., 49.] nm>)>
+    <Spectrum1D(flux=[0.22733602246716966 ... 0.005022333717131788] Jy (shape=(49,), mean=0.48620 Jy); spectral_axis=<SpectralAxis [ 1.  2.  3. ... 47. 48. 49.] nm> (length=49))>
 
 Resampling
 ----------
@@ -374,14 +338,14 @@ the ``spectral_axis``. Therefore one can use a construct like this:
     >>> spectrum  # doctest: +FLOAT_CMP
     <Spectrum1D(flux=<Quantity [ 30.47170798, 103.99841062,  75.04511958,  94.05647164,
                195.10351887, 130.21795069,  12.78404032,  31.62425923,
-                 1.68011575,  85.30439276] Jy>, spectral_axis=<SpectralAxis [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.] um>)>
+                 1.68011575,  85.30439276] Jy> (shape=(10,), mean=76.02860 Jy); spectral_axis=<SpectralAxis [0. 1. 2. ... 7. 8. 9.] um> (length=10))>
 
     >>> shift = 12300 * u.AA
     >>> new_spec = Spectrum1D(spectral_axis=spectrum.spectral_axis + shift, flux=spectrum.flux)
     >>> new_spec  # doctest: +FLOAT_CMP
     <Spectrum1D(flux=<Quantity [ 30.47170798, 103.99841062,  75.04511958,  94.05647164,
                195.10351887, 130.21795069,  12.78404032,  31.62425923,
-                 1.68011575,  85.30439276] Jy>, spectral_axis=<SpectralAxis [ 1.23,  2.23,  3.23,  4.23,  5.23,  6.23,  7.23,  8.23,  9.23, 10.23] um>)>
+                 1.68011575,  85.30439276] Jy> (shape=(10,), mean=76.02860 Jy); spectral_axis=<SpectralAxis [ 1.23  2.23  3.23 ...  8.23  9.23 10.23] um> (length=10))>
 
 Replacing a region
 ------------------
@@ -403,7 +367,7 @@ with the spline knots:
     >>> spline_knots = [3.5, 4.7, 6.8, 7.1] * u.AA
     >>> result = model_replace(input_spectrum, None, model=spline_knots)
     >>> result
-    <Spectrum1D(flux=<Quantity [ 2.,  4.,  6.,  8., 10., 12., 14., 16., 18., 20.] mJy>, spectral_axis=<SpectralAxis [ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.] Angstrom>)>
+    <Spectrum1D(flux=<Quantity [ 2.,  4.,  6.,  8., 10., 12., 14., 16., 18., 20.] mJy> (shape=(10,), mean=11.00000 mJy); spectral_axis=<SpectralAxis [ 1.  2.  3. ...  8.  9. 10.] Angstrom> (length=10))>
 
 The default behavior is to keep the data outside the replaced region unchanged.
 Alternatively, the spectrum outside the replaced region can be filled with zeros:
@@ -413,8 +377,7 @@ Alternatively, the spectrum outside the replaced region can be filled with zeros
     >>> spline_knots = [3.5, 4.7, 6.8, 7.1] * u.AA
     >>> result = model_replace(input_spectrum, None, model=spline_knots, extrapolation_treatment='zero_fill')
     >>> result
-    <Spectrum1D(flux=<Quantity [ 0.,  0.,  0.,  8., 10., 12., 14.,  0.,  0.,  0.] mJy>,
-        spectral_axis=<SpectralAxis [ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.] Angstrom>)>
+    <Spectrum1D(flux=<Quantity [ 0.,  0.,  0.,  8., 10., 12., 14.,  0.,  0.,  0.] mJy> (shape=(10,), mean=4.40000 mJy); spectral_axis=<SpectralAxis [ 1.  2.  3. ...  8.  9. 10.] Angstrom> (length=10))>
 
 One can define the spline knots by providing an instance of `~specutils.SpectralRegion`,
 and the number of knots to be evenly spread along the region:
@@ -425,8 +388,7 @@ and the number of knots to be evenly spread along the region:
     >>> region = SpectralRegion(3.5*u.AA, 7.1*u.AA)
     >>> result = model_replace(input_spectrum, region, model=4)
     >>> result
-    <Spectrum1D(flux=<Quantity [ 2.,  4.,  6.,  8., 10., 12., 14., 16., 18., 20.] mJy>,
-        spectral_axis=<SpectralAxis [ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.] Angstrom>)>
+    <Spectrum1D(flux=<Quantity [ 2.,  4.,  6.,  8., 10., 12., 14., 16., 18., 20.] mJy> (shape=(10,), mean=11.00000 mJy); spectral_axis=<SpectralAxis [ 1.  2.  3. ...  8.  9. 10.] Angstrom> (length=10))>
 
 A model fitted over the region can also be used to replace the spectrum flux values:
 
@@ -442,7 +404,7 @@ A model fitted over the region can also be used to replace the spectrum flux val
     >>> result = model_replace(input_spectrum, region, model=fitted_model)
     >>> result  # doctest: +FLOAT_CMP
     <Spectrum1D(flux=<Quantity [1.        , 1.1       , 0.9       , 4.40801804, 9.58271877,
-               5.61238054, 0.88556096, 1.        , 1.2       , 1.1       ] mJy>, spectral_axis=<SpectralAxis [ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.] Angstrom>)>
+               5.61238054, 0.88556096, 1.        , 1.2       , 1.1       ] mJy> (shape=(10,), mean=2.67887 mJy); spectral_axis=<SpectralAxis [ 1.  2.  3. ...  8.  9. 10.] Angstrom> (length=10))>
 
 Reference/API
 -------------
diff --git a/docs/spectral_cube.rst b/docs/spectral_cube.rst
index 42fd2bd26..8a8d4bd8b 100644
--- a/docs/spectral_cube.rst
+++ b/docs/spectral_cube.rst
@@ -49,31 +49,11 @@ Print the contents of 3 spectral axis points in a 3x3 spaxel array:
 .. code-block:: python
 
     >>> sc[30:33,30:33,2000:2003] # doctest: +REMOTE_DATA
-    <Spectrum1D(flux=<Quantity [[[0.48920232, 0.4987253 , 0.5098349 ],
-                [0.493365  , 0.4964812 , 0.5223962 ],
-                [0.49446177, 0.4909543 , 0.5304416 ]],
-    <BLANKLINE>
-               [[0.53114057, 0.53538376, 0.5467784 ],
-                [0.53761804, 0.533159  , 0.554437  ],
-                [0.5470889 , 0.54905874, 0.57109433]],
-    <BLANKLINE>
-               [[0.5599331 , 0.554316  , 0.5618426 ],
-                [0.5763055 , 0.5668046 , 0.5774939 ],
-                [0.59571505, 0.60118765, 0.59942234]]] 1e-17 erg / (Angstrom s spaxel cm2)>, spectral_axis=<SpectralAxis
+    <Spectrum1D(flux=[[[0.4892023205757141 ... 0.5994223356246948]]] 1e-17 erg / (Angstrom s spaxel cm2) (shape=(3, 3, 3), mean=0.54165 1e-17 erg / (Angstrom s spaxel cm2)); spectral_axis=<SpectralAxis
        (observer to target:
           radial_velocity=0.0 km / s
           redshift=0.0)
-      [5.73984286e-07, 5.74116466e-07, 5.74248676e-07] m>, uncertainty=InverseVariance([[[4324.235 , 4326.87  , 4268.985 ],
-                      [5128.3867, 5142.5005, 4998.457 ],
-                      [4529.9097, 4545.8345, 4255.305 ]],
-    <BLANKLINE>
-                     [[4786.163 , 4811.216 , 4735.3135],
-                      [4992.71  , 5082.1294, 4927.881 ],
-                      [4992.9683, 5046.971 , 4798.005 ]],
-    <BLANKLINE>
-                     [[4831.2236, 4887.096 , 4806.84  ],
-                      [3895.8677, 4027.9104, 3896.0195],
-                      [4521.258 , 4630.997 , 4503.0396]]]))>
+      [5.73984286e-07 5.74116466e-07 5.74248676e-07] m> (length=3); uncertainty=InverseVariance)>
 
 
 Spectral slab extraction
@@ -90,31 +70,11 @@ spectral regions from the cube.
     >>> ss.shape  # doctest: +REMOTE_DATA
     (74, 74, 3)
     >>> ss[30:33,30:33,::] # doctest: +REMOTE_DATA
-    <Spectrum1D(flux=<Quantity [[[0.6103081 , 0.95697385, 1.0791174 ],
-                [0.5663384 , 0.8872061 , 1.0814004 ],
-                [0.520966  , 0.7819859 , 1.024845  ]],
-    <BLANKLINE>
-               [[0.64514536, 0.96376216, 1.083235  ],
-                [0.6112465 , 0.89025146, 1.058679  ],
-                [0.56316894, 0.77895504, 0.99165994]],
-    <BLANKLINE>
-               [[0.65954393, 0.9084677 , 0.9965009 ],
-                [0.6255246 , 0.84401435, 0.9930112 ],
-                [0.59066033, 0.762025  , 0.9361185 ]]] 1e-17 erg / (Angstrom s spaxel cm2)>, spectral_axis=<SpectralAxis
+    <Spectrum1D(flux=[[[0.6103081107139587 ... 0.936118483543396]]] 1e-17 erg / (Angstrom s spaxel cm2) (shape=(3, 3, 3), mean=0.83004 1e-17 erg / (Angstrom s spaxel cm2)); spectral_axis=<SpectralAxis
        (observer to target:
           radial_velocity=0.0 km / s
           redshift=0.0)
-      [5.00034537e-07, 5.00149688e-07, 5.00264865e-07] m>, uncertainty=InverseVariance([[[3449.242 , 2389.292 , 2225.105 ],
-                      [4098.7485, 2965.88  , 2632.497 ],
-                      [3589.92  , 2902.7622, 2292.3823]],
-    <BLANKLINE>
-                     [[3563.3342, 2586.58  , 2416.039 ],
-                      [4090.8855, 3179.1702, 2851.823 ],
-                      [4158.919 , 3457.0115, 2841.1965]],
-    <BLANKLINE>
-                     [[3684.6013, 3056.2   , 2880.6592],
-                      [3221.7888, 2801.3518, 2525.541 ],
-                      [3936.68  , 3461.534 , 3047.6135]]]))>
+      [5.00034537e-07 5.00149688e-07 5.00264865e-07] m> (length=3); uncertainty=InverseVariance)>
 
 
 Spectral Bounding Region
diff --git a/docs/spectrum1d.rst b/docs/spectrum1d.rst
index f29f864f2..f844a5fe5 100644
--- a/docs/spectrum1d.rst
+++ b/docs/spectrum1d.rst
@@ -64,7 +64,7 @@ a loader.
     >>> import urllib
     >>> specs = urllib.request.urlopen('https://data.sdss.org/sas/dr14/sdss/spectro/redux/26/spectra/0751/spec-0751-52251-0160.fits') # doctest: +REMOTE_DATA
     >>> Spectrum1D.read(specs, format="SDSS-III/IV spec") # doctest: +REMOTE_DATA
-    <Spectrum1D(flux=<Quantity [30.596626,...]...>
+    <Spectrum1D(flux=[30.59662628173828 ... 51.70271682739258] 1e-17 erg / (Angstrom s cm2) (shape=(3841,), mean=51.88042 1e-17 erg / (Angstrom s cm2)); spectral_axis=<SpectralAxis [3799.2686 3800.1426 3801.0188 ... 9193.905  9196.0205 9198.141 ] Angstrom> (length=3841); uncertainty=InverseVariance)>
 
 Note that the same spectrum could be more conveniently downloaded via
 astroquery, if the user has that package installed:
@@ -74,7 +74,7 @@ astroquery, if the user has that package installed:
      >>> from astroquery.sdss import SDSS  # doctest: +REMOTE_DATA
      >>> specs = SDSS.get_spectra(plate=751, mjd=52251, fiberID=160, data_release=14)  # doctest: +REMOTE_DATA
      >>> Spectrum1D.read(specs[0], format="SDSS-III/IV spec")  # doctest: +REMOTE_DATA
-     <Spectrum1D(flux=<Quantity [30.596626,...]...>
+     <Spectrum1D(flux=[30.59662628173828 ... 51.70271682739258] 1e-17 erg / (Angstrom s cm2) (shape=(3841,), mean=51.88042 1e-17 erg / (Angstrom s cm2)); spectral_axis=<SpectralAxis [3799.2686 3800.1426 3801.0188 ... 9193.905  9196.0205 9198.141 ] Angstrom> (length=3841); uncertainty=InverseVariance)>
 
 
 List of Loaders
@@ -326,69 +326,11 @@ value will apply to the lower bound input.
     >>> lower = [SpectralCoord(4.9, unit=u.um), SkyCoord(ra=205, dec=26, unit=u.deg)]
     >>> upper = [SpectralCoord(4.9, unit=u.um), SkyCoord(ra=205.5, dec=27.5, unit=u.deg)]
     >>> spec.crop(lower, upper)  # doctest: +IGNORE_WARNINGS +FLOAT_CMP
-    <Spectrum1D(flux=<Quantity [[[0.70861236],
-                [0.5663815 ],
-                [0.0606386 ],
-                [0.13811995],
-                [0.8974065 ]],
-    <BLANKLINE>
-               [[0.97618597],
-                [0.870499  ],
-                 [0.01522275],
-                 [0.59180312],
-                 [0.29160346]],
-    <BLANKLINE>
-               [[0.13274479],
-                [0.99381789],
-                [0.767089  ],
-                [0.73765335],
-                [0.36185756]],
-    <BLANKLINE>
-               [[0.09635759],
-                [0.38060021],
-                [0.63263223],
-                [0.60785285],
-                [0.4914904 ]],
-    <BLANKLINE>
-               [[0.83173506],
-                [0.1679683 ],
-                [0.39415721],
-                [0.25540459],
-                [0.39779061]],
-    <BLANKLINE>
-               [[0.27625437],
-                [0.90381118],
-                [0.65453332],
-                [0.64141404],
-                [0.36941242]],
-    <BLANKLINE>
-               [[0.90620901],
-                [0.41437874],
-                [0.71279457],
-                [0.41105785],
-                [0.7459662 ]],
-    <BLANKLINE>
-               [[0.57160214],
-                [0.45227951],
-                [0.29112881],
-                [0.44654224],
-                [0.65356417]],
-    <BLANKLINE>
-               [[0.60171961],
-                [0.46916617],
-                [0.3919754 ],
-                [0.01304226],
-                [0.32085301]],
-    <BLANKLINE>
-               [[0.00470581],
-                [0.54743546],
-                [0.24740782],
-                [0.77903598],
-                [0.63457146]]] Jy>, spectral_axis=<SpectralAxis
-       (observer to target:
-          radial_velocity=0.0 km / s
-          redshift=0.0)
-      [4.90049987e-06] m>)>
+    <Spectrum1D(flux=[[[0.708612359963129 ... 0.6345714580773677]]] Jy (shape=(10, 5, 1), mean=0.49653 Jy); spectral_axis=<SpectralAxis
+        (observer to target:
+           radial_velocity=0.0 km / s
+           redshift=0.0)
+      [4.90049987e-06] m> (length=1))>
 
 Collapsing
 ----------
@@ -415,11 +357,11 @@ spectral axis, or 'spatial', which will collapse along all non-spectral axes.
 
     >>> spec.mean(axis='spatial')  # doctest: +FLOAT_CMP
     <Spectrum1D(flux=<Quantity [0.37273938, 0.53843905, 0.61351648, 0.57311623, 0.44339915,
-               0.66084728, 0.45881921, 0.38715911, 0.39967185, 0.53257671] Jy>, spectral_axis=<SpectralAxis
+               0.66084728, 0.45881921, 0.38715911, 0.39967185, 0.53257671] Jy> (shape=(10,), mean=0.49803 Jy); spectral_axis=<SpectralAxis
        (observer to target:
           radial_velocity=0.0 km / s
           redshift=0.0)
-      [5000., 5001., 5002., 5003., 5004., 5005., 5006., 5007., 5008., 5009.] Angstrom>)>
+      [5000. 5001. 5002. ... 5007. 5008. 5009.] Angstrom> (length=10))>
 
 Note that in this case, the result of the collapse operation is a
 `~specutils.Spectrum1D` rather than an `astropy.units.Quantity`, because the
diff --git a/specutils/manipulation/model_replace.py b/specutils/manipulation/model_replace.py
index 4d86af725..54f627ace 100644
--- a/specutils/manipulation/model_replace.py
+++ b/specutils/manipulation/model_replace.py
@@ -65,7 +65,7 @@ def model_replace(spectrum, replace_region, model=10, extrapolation_treatment='d
     >>> input_spectrum = Spectrum1D(spectral_axis=wave_val * u.AA, flux=flux_val * u.mJy)
     >>> spline_knots = [3.5, 4.7, 6.8, 7.1] * u.AA
     >>> model_replace(input_spectrum, None, spline_knots)
-    <Spectrum1D(flux=<Quantity [ 2.,  4.,  6.,  8., 10., 12., 14., 16., 18., 20.] mJy>, spectral_axis=<SpectralAxis [ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.] Angstrom>)>
+    <Spectrum1D(flux=<Quantity [ 2.,  4.,  6.,  8., 10., 12., 14., 16., 18., 20.] mJy> (shape=(10,), mean=11.00000 mJy); spectral_axis=<SpectralAxis [ 1.  2.  3. ...  8.  9. 10.] Angstrom> (length=10))>
 
     """
     if extrapolation_treatment not in ('data_fill', 'zero_fill'):
diff --git a/specutils/manipulation/resample.py b/specutils/manipulation/resample.py
index ab1db3cf9..e208d2d14 100644
--- a/specutils/manipulation/resample.py
+++ b/specutils/manipulation/resample.py
@@ -76,7 +76,8 @@ class FluxConservingResampler(ResamplerBase):
     >>> resample_grid = [1, 5, 9, 13, 14, 17, 21, 22, 23]  *u.nm
     >>> fluxc_resample = FluxConservingResampler()
     >>> fluxc_resample(input_spectra, resample_grid)  # doctest: +FLOAT_CMP
-    <Spectrum1D(flux=<Quantity [  1. ,  3.  ,  6.  ,  7.  ,  6.25, 10.  , 20.  ,   nan,   nan] mJy>, spectral_axis=<SpectralAxis [ 1.,  5.,  9., 13., 14., 17., 21., 22., 23.] nm>)>
+    <Spectrum1D(flux=<Quantity [ 1.  ,  3.  ,  6.  ,  7.  ,  6.25, 10.  , 20.  ,   nan,   nan] mJy> (shape=(9,), mean=7.60714 mJy); spectral_axis=<SpectralAxis [ 1.  5.  9. ... 21. 22. 23.] nm> (length=9))>
+
     """
 
     def _fluxc_resample(self, input_bin_centers, output_bin_centers,
@@ -332,7 +333,7 @@ class LinearInterpolatedResampler(ResamplerBase):
     >>> resample_grid = [1, 5, 9, 13, 14, 17, 21, 22, 23] * u.nm
     >>> fluxc_resample = LinearInterpolatedResampler()
     >>> fluxc_resample(input_spectra, resample_grid)  # doctest: +FLOAT_CMP
-    <Spectrum1D(flux=<Quantity [ nan, 3.5 , 5.5 , 6.75, 6.5 , 9.5 ,  nan,  nan,  nan] mJy>, spectral_axis=<SpectralAxis [ 1.,  5.,  9., 13., 14., 17., 21., 22., 23.] nm>)>
+    <Spectrum1D(flux=<Quantity [ nan, 3.5 , 5.5 , 6.75, 6.5 , 9.5 ,  nan,  nan,  nan] mJy> (shape=(9,), mean=6.35000 mJy); spectral_axis=<SpectralAxis [ 1.  5.  9. ... 21. 22. 23.] nm> (length=9))>
 
     """
     def __init__(self, extrapolation_treatment='nan_fill'):
@@ -417,7 +418,7 @@ class SplineInterpolatedResampler(ResamplerBase):
     >>> fluxc_resample = SplineInterpolatedResampler()
     >>> fluxc_resample(input_spectra, resample_grid)  # doctest: +FLOAT_CMP
     <Spectrum1D(flux=<Quantity [       nan, 3.98808594, 6.94042969, 6.45869141, 5.89921875,
-               7.29736328,        nan,        nan,        nan] mJy>, spectral_axis=<SpectralAxis [ 1.,  5.,  9., 13., 14., 17., 21., 22., 23.] nm>)>
+               7.29736328,        nan,        nan,        nan] mJy> (shape=(9,), mean=6.11676 mJy); spectral_axis=<SpectralAxis [ 1.  5.  9. ... 21. 22. 23.] nm> (length=9))>
 
     """
     def __init__(self, extrapolation_treatment='nan_fill'):
diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py
index 24def21f3..7702b58f5 100644
--- a/specutils/spectra/spectrum1d.py
+++ b/specutils/spectra/spectrum1d.py
@@ -724,46 +724,50 @@ def __rsub__(self, other):
         return -1 * (self - other)
 
     def _format_array_summary(self, label, array):
-        if len(array) == 1:
-            mean = np.mean(array)
-            s = "{:17} [ {:.5} ],  mean={:.5}"
-            return s.format(label+':', array[0], array[-1], mean)
-        elif len(array) > 1:
-            mean = np.mean(array)
-            s = "{:17} [ {:.5}, ..., {:.5} ],  mean={:.5}"
-            return s.format(label+':', array[0], array[-1], mean)
+        array_str = np.array2string(array, threshold=8, prefix=label)
+        if len(array) >= 1:
+            mean = np.nanmean(array)
+            s = f"{label}{array_str} {array.unit},  mean={mean:.5f}"
+            return s
         else:
             return "{:17} [ ],  mean= n/a".format(label+':')
 
     def __str__(self):
         result = "Spectrum1D "
-        # Handle case of single value flux
-        if self.flux.ndim == 0:
-            result += "(length=1)\n"
-            return result + "flux:   {}".format(self.flux)
-
-        # Handle case of multiple flux arrays
         result += "(length={})\n".format(len(self.spectral_axis))
-        if self.flux.ndim > 1:
-            for i, flux in enumerate(self.flux):
-                label = 'flux{:2}'.format(i)
-                result += self._format_array_summary(label, flux) + '\n'
-        else:
-            result += self._format_array_summary('flux', self.flux) + '\n'
+
+        # Add Flux information
+        result += self._format_array_summary('Flux=', self.flux) + '\n'
+
         # Add information about spectral axis
-        result += self._format_array_summary('spectral axis', self.spectral_axis)
+        result += self._format_array_summary('Spectral Axis=', self.spectral_axis)
+
         # Add information about uncertainties if available
         if self.uncertainty:
-            result += "\nuncertainty:      [ {}, ..., {} ]".format(
-                self.uncertainty[0], self.uncertainty[-1])
+            result += (f'\nUncertainty={type(self.uncertainty).__name__} '
+                       f'({np.array2string(self.uncertainty.array, threshold=8)}'
+                       f' {self.uncertainty.unit})')
+
         return result
 
     def __repr__(self):
-        inner_str = "flux={}, spectral_axis={}".format(repr(self.flux),
-                                                       repr(self.spectral_axis))
+        flux_str  = "flux="
+        if (self.flux.ndim == 1 and self.flux.size <= 10) or self.flux.size <= 20:
+            flux_str += repr(self.flux)
+        else:
+            flux_summary = f"{self.flux.value.flat[0]} ... {self.flux.value.flat[-1]}"
+            flux_str = flux_str + "[" * self.flux.ndim + flux_summary + "]" * self.flux.ndim
+            flux_str += f" {self.flux.unit}"
+
+        flux_str += f" (shape={self.flux.shape}, mean={np.nanmean(self.flux):.5f}); "
+        spectral_axis_str = (repr(self.spectral_axis).split("[")[0] +
+                             np.array2string(self.spectral_axis, threshold=8) +
+                             f" {self.spectral_axis.unit}>")
+        spectral_axis_str = f"spectral_axis={spectral_axis_str} (length={len(self.spectral_axis)})"
+        inner_str = (flux_str + spectral_axis_str)
 
         if self.uncertainty is not None:
-            inner_str += ", uncertainty={}".format(repr(self.uncertainty))
+            inner_str += f"; uncertainty={self.uncertainty.__class__.__name__}"
 
         result = "<Spectrum1D({})>".format(inner_str)
 
diff --git a/specutils/tests/test_spectrum1d.py b/specutils/tests/test_spectrum1d.py
index eea2a43ad..470f6d516 100644
--- a/specutils/tests/test_spectrum1d.py
+++ b/specutils/tests/test_spectrum1d.py
@@ -437,7 +437,7 @@ def test_repr():
     result = repr(spec_with_unc)
     assert result.startswith('<Spectrum1D(flux=<Quantity [')
     assert 'spectral_axis=<SpectralAxis' in result
-    assert 'uncertainty=StdDevUncertainty(' in result
+    assert 'uncertainty=StdDevUncertainty' in result
 
 
 def test_str():
@@ -449,42 +449,39 @@ def test_str():
     lines = result.split('\n')
     flux = spec.flux
     sa = spec.spectral_axis
-    assert len(lines) == 3
-    assert lines[1].startswith('flux:')
-    assert '[ {:.5}, ..., {:.5} ]'.format(flux[0], flux[-1]) in lines[1]
-    assert 'mean={:.5}'.format(np.mean(flux)) in lines[1]
-    assert lines[2].startswith('spectral axis:')
-    assert '[ {:.5}, ..., {:.5} ]'.format(sa[0], sa[-1]) in lines[2]
-    assert 'mean={:5}'.format(np.mean(sa)) in lines[2]
+    assert len(lines) == 4
+    assert lines[1].startswith('Flux=')
+    assert 'mean={:.5f}'.format(np.nanmean(flux)) in lines[2]
+    assert lines[3].startswith('Spectral Axis=')
+    assert f'mean={np.nanmean(sa):.5f}' in lines[3]
 
     # Test string representation with uncertainty
     spec_with_uncertainty = Spectrum1D(spectral_axis=np.linspace(100, 1000, 10) * u.nm,
                                        flux=np.random.random(10) * u.Jy,
-                                       uncertainty=np.random.random(10))
+                                       uncertainty=StdDevUncertainty(np.random.random(10)))
     result = str(spec_with_uncertainty)
     lines = result.split('\n')
     unc = spec_with_uncertainty.uncertainty
     print(spec_with_uncertainty)
-    assert len(lines) == 4
-    assert lines[3].startswith('uncertainty')
-    assert '[ {}, ..., {} ]'.format(unc[0], unc[-1]) in lines[3]
+    assert len(lines) == 5
+    assert lines[4].startswith('Uncertainty')
+    assert f'StdDevUncertainty ([{unc.array[0]:.7f}' in lines[4]
 
     # Test string representation with multiple flux
     spec_multi_flux = Spectrum1D(spectral_axis=np.linspace(100, 1000, 10) * u.nm,
                                  flux=np.random.random((3,10)) * u.Jy)
     result = str(spec_multi_flux)
     lines = result.split('\n')
-    assert len(lines) == 5
-    for i, line in enumerate(lines[1:4]):
-        assert line.startswith('flux{:2}:'.format(i))
+    assert len(lines) == 8
+    assert lines[1].startswith('Flux=')
 
     # Test string representation with single-dimensional flux
     spec_single_flux = Spectrum1D([1] * u.Jy, [0] * u.nm)
     result = str(spec_single_flux)
     assert result == \
 """Spectrum1D (length=1)
-flux:             [ 1.0 Jy ],  mean=1.0 Jy
-spectral axis:    [ 0.0 nm ],  mean=0.0 nm"""
+Flux=[1.] Jy,  mean=1.00000 Jy
+Spectral Axis=[0.] nm,  mean=0.00000 nm"""
 
 
 def test_equivalencies():