Merge pull request #87 from ninahakansson/split_files

Split VGAC files at midnight

Author: ninahakansson

.github/workflows/ci.yaml

@@ -10,10 +10,10 @@ jobs:
       fail-fast: true
       matrix:
         os: ["ubuntu-latest", "macos-latest"]
-        python-version: ["3.9", "3.10", "3.11"]
+        python-version: ["3.9", "3.11", "3.12"]
         experimental: [false]
         include:
-          - python-version: "3.10"
+          - python-version: "3.12"
             os: "ubuntu-latest"
             experimental: true
 

bin/vgac2pps.py

@@ -53,8 +53,10 @@
     parser.add_argument('-on', '--orbit_number', type=int, nargs='?',
                         required=False, default=0,
                         help="Orbit number (default is 00000).")
-
+    parser.add_argument('--don_split_files_at_midnight', action='store_true',
+                        help="Don't split files at midnight, keep as one level1c file.")
     options = parser.parse_args()
     process_one_scene(options.files, options.out_dir, engine=options.nc_engine,
                       all_channels=options.all_channels, pps_channels=options.pps_channels,
-                      orbit_n=options.orbit_number, as_noaa19=options.as_noaa19, avhrr_channels=options.avhrr_channels)
+                      orbit_n=options.orbit_number, as_noaa19=options.as_noaa19, avhrr_channels=options.avhrr_channels,
+                      split_files_at_midnight = not options.don_split_files_at_midnight)

continuous_integration/environment.yaml

@@ -7,7 +7,7 @@ dependencies:
   - h5py
   - python-geotiepoints
   - mock
-  - numpy
+  - numpy<2.0.0
   - satpy>0.41.1
   - pyspectral
   - h5netcdf

level1c4pps/__init__.py

@@ -553,8 +553,9 @@ def compose_filename(scene, out_path, instrument, band=None):
         end_time = band.attrs['end_time']
     platform_name = scene.attrs['platform']
     orbit_number = int(scene.attrs['orbit_number'])
+    out_path_with_dates = start_time.strftime(out_path)
     filename = os.path.join(
-        out_path,
+        out_path_with_dates,
         "S_NWC_{:s}_{:s}_{:05d}_{:s}Z_{:s}Z.nc".format(
             instrument,
             platform_name_to_use_in_filename(platform_name),

level1c4pps/eumgacfdr2pps_lib.py

@@ -35,9 +35,9 @@
                          logger,
                          get_header_attrs, convert_angles)
 from satpy.utils import debug_on
-from distutils.version import LooseVersion
+from packaging.version import Version
 
-if LooseVersion(satpy.__version__) < LooseVersion('0.24.0'):
+if Version(satpy.__version__) < Version('0.24.0'):
     debug_on()
     raise ImportError("'eumgac2pps' writer requires satpy 0.24.0 or greater")
 # import xarray as xr

level1c4pps/gac2pps_lib.py

@@ -39,8 +39,17 @@
                          get_header_attrs, convert_angles)
 import logging
 
+from packaging.version import Version
 logger = logging.getLogger('gac2pps')
 
+if Version(np.__version__) >= Version('2.0.0'):
+    if Version(pygac.__version__) == Version('1.7.3'):
+        raise ImportError("pygac 1.7.3 requires numpy < 2.0.0")
+    else:
+        logger.warning("pygac 1.7.3 requires numpy < 2.0.0 or greater")
+
+
+
 BANDNAMES = ['1', '2', '3', '3a', '3b', '4', '5']
 
 REFL_BANDS = ['1', '2', '3a']

level1c4pps/slstr2pps_lib.py

@@ -36,9 +36,9 @@
 import pyspectral  # testing that pyspectral is available # noqa: F401
 import logging
 from satpy.utils import debug_on
+from packaging.version import Version
 
-from distutils.version import LooseVersion
-if LooseVersion(satpy.__version__) < LooseVersion('0.22.1'):
+if Version(satpy.__version__) < Version('0.22.1'):
     raise ImportError("'slstr2pps' requires satpy 0.22.1 or greater")
 
 debug_on()

level1c4pps/tests/VGAC_VNPP02MOD_A2012365_2304_n06095_K005.nc

@@ -156,3 +156,34 @@ def test_process_one_scene_n19(self):
 
         np.testing.assert_equal(pps_nc.__dict__["platform"], "vgac20")
         self.assertTrue(np.abs(pps_nc.variables['image1'][0,0,0] - pps_nc_viirs.variables['image1'][0,0,0])>0.01)
+        
+    def test_process_one_scene_midnight(self):
+        """Test process one scene for one example file."""
+
+        vgac2pps.process_one_scene(
+            ['./level1c4pps/tests/VGAC_VNPP02MOD_A2012365_2304_n06095_K005.nc'],
+            out_path='./level1c4pps/tests/',
+        )
+        filename = './level1c4pps/tests/S_NWC_viirs_npp_00000_20121230T2359563Z_20121230T2359599Z.nc'
+        # written with hfnetcdf read with NETCDF4 ensure compatability
+        pps_nc = netCDF4.Dataset(filename, 'r', format='NETCDF4')  # Check compatability implicitly
+
+        for key in ['start_time', 'end_time', 'history', 'instrument',
+                    'orbit_number', 'platform',
+                    'sensor', 'source']:
+            if key not in pps_nc.__dict__.keys():
+                print("Missing in attributes:", key)
+            self.assertTrue(key in pps_nc.__dict__.keys())
+
+        expected_vars = ['satzenith', 'azimuthdiff',
+                         'satazimuth', 'sunazimuth', 'sunzenith',
+                         'lon', 'lat',
+                         'image1', 'image2', 'image3', 'image4', 'image5',
+                         'image6', 'image7', 'image8', 'image9',
+                         'scanline_timestamps', 'time', 'time_bnds']
+        for var in expected_vars:
+            self.assertTrue(var in pps_nc.variables.keys())
+
+        print(pps_nc.variables['image1'].shape)
+
+        np.testing.assert_equal(pps_nc.variables['image1'].shape, (1, 7, 801))

level1c4pps/tests/test_vgac2pps.py

@@ -30,6 +30,7 @@
 from level1c4pps import (get_encoding, compose_filename,
                          set_header_and_band_attrs_defaults,
                          rename_latitude_longitude,
+                         dt64_to_datetime,
                          update_angle_attributes, get_header_attrs,
                          convert_angles)
 import pyspectral  # testing that pyspectral is available # noqa: F401
@@ -59,7 +60,7 @@
 
 
 REFL_BANDS = ["M01", "M02", "M03", "M04", "M05", "M06", "M07", "M08",
-              "M09"]
+              "M09", "M10", "M11", "I01", "I02", "I03"]
 
 MBAND_PPS = ["M05", "M07", "M09", "M10", "M11", "M12", "M14", "M15", "M16"]
 
@@ -172,12 +173,83 @@ def set_header_and_band_attrs(scene, orbit_n=0):
         scene[band].attrs['sun_zenith_angle_correction_applied'] = 'True'
     return nimg
 
+def midnight_scene(scene):
+    """Check if scene passes midnight."""
+    start_date = scene["M05"].attrs["start_time"].strftime("%Y%m%d")
+    end_date = scene["M05"].attrs["end_time"].strftime("%Y%m%d")
+    if start_date == end_date:
+        return False
+    return True
+
+
+def get_midnight_line_nr(scene):
+    """Find midnight_line, start_time and new end_time."""
+    start_date = scene["M05"].attrs["start_time"].strftime("%Y-%m-%d")
+    end_date = scene["M05"].attrs["end_time"].strftime("%Y-%m-%d")
+    start_fine_search = len(scene['scanline_timestamps']) - 1  # As default start the fine search from end of time array
+    for ind in range(0, len(scene['scanline_timestamps']), 100):
+        # Search from the beginning in large chunks (100) and break when we
+        # pass midnight.
+        dt_obj = dt64_to_datetime(scene['scanline_timestamps'].values[:][ind])
+        date_linei = dt_obj.strftime("%Y-%m-%d")
+        if date_linei == end_date:
+            # We just passed midnight stop and search backwards for exact line.
+            start_fine_search = ind
+            break
+    for indj in range(start_fine_search, start_fine_search - 100, -1):
+        # Midnight is in one of the previous 100 lines.
+        dt_obj = dt64_to_datetime(scene['scanline_timestamps'].values[:][indj])
+        date_linei = dt_obj.strftime("%Y-%m-%d") 
+        if date_linei == start_date:
+            # We just passed midnight this is the last line for previous day.
+            midnight_linenr = indj
+            break
+    return  midnight_linenr
+
+
+
+def set_exact_time_and_crop(scene, start_line, end_line, time_key='scanline_timestamps'):
+    """Crop datasets and update start_time end_time objects."""
+    if start_line is None:
+        start_line = 0
+    if end_line is None:
+        end_line = len(scene[time_key]) - 1
+    start_time_dt64 = scene[time_key].values[start_line]
+    end_time_dt64 = scene[time_key].values[end_line]
+    start_time = dt64_to_datetime(start_time_dt64)
+    end_time = dt64_to_datetime(end_time_dt64)
+    for ds in BANDNAMES + ANGLE_NAMES + ['latitude', 'longitude', 'scanline_timestamps']:
+        if ds in scene and 'nscn' in scene[ds].dims:
+            scene[ds] = scene[ds].isel(nscn=slice(start_line, end_line + 1))
+            try:
+                # Update scene attributes to get the filenames right
+                scene[ds].attrs['start_time'] = start_time
+                scene[ds].attrs['end_time'] = end_time
+            except TypeError:
+                pass
+    if start_time_dt64 != scene[time_key].values[0]:
+        raise ValueError
+    if end_time_dt64 != scene[time_key].values[-1]:
+        raise ValueError
+        
+def split_scene_at_midnight(scene):
+    """Split scenes at midnight."""
+    if midnight_scene(scene):
+        midnight_linenr = get_midnight_line_nr(scene)
+        scene1 = scene.copy()
+        scene2 = scene.copy()
+        set_exact_time_and_crop(scene1, None, midnight_linenr)
+        set_exact_time_and_crop(scene2, midnight_linenr + 1, None)
+        return [scene1, scene2]
+    return [scene]
+
 
 def process_one_scene(scene_files, out_path, engine='h5netcdf',
-                      all_channels=False, pps_channels=False, orbit_n=0, as_noaa19=False, avhrr_channels=False):
+                      all_channels=False, pps_channels=False, orbit_n=0, as_noaa19=False, avhrr_channels=False,
+                      split_files_at_midnight=True):
     """Make level 1c files in PPS-format."""
     tic = time.time()
-    scn_ = Scene(
+    scn_in = Scene(
         reader='viirs_vgac_l1c_nc',
         filenames=scene_files)
 
@@ -192,41 +264,46 @@ def process_one_scene(scene_files, out_path, engine='h5netcdf',
     if avhrr_channels:
         MY_MBAND = MBAND_AVHRR
 
-    scn_.load(MY_MBAND
-              + ANGLE_NAMES
-              # + ['M12_LUT', 'M13_LUT', 'M15_LUT', 'M16_LUT']
-              + ['latitude', 'longitude', 'scanline_timestamps'])
-
-    # one ir channel
-    irch = scn_['M15']
-
-    # Set header and band attributes
-    set_header_and_band_attrs(scn_, orbit_n=orbit_n)
-
-    # Rename longitude, latitude to lon, lat.
-    rename_latitude_longitude(scn_)
-
-    # Convert angles to PPS
-    convert_angles(scn_, delete_azimuth=False)
-    update_angle_attributes(scn_, irch)
-
-    # Adjust to noaa19 with sbafs from KG
-    sensor = "viirs"
-    if as_noaa19:
-        sensor = "avhrr"
-        convert_to_noaa19(scn_)
-
-    filename = compose_filename(scn_, out_path, instrument=sensor, band=irch)
-    encoding = get_encoding_viirs(scn_)
-
-    scn_.save_datasets(writer='cf',
-                       filename=filename,
-                       header_attrs=get_header_attrs(scn_, band=irch, sensor=sensor),
-                       engine=engine,
-                       include_lonlats=False,
-                       flatten_attrs=True,
-                       encoding=encoding)
-    print("Saved file {:s} after {:3.1f} seconds".format(
-        os.path.basename(filename),
-        time.time()-tic))
-    return filename
+    scn_in.load(MY_MBAND
+                + ANGLE_NAMES
+                # + ['M12_LUT', 'M13_LUT', 'M15_LUT', 'M16_LUT']
+                + ['latitude', 'longitude', 'scanline_timestamps'])
+    if split_files_at_midnight:
+        scenes = split_scene_at_midnight(scn_in)
+    else:
+        scenes = [scn_in]
+    filenames = []
+    for scn_ in scenes:
+        # one ir channel
+        irch = scn_['M15']
+
+        # Set header and band attributes
+        set_header_and_band_attrs(scn_, orbit_n=orbit_n)
+
+        # Rename longitude, latitude to lon, lat.
+        rename_latitude_longitude(scn_)
+    
+        # Convert angles to PPS
+        convert_angles(scn_, delete_azimuth=False)
+        update_angle_attributes(scn_, irch)
+        # Adjust to noaa19 with sbafs from KG
+        sensor = "viirs"
+        if as_noaa19:
+            sensor = "avhrr"
+            convert_to_noaa19(scn_)
+
+        filename = compose_filename(scn_, out_path, instrument=sensor, band=irch)
+        encoding = get_encoding_viirs(scn_)
+
+        scn_.save_datasets(writer='cf',
+                           filename=filename,
+                           header_attrs=get_header_attrs(scn_, band=irch, sensor=sensor),
+                           engine=engine,
+                           include_lonlats=False,
+                           flatten_attrs=True,
+                           encoding=encoding)
+        print("Saved file {:s} after {:3.1f} seconds".format(
+            os.path.basename(filename),
+            time.time()-tic))
+        filenames.append(filename)
+    return filenames

level1c4pps/vgac2pps_lib.py

@@ -70,7 +70,7 @@
       data_files=[],
       zip_safe=False,
       use_scm_version=True,
-      python_requires='>=2.7,!=3.0.*,!=3.1.*,!=3.2.*,!=3.3.*',
+      python_requires='>=3.7',
       install_requires=requires,
       test_suite='level1c4pps.tests.suite',
       )