diff --git a/src/gf3d/bin/gf3d_exec.py b/src/gf3d/bin/gf3d_exec.py
index b70492a..13ec9d6 100644
--- a/src/gf3d/bin/gf3d_exec.py
+++ b/src/gf3d/bin/gf3d_exec.py
@@ -114,15 +114,25 @@ def query_stations(databasename: str, debug: bool, local: bool):
             GFM = GFManager(db_files)
             GFM.load_header_variables()
 
-            for net, sta, lat, lon, bur in zip(
-                    GFM.networks,
-                    GFM.stations,
-                    GFM.latitudes,
-                    GFM.longitudes,
-                    GFM.burials):
-                print(f"[{net},{sta},{lat},{lon},{bur}]")
+            # Get station info:
+            stations = GFM.get_stations()
+            Nsta = len(stations)
 
+            # Print copyable string array
+            endofline = '\n'
+
+            print('[', end="")
+            for _i, (net, sta, lat, lon, bur) in enumerate(stations):
+
+                # Remove newline character for the last station
+                if _i == Nsta-1:
+                    endofline = ''
+                print(f"['{net}', '{sta}', {lat}, {lon}, {bur}],",
+                      end=endofline)
+
+            print(']')
     else:
+
         from gf3d.client import GF3DClient
         gfcl = GF3DClient(databasename, debug=debug)
         stations = gfcl.stations_avail()
diff --git a/src/gf3d/seismograms.py b/src/gf3d/seismograms.py
index b1f649f..f20c3eb 100644
--- a/src/gf3d/seismograms.py
+++ b/src/gf3d/seismograms.py
@@ -800,6 +800,43 @@ def load_header_variables(self):
                             'Dont have buffer elements for doing the ')
                         self.do_adjacency_search = False
 
+    def get_stations(self):
+
+        # Check if KDtree is loaded
+        if not self.header:
+            self.load_header_variables()
+
+        if self.subset:
+            networks = self.networks
+            stations = self.stations
+            latitudes = self.latitudes
+            longitudes = self.longitudes
+            burials = self.burials
+
+        else:
+
+            networks = []
+            stations = []
+            latitudes = []
+            longitudes = []
+            burials = []
+            # Get number of files
+            logger.info("Opening h5 files ...")
+            with contextlib.ExitStack() as stack:
+                # Open Stack of files
+                dbs = [stack.enter_context(h5py.File(fname, 'r'))
+                       for fname in self.db]
+
+                for db in dbs:
+                    networks.append(db['Network'][()].decode())
+                    stations.append(db['Station'][()].decode())
+                    latitudes.append(db['latitude'][()])
+                    longitudes.append(db['longitude'][()])
+                    burials.append(db['burial'][()])
+
+        return [sta_tup for sta_tup in zip(
+            networks, stations, latitudes, longitudes, burials)]
+
     def get_elements(self, lat, lon, depth, dist_in_km=125.0, NGLL=5, threading: bool = True):
 
         if self.subset:
@@ -1805,7 +1842,8 @@ def get_mt_frechet_station(
 
         # Heaviside STF to reproduce SPECFEM stf
         _, stf_r = create_stf(0, 400.0, self.header['nsteps'],
-                              self.header['dt'], hdur_diff, cutoff=None, gaussian=False, lpfilter='butter')
+                              self.header['dt'], hdur_diff, cutoff=None,
+                              gaussian=False, lpfilter='butter')
 
         STF_R = fft.fft(stf_r, n=NP2)