diff --git a/modules/icing/pirep_goes.py b/modules/icing/pirep_goes.py
index 204c47fec55ab5572d543fe1d2dded15d245ea1e..4c63fe76d47a8817462cabc532aa282c1f2f8614 100644
--- a/modules/icing/pirep_goes.py
+++ b/modules/icing/pirep_goes.py
@@ -800,6 +800,12 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask,
 def fov_extract(outfile='/Users/tomrink/fovs_out.h5'):
     ice_times = []
     icing_int_s = []
+    ice_lons = []
+    ice_lats = []
+
+    no_ice_times = []
+    no_ice_lons = []
+    no_ice_lats = []
 
     h5_s_icing = []
     h5_s_no_icing = []
@@ -819,6 +825,8 @@ def fov_extract(outfile='/Users/tomrink/fovs_out.h5'):
         num_obs = len(times)
 
         icing_int = f['icing_intensity'][:]
+        lons = f['longitude'][:]
+        lats = f['latitude'][:]
 
         cld_mask = f['cloud_mask'][:, 10:30, 10:30]
         cld_mask = cld_mask.reshape((num_obs, -1))
@@ -843,8 +851,10 @@ def fov_extract(outfile='/Users/tomrink/fovs_out.h5'):
                 dat = f[ds_name][i, 10:30, 10:30].flatten()
                 icing_data_dct[ds_name].append(dat[k_idxs])
 
-            ice_flags = np.full(len(k_idxs), icing_int[i])
-            icing_int_s.append(ice_flags)
+            icing_int_s.append(np.full(len(k_idxs), icing_int[i]))
+            ice_times.append(np.full(len(k_idxs), times[i]))
+            ice_lons.append(np.full(len(k_idxs), lons[i]))
+            ice_lats.append(np.full(len(k_idxs), lats[i]))
 
         print(fname)
 
@@ -853,6 +863,9 @@ def fov_extract(outfile='/Users/tomrink/fovs_out.h5'):
         icing_data_dct[ds_name] = np.concatenate(lst)
 
     icing_int_s = np.concatenate(icing_int_s)
+    ice_times = np.concatenate(ice_times)
+    ice_lons = np.concatenate(ice_lons)
+    ice_lats = np.concatenate(ice_lats)
 
     num_no_ice = 0
     for fidx in range(len(no_icing_files)):
@@ -862,6 +875,8 @@ def fov_extract(outfile='/Users/tomrink/fovs_out.h5'):
 
         times = f['time']
         num_obs = len(times)
+        lons = f['longitude']
+        lats = f['latitude']
 
         cld_mask = f['cloud_mask'][:, 10:30, 10:30]
         cld_mask = cld_mask.reshape((num_obs, -1))
@@ -881,6 +896,9 @@ def fov_extract(outfile='/Users/tomrink/fovs_out.h5'):
             else:
                 k_idxs = k_idxs[0:len(k_idxs)]
             num_no_ice += len(k_idxs)
+            no_ice_times.append(np.full(len(k_idxs), times[i]))
+            no_ice_lons.append(np.full(len(k_idxs), lons[i]))
+            no_ice_lats.append(np.full(len(k_idxs), lats[i]))
 
             for ds_name in train_params_day:
                 dat = f[ds_name][i, 10:30, 10:30].flatten()
@@ -892,19 +910,29 @@ def fov_extract(outfile='/Users/tomrink/fovs_out.h5'):
         lst = no_icing_data_dct[ds_name]
         no_icing_data_dct[ds_name] = np.concatenate(lst)
     no_icing_int_s = np.full(num_no_ice, -1)
-
-    ds_indexes = np.arange(num_ice + num_no_ice)
-    np.random.shuffle(ds_indexes)
+    no_ice_times = np.concatenate(no_ice_times)
+    no_ice_lons = np.concatenate(no_ice_lons)
+    no_ice_lats = np.concatenate(no_ice_lats)
 
     icing_intensity = np.concatenate([icing_int_s, no_icing_int_s])
+    icing_times = np.concatenate([ice_times, no_ice_times])
+    icing_lons = np.concatenate([ice_lons, no_ice_lons])
+    icing_lats = np.concatenate([ice_lats, no_ice_lats])
 
     data_dct = {}
     for ds_name in train_params_day:
         data_dct[ds_name] = np.concatenate([icing_data_dct[ds_name], no_icing_data_dct[ds_name]])
 
+    # apply shuffle indexes
+    ds_indexes = np.arange(num_ice + num_no_ice)
+    np.random.shuffle(ds_indexes)
+
     for ds_name in train_params_day:
         data_dct[ds_name] = data_dct[ds_name][ds_indexes]
     icing_intensity = icing_intensity[ds_indexes]
+    icing_times = icing_times[ds_indexes]
+    icing_lons = icing_lons[ds_indexes]
+    icing_lats = icing_lats[ds_indexes]
 
     h5f_expl = h5py.File(a_clvr_file, 'r')
     h5f_out = h5py.File(outfile, 'w')
@@ -917,6 +945,18 @@ def fov_extract(outfile='/Users/tomrink/fovs_out.h5'):
     icing_int_ds = h5f_out.create_dataset('icing_intensity', data=icing_intensity, dtype='i4')
     icing_int_ds.attrs.create('long_name', data='From PIREP. -1:No Icing, 1:Trace, 2:Light, 3:Light Moderate, 4:Moderate, 5:Moderate Severe, 6:Severe')
 
+    time_ds = h5f_out.create_dataset('time', data=icing_times, dtype='f4')
+    time_ds.attrs.create('units', data='seconds since 1970-1-1 00:00:00')
+    time_ds.attrs.create('long_name', data='PIREP time')
+
+    lon_ds = h5f_out.create_dataset('longitude', data=icing_lons, dtype='f4')
+    lon_ds.attrs.create('units', data='degrees_east')
+    lon_ds.attrs.create('long_name', data='PIREP longitude')
+
+    lat_ds = h5f_out.create_dataset('latitude', data=icing_lats, dtype='f4')
+    lat_ds.attrs.create('units', data='degrees_north')
+    lat_ds.attrs.create('long_name', data='PIREP latitude')
+
     # copy relevant attributes
     for ds_name in train_params_day:
         h5f_ds = h5f_out[ds_name]
@@ -937,6 +977,13 @@ def fov_extract(outfile='/Users/tomrink/fovs_out.h5'):
 
 def box_extract(outfile='/Users/tomrink/box_out.h5'):
     icing_int_s = []
+    ice_time_s = []
+    no_ice_time_s = []
+    ice_lon_s = []
+    no_ice_lon_s = []
+    ice_lat_s = []
+    no_ice_lat_s = []
+
 
     h5_s_icing = []
     h5_s_no_icing = []
@@ -951,6 +998,8 @@ def box_extract(outfile='/Users/tomrink/box_out.h5'):
 
         times = f['time'][:]
         num_obs = len(times)
+        lons = f['longitude']
+        lats = f['latitude']
 
         icing_int = f['icing_intensity'][:]
 
@@ -959,6 +1008,9 @@ def box_extract(outfile='/Users/tomrink/box_out.h5'):
                 dat = f[ds_name][i, 12:28, 12:28]
                 icing_data_dct[ds_name].append(dat)
             icing_int_s.append(icing_int[i])
+            ice_time_s.append(times[i])
+            ice_lon_s.append(lons[i])
+            ice_lat_s.append(lats[i])
 
         print(fname)
 
@@ -966,6 +1018,9 @@ def box_extract(outfile='/Users/tomrink/box_out.h5'):
         lst = icing_data_dct[ds_name]
         icing_data_dct[ds_name] = np.stack(lst, axis=0)
     icing_int_s = np.array(icing_int_s)
+    ice_time_s = np.array(ice_time_s)
+    ice_lon_s = np.array(ice_lon_s)
+    ice_lat_s = np.array(ice_lat_s)
     num_ice = icing_int_s.shape[0]
 
     num_no_ice = 0
@@ -976,12 +1031,17 @@ def box_extract(outfile='/Users/tomrink/box_out.h5'):
 
         times = f['time']
         num_obs = len(times)
+        lons = f['longitude']
+        lats = f['latitude']
 
         for i in range(num_obs):
             for ds_name in train_params_day:
                 dat = f[ds_name][i, 12:28, 12:28]
                 no_icing_data_dct[ds_name].append(dat)
             num_no_ice += 1
+            no_ice_time_s.append(times[i])
+            no_ice_lon_s.append(lons[i])
+            no_ice_lat_s.append(lats[i])
 
         print(fname)
 
@@ -989,19 +1049,29 @@ def box_extract(outfile='/Users/tomrink/box_out.h5'):
         lst = no_icing_data_dct[ds_name]
         no_icing_data_dct[ds_name] = np.stack(lst, axis=0)
     no_icing_int_s = np.full(num_no_ice, -1)
-
-    ds_indexes = np.arange(num_ice + num_no_ice)
-    np.random.shuffle(ds_indexes)
+    no_ice_time_s = np.array(no_ice_time_s)
+    no_ice_lon_s = np.array(no_ice_lon_s)
+    no_ice_lat_s = np.array(no_ice_lat_s)
 
     icing_intensity = np.concatenate([icing_int_s, no_icing_int_s])
+    icing_times = np.concatenate([ice_time_s, no_ice_time_s])
+    icing_lons = np.concatenate([ice_lon_s, no_ice_lon_s])
+    icing_lats = np.concatenate([ice_lat_s, no_ice_lat_s])
 
     data_dct = {}
     for ds_name in train_params_day:
         data_dct[ds_name] = np.concatenate([icing_data_dct[ds_name], no_icing_data_dct[ds_name]])
 
+    # Do shuffle
+    ds_indexes = np.arange(num_ice + num_no_ice)
+    np.random.shuffle(ds_indexes)
+
     for ds_name in train_params_day:
         data_dct[ds_name] = data_dct[ds_name][ds_indexes]
     icing_intensity = icing_intensity[ds_indexes]
+    icing_times = icing_times[ds_indexes]
+    icing_lons = icing_lons[ds_indexes]
+    icing_lats = icing_lats[ds_indexes]
 
     h5f_expl = h5py.File(a_clvr_file, 'r')
     h5f_out = h5py.File(outfile, 'w')
@@ -1014,6 +1084,18 @@ def box_extract(outfile='/Users/tomrink/box_out.h5'):
     icing_int_ds = h5f_out.create_dataset('icing_intensity', data=icing_intensity, dtype='i4')
     icing_int_ds.attrs.create('long_name', data='From PIREP. -1:No Icing, 1:Trace, 2:Light, 3:Light Moderate, 4:Moderate, 5:Moderate Severe, 6:Severe')
 
+    time_ds = h5f_out.create_dataset('time', data=icing_times, dtype='f4')
+    time_ds.attrs.create('units', data='seconds since 1970-1-1 00:00:00')
+    time_ds.attrs.create('long_name', data='PIREP time')
+
+    lon_ds = h5f_out.create_dataset('longitude', data=icing_lons, dtype='f4')
+    lon_ds.attrs.create('units', data='degrees_east')
+    lon_ds.attrs.create('long_name', data='PIREP longitude')
+
+    lat_ds = h5f_out.create_dataset('latitude', data=icing_lats, dtype='f4')
+    lat_ds.attrs.create('units', data='degrees_north')
+    lat_ds.attrs.create('long_name', data='PIREP latitude')
+
     # copy relevant attributes
     for ds_name in train_params_day:
         h5f_ds = h5f_out[ds_name]