diff --git a/pyproject.toml b/pyproject.toml
index deff8ef7..c5642263 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -18,7 +18,7 @@ dependencies = [
     "scipy>=1.5.0",
     "matplotlib>=3.2.0",
     "astropy>=5.1.0",
-    "photutils>=1.0.0",
+    "photutils>=1.10.0",
     "poppy>=1.0.0",
     "pysiaf>=0.19.1",
     "soc_roman_tools>=0.1.0",
diff --git a/requirements.txt b/requirements.txt
index 33754585..c280f922 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,10 +1,10 @@
-astropy==5.3.4
-ipython==8.17.2
-matplotlib==3.8.1
-numpy==1.26.1
-photutils==1.9.0
+astropy==6.0.0
+ipython==8.18.1
+matplotlib==3.8.2
+numpy==1.26.2
+photutils==1.10.0
 poppy>=1.0.0
 pysiaf==0.21.0
-scipy==1.11.3
-synphot==1.2.1
+scipy==1.11.4
+synphot==1.3.post0
 astroquery==0.4.6
diff --git a/webbpsf/gridded_library.py b/webbpsf/gridded_library.py
index 6a852797..d1b95040 100644
--- a/webbpsf/gridded_library.py
+++ b/webbpsf/gridded_library.py
@@ -1,13 +1,13 @@
-from collections import OrderedDict
 import itertools
 import os
+from collections import OrderedDict
 
 import astropy.convolution
+import numpy as np
 from astropy.io import fits
 from astropy.nddata import NDData
-import numpy as np
-import poppy
 
+import poppy
 import webbpsf.detectors
 
 
@@ -399,7 +399,8 @@ def create_grid(self):
             model_list.append(model)
 
             if self.save is True:
-                self.writeto(psf_arr, meta, det)
+                # model.data is sorted with xygrid, psf_arr is not. Store the sorted data
+                self.writeto(model.data, meta, det)
 
         # If only 1 detector, only return that 1 object. Else, return list of objects
         if len(self.detector_list) == 1:
diff --git a/webbpsf/tests/test_psfgrid.py b/webbpsf/tests/test_psfgrid.py
index a41f1800..013d62b4 100644
--- a/webbpsf/tests/test_psfgrid.py
+++ b/webbpsf/tests/test_psfgrid.py
@@ -33,7 +33,7 @@ def test_compare_to_calc_psf_oversampled():
 
     # Pull one of the PSFs out of the grid
     psfnum = 1
-    loc = grid.meta["grid_xypos"][psfnum]
+    loc = grid.grid_xypos[psfnum]
     locy = int(float(loc[1]) - 0.5)
     locx = int(float(loc[0]) - 0.5)
     gridpsf = grid.data[psfnum, :, :]
@@ -68,7 +68,7 @@ def test_compare_to_calc_psf_detsampled():
 
     # Pull one of the PSFs out of the grid
     psfnum = 1
-    loc = grid.meta["grid_xypos"][psfnum]
+    loc = grid.grid_xypos[psfnum]
     locy = int(float(loc[1]))
     locx = int(float(loc[0]))
     gridpsf = grid.data[psfnum, :, :]
@@ -142,8 +142,8 @@ def test_one_psf():
     scalefactor = oversample**2 # normalization as used internally in GriddedPSFModel; see #302
 
 
-    assert grid1.meta["grid_xypos"] == [(1023, 1023)], "Center position not as expected"  # the default is the center of the NIS aperture
-    assert grid2.meta["grid_xypos"] == [(10, 0)], "Corner position not as expected" # it's in (x,y)
+    assert np.all(grid1.grid_xypos == [[1023, 1023]]), "Center position not as expected"  # the default is the center of the NIS aperture
+    assert np.all(grid2.grid_xypos == [(10, 0)]), "Corner position not as expected" # it's in (x,y)
     # check for near-equality of the PSFs computed both ways,
     #  but ignore the outer few pixels rows and columns, for which boundary wrapping leads to imperfect equality
     #  depending on the order of the convolutions. Ignore 2 rows/cols on either side based on oversample value
@@ -209,7 +209,7 @@ def test_saving(tmpdir):
         # Check meta data
         model = utils.to_griddedpsfmodel(infile)
         assert model.meta.keys() == grid.meta.keys()
-        assert model.meta["grid_xypos"] == grid.meta["grid_xypos"]
+        assert np.all(model.grid_xypos == grid.grid_xypos)
         assert model.meta["oversampling"] == grid.meta["oversampling"]
 
     # Remove temporary directory
@@ -249,7 +249,7 @@ def test_wfi():
 
     # Pull one of the PSFs out of the grid
     psfnum = 1
-    loc = grid.meta["grid_xypos"][psfnum]
+    loc = grid.grid_xypos[psfnum]
     locy = int(float(loc[1])-0.5)
     locx = int(float(loc[0])-0.5)
     gridpsf = grid.data[psfnum, :, :]
@@ -263,6 +263,4 @@ def test_wfi():
 
     # Compare to make sure they are in fact the same PSF
     assert gridpsf.shape == calcpsf.shape, "Shape mismatch"
-    assert np.allclose(gridpsf, convpsf*scalefactor), "Data values not as expected"
-
-
+    assert np.allclose(gridpsf, convpsf*scalefactor), "Data values not as expected"
\ No newline at end of file