diff --git a/doc/src/Section_howto.txt b/doc/src/Section_howto.txt
index 3ca9c3d1990ecc21f978a0016fb8b31517fc5af6..78478c3fe6ffb85593b6a1634110db8bb5f5a0b1 100644
--- a/doc/src/Section_howto.txt
+++ b/doc/src/Section_howto.txt
@@ -1989,6 +1989,11 @@ Both methods are thus a means to extract or assign (overwrite) any
 peratom quantities within LAMMPS.  See the extract() method in the
 src/atom.cpp file for a list of valid per-atom properties.  New names
 could easily be added if the property you want is not listed.
+A special treatment is applied for accessing image flags via the
+"image" property. Image flags are stored in a packed format with all
+three image flags stored in a single integer. When signaling to access
+the image flags as 3 individual values per atom instead of 1, the data
+is transparently packed or unpacked by the library interface.
 
 The lammps_create_atoms() function takes a list of N atoms as input
 with atom types and coords (required), an optionally atom IDs and
diff --git a/doc/src/Section_python.txt b/doc/src/Section_python.txt
index b994a564027e8b86b3230cdfc67c69f9b754c9c2..ff26d18e069d818ff5a3cead58180545a9b2690d 100644
--- a/doc/src/Section_python.txt
+++ b/doc/src/Section_python.txt
@@ -698,7 +698,12 @@ method in the src/atom.cpp file for a list of valid names.  Again, new
 names could easily be added if the property you want is missing.  The
 vector can be used via normal Python subscripting.  If atom IDs are
 not consecutively ordered within LAMMPS, a None is returned as
-indication of an error.
+indication of an error. A special treatment is applied for image flags
+stored in the "image" property. All three image flags are stored in
+a packed format in a single integer, so count would be 1 to retrieve
+that integer, however also a count value of 3 can be used and then
+the image flags will be unpacked into 3 individual integers, ordered
+in a similar fashion as coordinates. 
 
 Note that the data structure gather_atoms("x") returns is different
 from the data structure returned by extract_atom("x") in four ways.
@@ -727,6 +732,10 @@ corresponding properties for each atom inside LAMMPS.  This requires
 LAMMPS to have its "map" option enabled; see the
 "atom_modify"_atom_modify.html command for details.  If it is not, or
 if atom IDs are not consecutively ordered, no coordinates are reset.
+Similar as for gather_atoms() a special treatment is applied for image
+flags, which can be provided in packed (count = 1) or unpacked (count = 3)
+format and in the latter case, they will be packed before applied to
+atoms.
 
 The array of coordinates passed to scatter_atoms() must be a ctypes
 vector of ints or doubles, allocated and initialized something like
diff --git a/src/library.cpp b/src/library.cpp
index f9c5ad4892f1bb6d383dc4f6fd056f49b8168b0e..892818ff9ca873509226b13ab4b67ebd9b4e4b3a 100644
--- a/src/library.cpp
+++ b/src/library.cpp
@@ -773,7 +773,9 @@ void lammps_gather_atoms(void *ptr, char *name,
     if (type == 0) {
       int *vector = NULL;
       int **array = NULL;
-      if (count == 1) vector = (int *) vptr;
+      const int imgunpack = (count == 3) && (strcmp(name,"image") == 0);
+
+      if ((count == 1) || imgunpack) vector = (int *) vptr;
       else array = (int **) vptr;
 
       int *copy;
@@ -786,7 +788,15 @@ void lammps_gather_atoms(void *ptr, char *name,
       if (count == 1)
         for (i = 0; i < nlocal; i++)
           copy[tag[i]-1] = vector[i];
-      else
+      else if (imgunpack) {
+        for (i = 0; i < nlocal; i++) {
+          offset = count*(tag[i]-1);
+          const int image = vector[i];
+          copy[offset++] = (image & IMGMASK) - IMGMAX;
+          copy[offset++] = ((image >> IMGBITS) & IMGMASK) - IMGMAX;
+          copy[offset++] = ((image >> IMG2BITS) & IMGMASK) - IMGMAX;
+        }
+      } else
         for (i = 0; i < nlocal; i++) {
           offset = count*(tag[i]-1);
           for (j = 0; j < count; j++)
@@ -874,7 +884,9 @@ void lammps_scatter_atoms(void *ptr, char *name,
     if (type == 0) {
       int *vector = NULL;
       int **array = NULL;
-      if (count == 1) vector = (int *) vptr;
+      const int imgpack = (count == 3) && (strcmp(name,"image") == 0);
+
+      if ((count == 1) || imgpack) vector = (int *) vptr;
       else array = (int **) vptr;
       int *dptr = (int *) data;
 
@@ -882,6 +894,15 @@ void lammps_scatter_atoms(void *ptr, char *name,
         for (i = 0; i < natoms; i++)
           if ((m = lmp->atom->map(i+1)) >= 0)
             vector[m] = dptr[i];
+      } else if (imgpack) {
+        for (i = 0; i < natoms; i++)
+          if ((m = lmp->atom->map(i+1)) >= 0) {
+            offset = count*i;
+            int image = dptr[offset++] + IMGMAX;
+            image += (dptr[offset++] + IMGMAX) << IMGBITS;
+            image += (dptr[offset++] + IMGMAX) << IMG2BITS;
+            vector[m] = image;
+          }
       } else {
         for (i = 0; i < natoms; i++)
           if ((m = lmp->atom->map(i+1)) >= 0) {