new: read mfrsr raw files

update: some adjustements to the spectal_irradiance module

atmPy/data_archives/NOAA_ESRL_GMD_GRAD/surfrad/surfrad.py

@@ -11,6 +11,9 @@
 import atmPy.radiation.observations.spectral_irradiance as atmradobs
 import atmPy.data_archives.NOAA_ESRL_GMD_GRAD.cal_facility.lab as atmcucf
 
+import subprocess
+import io
+
 # from atmPy.general import measurement_site as _measurement_site
 # from atmPy.radiation import solar as _solar
 
@@ -253,6 +256,102 @@ class FilterFunctionError(Exception):
                          }
 
 
+
+class FileCorruptError(Exception):
+    """Custom exception raised when a file is detected as corrupt."""
+    def __init__(self, filename, message="File is corrupt or unreadable."):
+        self.filename = filename
+        self.message = f"{message} ({filename})"
+        super().__init__(self.message)
+
+def read_raw(path2file, read_header_only = False):
+    """
+    Opens a MFRSR or MFR raw file. Those are the files that have an extension 
+    like .mtm, .xmd, .rsr. Surfrad files are stored here:
+    /nfs/grad/Inst/MFR/SURFRAD/{site}/mfrsr/raw/
+
+    Parameters
+    ----------
+    path2file : TYPE
+        DESCRIPTION.
+    read_header_only : TYPE, optional
+        DESCRIPTION. The default is False.
+
+    Returns
+    -------
+    TYPE
+        DESCRIPTION.
+
+    """
+    command = f"tu -d joe -H {path2file.absolute()}"
+    # print(command)
+    # out = subprocess.check_output(command, shell=True).decode()
+
+
+    process = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
+
+    out, err = process.communicate()
+
+    out = out.decode()
+    err = err.decode()
+
+    if err.split()[0] == 'Error':
+        raise FileCorruptError(filename = path2file.as_posix(), message = err)
+
+
+    #### the header
+    header = out.split('\n')[0]    
+    if read_header_only:
+        return header
+
+    #### format the data
+    df = _pd.read_csv(io.StringIO(out), sep = r'\s+', header = None, skiprows=1)
+    df.index = df.apply(lambda row: _pd.to_datetime('19000101') + _pd.to_timedelta(row[0] - 1, 'd'), axis = 1)
+    df.index.name = 'datetime'
+    df = df.where(df!=-9999)
+
+
+    #### assign collumns
+    attrs = dict(site_longitude = - float(header.split()[4]),
+                 site_latitude = float(header.split()[3]),
+                 site_elevation = 0,
+                 site = 'TMP',
+                 site_name = 'unknown',
+                 calibrated_irradiance = False,
+                 calibrated_cosine = False,
+                 info = 'This is the raw file',
+                 file_type = int(header.split()[0]),
+                 serial_no = header.split()[1],
+
+                )
+
+    di = 7
+    si = 2
+    alltime = df.iloc[:,si: si+di]
+    si = 11
+    global_horizontal = df.iloc[:,si: si+di]
+    si =  18
+    diffuse_horizontal = df.iloc[:,si: si+di]
+    si =  25
+    direct = df.iloc[:,si: si+di]
+
+    for dft in [alltime, global_horizontal, diffuse_horizontal, direct]:
+        dft.columns = range(7)
+        dft.columns.name = 'channel'
+
+    ds = _xr.Dataset()
+    # ds['sun_position'] = sun_position
+    ds['alltime'] = alltime
+    ds['global_horizontal'] = global_horizontal
+    ds['diffuse_horizontal'] = diffuse_horizontal
+    ds['direct_normal'] = direct
+    ds.attrs = attrs
+
+    obs = atmradobs.CombinedGlobalDiffuseDirect(ds)
+    return obs
+
+
+
 def _path2files(path2base_folder = '/nfs/grad/surfrad/aod/', site = 'bon', window = ('2020-07-31 00:00:00', '2020-10-30 23:00:00'), product = 'aod'):
     path2base_folder = _pl.Path(path2base_folder)
 
@@ -298,7 +397,7 @@ def _read_data(fname, UTC = False, header=None):
 
     # dateparse = lambda x: _pd.datetime.strptime(x, "%d:%m:%Y %H:%M:%S")
     df = _pd.read_csv(fname, skiprows=header['header_size'],
-                      sep ='\s+',
+                      sep =r'\s+',
                      # delim_whitespace=True,
                      #                      na_values=['N/A'],
                      #                   parse_dates={'times': [0, 1]},
@@ -513,7 +612,7 @@ def read_readme(p2r):
     def read_data(p2f, colnames):
         df = _pd.read_csv(p2f, 
                           # delim_whitespace=True,
-                      sep ='\s+',
+                      sep =r'\s+',
                     names = colnames,
                    )
 
@@ -809,7 +908,7 @@ def read_ccc(p2f, site = None, path2cal_file = None, verbose = False, raisefilte
         timezoneadjust = 0
     df = _pd.read_csv(p2f,
                 names =cols,
-                      sep ='\s+',
+                      sep =r'\s+',
                 # delim_whitespace = True,
                 skiprows = skiprows)
     # return df

atmPy/radiation/observations/spectral_irradiance.py

@@ -15,7 +15,7 @@
 import atmPy.data_archives.NOAA_ESRL_GMD_GRAD.surfrad.surfrad as atmsrf
 import atmPy.data_archives.NOAA_ESRL_GMD_GRAD.baseline.baseline as atmbl
 import sqlite3
-
+import matplotlib.pyplot as _plt
 import copy
 
 class RadiationDatabase(object):
@@ -392,9 +392,23 @@ def add_clearsky_parameters2database(self, database):
 
 
 class SolarIrradiation(object):
-    def __init__(self, dataset):
+    def __init__(self, dataset, site = None):
         self.dataset = self.unify_variables(dataset)
         self.clearsky = ClearSky(self)
+
+        if isinstance(site, type(None)):
+            assert('site' in dataset.attrs.keys()), 'If site is None, then the dataset has to have lat,lon,site, site_name, attributes'
+            self.site = atmms.Station(lat= dataset.attrs['site_latitude'], 
+                                      lon = dataset.attrs['site_longitude'], 
+                                      alt = dataset.attrs['site_elevation'], 
+                                      name = dataset.attrs['site_name'], 
+                                      abbreviation = dataset.attrs['site'],)
+        else:
+            self.site = site
+
+
+        self._sun_position = None
+
 
     def unify_variables(self, dataset):
         """Seach for variable names containing global, diffuse and direct and 
@@ -427,6 +441,11 @@ def register_file_in_database(self, database, overwrite = False):
                              product_name = ds.product_name, 
                              product_version = ds.product_version, 
                              overwrite = overwrite)
+    @property
+    def sun_position(self):
+        if isinstance(self._sun_position, type(None)):
+            self._sun_position = self.site.get_sun_position(self.dataset.datetime)
+        return self._sun_position
 
 class GlobalHorizontalIrradiation(SolarIrradiation):
     def __init__(self, dataset):
@@ -442,21 +461,21 @@ def __init__(self, dataset,
                  langley_fit_settings = None,
                  calibration_strategy = 'johns',
                  metdata = 'surfrad'):
-        super().__init__(dataset)
+        super().__init__(dataset, site = site)
         self.raw_data = dataset #this is not exactly raw, it is cosine corrected voltag readings, thus, un-calibrated irradiances
-        if isinstance(site, type(None)):
-            assert('site' in dataset.attrs.keys()), 'If site is None, then the dataset has to have lat,lon,site, site_name, attributes'
-            self.site = atmms.Station(lat= dataset.attrs['site_latitude'], 
-                                      lon = dataset.attrs['site_longitude'], 
-                                      alt = dataset.attrs['site_elevation'], 
-                                      name = dataset.attrs['site_name'], 
-                                      abbreviation = dataset.attrs['site'],)
-        else:
-            self.site = site
+        # if isinstance(site, type(None)):
+        #     assert('site' in dataset.attrs.keys()), 'If site is None, then the dataset has to have lat,lon,site, site_name, attributes'
+        #     self.site = atmms.Station(lat= dataset.attrs['site_latitude'], 
+        #                               lon = dataset.attrs['site_longitude'], 
+        #                               alt = dataset.attrs['site_elevation'], 
+        #                               name = dataset.attrs['site_name'], 
+        #                               abbreviation = dataset.attrs['site'],)
+        # else:
+        #     self.site = site
         self.langley_fit_settings = langley_fit_settings
         self.settings_calibration = calibration_strategy #
         self.settings_metdata = metdata
-        self._sun_position = None
+        # self._sun_position = None
         self._am = None
         self._pm = None
         self._transmission = None
@@ -928,11 +947,11 @@ def transmission(self):
                 self._transmission = raw_data/calib_interp_secorr
         return self._transmission
 
-    @property
-    def sun_position(self):
-        if isinstance(self._sun_position, type(None)):
-            self._sun_position = self.site.get_sun_position(self.raw_data.datetime)
-        return self._sun_position
+    # @property
+    # def sun_position(self):
+    #     if isinstance(self._sun_position, type(None)):
+    #         self._sun_position = self.site.get_sun_position(self.raw_data.datetime)
+    #     return self._sun_position
 
     @property
     def langley_am(self):
@@ -1057,6 +1076,30 @@ def __init__(self, dataset):
         self.global_horizontal_irradiation = GlobalHorizontalIrradiation(dataset)
         self.diffuse_horizontal_irradiation = DiffuseHorizontalIrradiation(dataset)
         self.direct_normal_irradiation = DirectNormalIrradiation(dataset)
+
+    def plot_overview(self, channel = 0, ax = None):
+
+        if isinstance(ax, type(None)):
+            f, a= _plt.subplots()    
+        else:
+            a = ax
+            f = a.get_figure()
+
+        dssel = self.dataset.sel(channel = 0)
+        dssel.alltime.plot(ax = a, label = 'alltime')
+
+        dssel.global_horizontal.plot(ax = a, label = 'global_horizontal')
+        dssel.diffuse_horizontal.plot(ax = a, label = 'diffuse_horizontal')
+        dssel.direct_normal.plot(ax = a, label = 'direct')
+
+        at = a.twinx()
+        self.sun_position.elevation.plot(ax = at, color = 'black', ls = '--')
+        at.set_ylim(top = 0.9, bottom = 0)
+
+        # a.set_xlim(left = pd.to_datetime('20220103 14:00:00'))
+        a.grid()
+        a.legend()
+        return f,a