diff --git a/modules/deeplearning/icing_cnn.py b/modules/deeplearning/icing_cnn.py
index 260f36637f41cd884420403a3bb691f174f6bdf8..6b32892a988720d5e52ca81503ff90af392b11ef 100644
--- a/modules/deeplearning/icing_cnn.py
+++ b/modules/deeplearning/icing_cnn.py
@@ -1013,16 +1013,20 @@ def run_restore_static(filename_l1b, filename_l2, ckpt_dir_s_path):
 
 def run_evaluate_static(h5f, ckpt_dir_s_path, prob_thresh=0.5, satellite='GOES16', domain='FD'):
     data_dct, ll, cc = make_for_full_domain_predict(h5f, name_list=train_params, domain=domain)
+    num_elems = len(cc)
+    num_lines = len(ll)
 
     ckpt_dir_s = os.listdir(ckpt_dir_s_path)
 
+    nav = get_navigation(satellite, domain)
+
     prob_s = []
     for ckpt in ckpt_dir_s:
         ckpt_dir = ckpt_dir_s_path + ckpt
         if not os.path.isdir(ckpt_dir):
             continue
         nn = IcingIntensityNN()
-        nn.setup_eval_pipeline(data_dct, len(ll) * len(cc))
+        nn.setup_eval_pipeline(data_dct, num_lines * num_elems)
         nn.build_model()
         nn.build_training()
         nn.build_evaluation()
@@ -1040,23 +1044,26 @@ def run_evaluate_static(h5f, ckpt_dir_s_path, prob_thresh=0.5, satellite='GOES16
         preds = np.where(probs > prob_thresh, 1, 0)
     else:
         preds = np.argmax(probs, axis=1)
+    preds_2d = preds.reshape((len(ll), len(cc)))
+    cc = np.array(cc)
+    ll = np.array(ll)
+    x_rad = cc * nav.CFAC + nav.COFF
+    y_rad = ll * nav.LFAC + nav.LOFF
 
     ll, cc = np.meshgrid(ll, cc, indexing='ij')
     cc = cc.flatten()
     ll = ll.flatten()
+    lon_s, lat_s = nav.lc_to_earth(cc, ll)
+    lons_2d = lon_s.reshape((num_lines, num_elems))
+    lats_2d = lat_s.reshape((num_lines, num_elems))
+
     ice_mask = preds == 1
     ice_cc = cc[ice_mask]
     ice_ll = ll[ice_mask]
 
-    nav = get_navigation(satellite, domain)
-    ice_lons = []
-    ice_lats = []
-    for k in range(ice_cc.shape[0]):
-        lon, lat = nav.lc_to_earth(ice_cc[k], ice_ll[k])
-        ice_lons.append(lon)
-        ice_lats.append(lat)
-
-    return ice_lons, ice_lats
+    ice_lons, ice_lats = nav.lc_to_earth(ice_cc, ice_ll)
+
+    return ice_lons, ice_lats, preds_2d, lons_2d, lats_2d, x_rad, y_rad
 
 
 if __name__ == "__main__":
diff --git a/modules/icing/pirep_goes.py b/modules/icing/pirep_goes.py
index f781fec0dc229d5d4c763947ae3169a658ba8ec5..af2f4e4cfc839644945fedbc1090ad58b7f92713 100644
--- a/modules/icing/pirep_goes.py
+++ b/modules/icing/pirep_goes.py
@@ -4,7 +4,7 @@ import pickle
 import matplotlib.pyplot as plt
 import os
 from util.util import get_time_tuple_utc, GenericException, add_time_range_to_filename, is_night, is_day, \
-    check_oblique, make_times, find_bin_index, get_timestamp, homedir
+    check_oblique, make_times, find_bin_index, get_timestamp, homedir, write_icing_file
 from util.plot import make_icing_image
 from aeolus.datasource import CLAVRx, CLAVRx_VIIRS, GOESL1B, CLAVRx_H08
 import h5py
@@ -1883,7 +1883,9 @@ def run_make_images(clvrx_dir='/Users/tomrink/data/clavrx/RadC/', ckpt_dir_s_pat
         else:
             keep_lons = None
             keep_lats = None
-        ice_lons, ice_lats = run_evaluate_static(h5f, ckpt_dir_s_path=ckpt_dir_s_path, prob_thresh=prob_thresh, satellite=satellite, domain=domain)
+
+        ice_lons, ice_lats, preds_2d, lons_2d, lats_2d, x_rad, y_rad = run_evaluate_static(h5f, ckpt_dir_s_path=ckpt_dir_s_path, prob_thresh=prob_thresh, satellite=satellite, domain=domain)
         make_icing_image(h5f, ice_lons, ice_lats, clvrx_str_time, satellite, domain, ice_lons_vld=keep_lons, ice_lats_vld=keep_lats, extent=extent)
+        write_icing_file(clvrx_str_time, preds_2d, x_rad, y_rad, lons_2d, lats_2d)
         print('Done: ', clvrx_str_time)
         h5f.close()
diff --git a/modules/util/util.py b/modules/util/util.py
index 9953dfcc255dd7899c7722e30cbd03db80645ac6..e474ad8d440c424aee2fc004859eabd4e4fe2c38 100644
--- a/modules/util/util.py
+++ b/modules/util/util.py
@@ -557,4 +557,57 @@ def make_for_full_domain_predict(h5f, name_list=None, satellite='GOES16', domain
         grd_dct[ds_name] = np.stack(grd_dct_n[ds_name])
 
     return grd_dct, ll, cc
-# -------------------------------------------------------------------------------------------
\ No newline at end of file
+# -------------------------------------------------------------------------------------------
+
+
+def write_icing_file(clvrx_str_time, preds, x, y, lons, lats):
+    outfile_name = homedir + 'icing_prediction_'+clvrx_str_time+'.h5'
+    h5f_out = h5py.File(outfile_name, 'w')
+
+    icing_pred_ds = h5f_out.create_dataset('icing_prediction', data=preds, dtype='i2')
+    icing_pred_ds.attrs.create('coordinates', data='y x')
+    icing_pred_ds.attrs.create('grid_mapping', data='Projection')
+
+    lon_ds = h5f_out.create_dataset('longitude', data=lons, dtype='f4')
+    lon_ds.attrs.create('units', data='degrees_east')
+    lon_ds.attrs.create('long_name', data='icing prediction longitude')
+
+    lat_ds = h5f_out.create_dataset('latitude', data=lats, dtype='f4')
+    lat_ds.attrs.create('units', data='degrees_north')
+    lat_ds.attrs.create('long_name', data='icing prediction latitude')
+
+    proj_ds = h5f_out.create_dataset('Projection', data=0, dtype='b')
+    proj_ds.attrs.create('long_name', data='Himawari Imagery Projection')
+    proj_ds.attrs.create('sweep_angle_axis', data='y')
+    proj_ds.attrs.create('units', data='radians')
+    proj_ds.attrs.create('semi_major_axis', data=6378.137)
+    proj_ds.attrs.create('semi_minor_axis', data=6356.7523)
+    proj_ds.attrs.create('inverse_flattening', data=298.257)
+    proj_ds.attrs.create('perspective_point_height', data=35785.863)
+    proj_ds.attrs.create('latitude_of_projection_origin', data=0.0)
+    proj_ds.attrs.create('longitude_of_projection_origin', data=140.7)
+    proj_ds.attrs.create('CFAC', data=20466275)
+    proj_ds.attrs.create('LFAC', data=20466275)
+    proj_ds.attrs.create('COFF', data=2750.5)
+    proj_ds.attrs.create('LOFF', data=2750.5)
+
+    if x is not None:
+        x_ds = h5f_out.create_dataset('x', data=x, dtype='f8')
+        x_ds.attrs.create('units', data='radians')
+        x_ds.attrs.create('standard_name', data='projection_x_coordinate')
+        x_ds.attrs.create('long_name', data='GOES PUG W-E fixed grid viewing angle')
+        x_ds.attrs.create('scale_factor', data=5.58879902955962e-05)
+        x_ds.attrs.create('add_offset', data=-0.153719917308037)
+        x_ds.attrs.create('CFAC', data=20466275)
+        x_ds.attrs.create('COFF', data=2750.5)
+
+        y_ds = h5f_out.create_dataset('y', data=y, dtype='f8')
+        y_ds.attrs.create('units', data='radians')
+        y_ds.attrs.create('standard_name', data='projection_y_coordinate')
+        y_ds.attrs.create('long_name', data='GOES PUG S-N fixed grid viewing angle')
+        y_ds.attrs.create('scale_factor', data=-5.58879902955962e-05)
+        y_ds.attrs.create('add_offset', data=0.153719917308037)
+        y_ds.attrs.create('LFAC', data=20466275)
+        y_ds.attrs.create('LOFF', data=2750.5)
+
+    h5f_out.close()