diff --git a/atmPy/data_archives/NOAA_ESRL_GMD_GRAD/surfrad/surfrad.py b/atmPy/data_archives/NOAA_ESRL_GMD_GRAD/surfrad/surfrad.py index b34a544..1ec0a58 100644 --- a/atmPy/data_archives/NOAA_ESRL_GMD_GRAD/surfrad/surfrad.py +++ b/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 diff --git a/atmPy/radiation/observations/spectral_irradiance.py b/atmPy/radiation/observations/spectral_irradiance.py index f4800dc..073c372 100644 --- a/atmPy/radiation/observations/spectral_irradiance.py +++ b/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