diff --git a/.coveragerc b/.coveragerc new file mode 100644 index 0000000000..527b44f8a6 --- /dev/null +++ b/.coveragerc @@ -0,0 +1,2 @@ +[run] +omit =./act/tests/*, ./act/*version*py diff --git a/.travis.yml b/.travis.yml index 7c0949d742..e4bf8742e2 100644 --- a/.travis.yml +++ b/.travis.yml @@ -8,22 +8,23 @@ env: matrix: include: - - python: 3.6 + - python: 3.7 env: - - PYTHON_VERSION="3.6" + - PYTHON_VERSION="3.7" - DOC_BUILD="true" - - python: 3.7 + - python: 3.8 sudo: yes dist: xenial env: - - PYTHON_VERSION="3.7" + - PYTHON_VERSION="3.8" - DOC_BUILD="true" install: - source continuous_integration/install.sh - pip install pytest-cov - pip install coveralls + - pip install metpy script: - - eval xvfb-run pytest --cov=act/ + - eval xvfb-run pytest --mpl --cov=act/ --cov-config=.coveragerc - flake8 --max-line-length=115 --ignore=F401,E402,W504,W605 after_success: - coveralls diff --git a/CREATING_ENVIRONMENTS.rst b/CREATING_ENVIRONMENTS.rst index a6451867a4..7f7d801f2a 100644 --- a/CREATING_ENVIRONMENTS.rst +++ b/CREATING_ENVIRONMENTS.rst @@ -74,7 +74,7 @@ do this step while the environment is activate:: Another way to create a conda environment is by doing it from scratch using the conda create command. An example of this:: - conda create -n act_env -c conda-forge python=3.7 numpy pandas astral + conda create -n act_env -c conda-forge python=3.7 numpy pandas scipy matplotlib dask xarray After activating the environment with:: diff --git a/MANIFEST.in b/MANIFEST.in index 8030c3d862..65598ad098 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -10,10 +10,13 @@ recursive-exclude * *.py[co] recursive-include act/plotting *.txt recursive-include docs *.rst conf.py Makefile make.bat -recursive-include act/tests/data *.cdf *.nc *.data +recursive-include act/tests/data * + include versioneer.py include act/_version.py +include act/utils/conf/de421.bsp + # If including data files in the package, add them like: # include path/to/data_file diff --git a/README.rst b/README.rst index 0fa318ed94..dd5b65ad8e 100644 --- a/README.rst +++ b/README.rst @@ -4,7 +4,7 @@ Atmospheric data Community Toolkit (ACT) |AnacondaCloud| |Travis| |Coveralls| -|CondaDownloads| |Zenodo| +|CondaDownloads| |Zenodo| |ARM| .. |AnacondaCloud| image:: https://anaconda.org/conda-forge/act-atmos/badges/version.svg :target: https://anaconda.org/conda-forge/act-atmos @@ -21,11 +21,15 @@ Atmospheric data Community Toolkit (ACT) .. |Coveralls| image:: https://coveralls.io/repos/github/ARM-DOE/ACT/badge.svg :target: https://coveralls.io/github/ARM-DOE/ACT +.. |ARM| image:: https://img.shields.io/badge/Sponsor-ARM-blue.svg?colorA=00c1de&colorB=00539c + :target: https://www.arm.gov/ -Python toolkit for working with atmospheric time-series datasets of varying dimensions. The toolkit is meant to have functions for every part of the scientific process; discovery, IO, quality control, corrections, retrievals, visualization, and analysis. This toolkit is meant to be a community platform for sharing code with the goal of reducing duplication of effort and better connecting the science community with programs such as the `Atmospheric Radiation Measurement (ARM) User Facility `_. Overarching development goals will be updated on a regular basis as part of the `Roadmap `_. +The Atmospheric data Community Toolkit (ACT) is an open source Python toolkit for working with atmospheric time-series datasets of varying dimensions. The toolkit is meant to have functions for every part of the scientific process; discovery, IO, quality control, corrections, retrievals, visualization, and analysis. It is meant to be a community platform for sharing code with the goal of reducing duplication of effort and better connecting the science community with programs such as the `Atmospheric Radiation Measurement (ARM) User Facility `_. Overarching development goals will be updated on a regular basis as part of the `Roadmap `_ . -* Free software: 3-clause BSD license +|act| + +.. |act| image:: ./docs/source/act_plots.png Important Links ~~~~~~~~~~~~~~~ @@ -34,19 +38,22 @@ Important Links * Examples: https://arm-doe.github.io/ACT/source/auto_examples/index.html * Issue Tracker: https://github.com/ARM-DOE/ACT/issues +Citing +~~~~~~ + +If you use ACT to prepare a publication, please cite the DOI listed in the badge above, which is updated with every version release to ensure that contributors get appropriate credit. DOI is provided through Zenodo. + Dependencies ~~~~~~~~~~~~ +* `xarray `_ * `NumPy `_ * `SciPy `_ * `matplotlib `_ -* `xarray `_ -* `astral `_ +* `skyfield `_ * `pandas `_ * `dask `_ * `Pint `_ -* `Cartopy `_ -* `Boto3 `_ * `PyProj `_ * `Proj `_ * `Six `_ @@ -55,7 +62,9 @@ Dependencies Optional Dependencies ~~~~~~~~~~~~~~~~~~~~~ -* `MPL2NC `_ For reading binary MPL data. +* `MPL2NC `_ Reading binary MPL data. +* `Cartopy `_ Mapping and geoplots +* `MetPy `_ >= V1.0 Skew-T plotting and some stabilities indices calculations Installation ~~~~~~~~~~~~ @@ -138,7 +147,7 @@ Testing After installation, you can launch the test suite from outside the source directory (you will need to have pytest installed):: - $ pytest --pyargs act + $ pytest --mpl --pyargs act In-place installs can be tested using the `pytest` command from within the source directory. diff --git a/act/discovery/get_armfiles.py b/act/discovery/get_armfiles.py index 031444e21e..8973111887 100644 --- a/act/discovery/get_armfiles.py +++ b/act/discovery/get_armfiles.py @@ -36,6 +36,11 @@ def download_data(username, token, datastream, current working directory with the same name as *datastream* to place the files in. + Returns + ------- + files : list + Returns list of files retrieved + Notes ----- This programmatic interface allows users to query and automate @@ -107,7 +112,12 @@ def download_data(username, token, datastream, output_dir = os.path.join(os.getcwd(), datastream) # not testing, response is successful and files were returned + if response_body_json is None: + print("ARM Data Live Webservice does not appear to be functioning") + return [] + num_files = len(response_body_json["files"]) + file_names = [] if response_body_json["status"] == "success" and num_files > 0: for fname in response_body_json['files']: if time is not None: @@ -125,6 +135,9 @@ def download_data(username, token, datastream, # create file and write bytes to file with open(output_file, 'wb') as open_bytes_file: open_bytes_file.write(urlopen(save_data_url).read()) + file_names.append(output_file) else: print("No files returned or url status error.\n" "Check datastream name, start, and end date.") + + return file_names diff --git a/act/io/__init__.py b/act/io/__init__.py index 90b93fe272..0bf33c5634 100644 --- a/act/io/__init__.py +++ b/act/io/__init__.py @@ -20,3 +20,4 @@ from . import armfiles from . import csvfiles from . import mpl +from . import noaagml diff --git a/act/io/armfiles.py b/act/io/armfiles.py index a96a131fe1..556a4d71e6 100644 --- a/act/io/armfiles.py +++ b/act/io/armfiles.py @@ -14,8 +14,9 @@ import numpy as np import urllib import json -from enum import Flag, auto import copy +import act.utils as utils +import warnings def read_netcdf(filenames, concat_dim='time', return_None=False, @@ -32,7 +33,7 @@ def read_netcdf(filenames, concat_dim='time', return_None=False, Name of file(s) to read. concat_dim : str Dimension to concatenate files along. Default value is 'time.' - return_none : bool, optional + return_None : bool, optional Catch IOError exception when file not found and return None. Default is False. combine : str @@ -134,6 +135,7 @@ def read_netcdf(filenames, concat_dim='time', return_None=False, arm_ds[var_name].astype(desired_time_precision), arm_ds[var_name].attrs)}) arm_ds[var_name] = temp_ds[var_name] + temp_ds.close() # If time_offset is in file try to convert base_time as well if var_name == 'time_offset': @@ -160,13 +162,17 @@ def read_netcdf(filenames, concat_dim='time', return_None=False, not np.issubdtype(arm_ds['time'].values.dtype, np.datetime64) and not type(arm_ds['time'].values[0]).__module__.startswith('cftime.')): # Use microsecond precision to create time since epoch. Then convert to datetime64 - time = (arm_ds['base_time'].values * 1000000 + - arm_ds['time'].values * 1000000.).astype('datetime64[us]') + if arm_ds['base_time'].values == arm_ds['time_offset'].values[0]: + time = arm_ds['time_offset'].values + else: + time = (arm_ds['base_time'].values + + arm_ds['time_offset'].values * 1000000.).astype('datetime64[us]') # Need to use a new Dataset creation to correctly index time for use with # .group and .resample methods in Xarray Datasets. temp_ds = xr.Dataset({'time': (arm_ds['time'].dims, time, arm_ds['time'].attrs)}) arm_ds['time'] = temp_ds['time'] + temp_ds.close() for att_name in ['units', 'ancillary_variables']: try: del arm_ds['time'].attrs[att_name] @@ -180,8 +186,17 @@ def read_netcdf(filenames, concat_dim='time', return_None=False, # Get file dates and times that were read in to the object filenames.sort() for f in filenames: - file_dates.append(f.split('.')[-3]) - file_times.append(f.split('.')[-2]) + # If Not ARM format, read in first time for infos + if len(f.split('/')[-1].split('.')) == 5: + file_dates.append(f.split('.')[-3]) + file_times.append(f.split('.')[-2]) + else: + if arm_ds['time'].size > 1: + dummy = arm_ds['time'].values[0] + else: + dummy = arm_ds['time'].values + file_dates.append(utils.numpy_to_arm_date(dummy)) + file_times.append(utils.numpy_to_arm_date(dummy, returnTime=True)) # Add attributes arm_ds.attrs['_file_dates'] = file_dates @@ -266,7 +281,7 @@ def create_obj_from_arm_dod(proc, set_dims, version='', fill_value=-9999., """ # Set base url to get DOD information - base_url = 'https://pcm.arm.gov/pcmserver/dods/' + base_url = 'https://pcm.arm.gov/pcm/api/dods/' # Get data from DOD api with urllib.request.urlopen(base_url + proc) as url: @@ -275,7 +290,9 @@ def create_obj_from_arm_dod(proc, set_dims, version='', fill_value=-9999., # Check version numbers and alert if requested version in not available keys = list(data['versions'].keys()) if version not in keys: - print(' '.join(['Version:', version, 'not available or not specified. Using Version:', keys[-1]])) + warnings.warn(' '.join(['Version:', version, + 'not available or not specified. Using Version:', keys[-1]]), + UserWarning) version = keys[-1] # Create empty xarray dataset @@ -351,9 +368,9 @@ def __init__(self, xarray_obj): self._obj = xarray_obj def write_netcdf(self, cleanup_global_atts=True, cleanup_qc_atts=True, - join_char='__', make_copy=True, + join_char='__', make_copy=True, cf_compliant=False, delete_global_attrs=['qc_standards_version', 'qc_method', 'qc_comment'], - FillValue=-9999, **kwargs): + FillValue=-9999, cf_convention='CF-1.8', **kwargs): """ This is a wrapper around Dataset.to_netcdf to clean up the Dataset before writing to disk. Some things are added to global attributes during ACT reading @@ -372,11 +389,17 @@ def write_netcdf(self, cleanup_global_atts=True, cleanup_qc_atts=True, Will use a single space a delimeter between values and join_char to replace white space between words. join_char : str - The character sting to use for replacing white spaces between words. + The character sting to use for replacing white spaces between words when converting + a list of strings to single character string attributes. make_copy : boolean Make a copy before modifying Dataset to write. For large Datasets this may add processing time and memory. If modifying the Dataset is OK try setting to False. + cf_compliant : boolean + Option to output file with additional attributes to make file Climate & Forecast + complient. May require runing .clean.cleanup() method on the object to fix other + issues first. This does the best it can but it may not be truely complient. You + should read the CF documents and try to make complient before writing to file. delete_global_attrs : list Optional global attributes to be deleted. Defaults to some standard QC attributes that are not needed. Can add more or set to None to not @@ -387,6 +410,8 @@ def write_netcdf(self, cleanup_global_atts=True, cleanup_qc_atts=True, so not a perfect fix. Set to None to leave Xarray to do what it wants. Set to a value to be the value used as _FillValue in the file and data array. This should then remove missing_value attribute from the file as well. + cf_convention : str + The Climate and Forecast convention string to add to Conventions attribute. **kwargs : keywords Keywords to pass through to Dataset.to_netcdf() @@ -447,4 +472,96 @@ def write_netcdf(self, cleanup_global_atts=True, cleanup_qc_atts=True, except KeyError: pass + # If requested update global attributes and variables attributes for required + # CF attributes. + if cf_compliant: + # Get variable names and standard name for each variable + var_names = list(write_obj.keys()) + standard_names = [] + for var_name in var_names: + try: + standard_names.append(write_obj[var_name].attrs['standard_name']) + except KeyError: + standard_names.append(None) + + # Check if time varible has axis and standard_name attribute + coord_name = 'time' + try: + write_obj[coord_name].attrs['axis'] + except KeyError: + try: + write_obj[coord_name].attrs['axis'] = 'T' + except KeyError: + pass + + try: + write_obj[coord_name].attrs['standard_name'] + except KeyError: + try: + write_obj[coord_name].attrs['standard_name'] = 'time' + except KeyError: + pass + + # Try to determine type of dataset by coordinate dimention named time + # and other factors + try: + write_obj.attrs['FeatureType'] + except KeyError: + dim_names = list(write_obj.dims) + FeatureType = None + if dim_names == ['time']: + FeatureType = "timeSeries" + elif len(dim_names) == 2 and 'time' in dim_names and 'bound' in dim_names: + FeatureType = "timeSeries" + elif len(dim_names) >= 2 and 'time' in dim_names: + for var_name in var_names: + dims = list(write_obj[var_name].dims) + if len(dims) == 2 and 'time' in dims: + prof_dim = list(set(dims) - set(['time']))[0] + if write_obj[prof_dim].values.size > 2: + FeatureType = "timeSeriesProfile" + break + + if FeatureType is not None: + write_obj.attrs['FeatureType'] = FeatureType + + # Add axis and positive attributes to variables with standard_name + # equal to 'altitude' + alt_variables = [var_names[ii] for ii, sn in enumerate(standard_names) if sn == 'altitude'] + for var_name in alt_variables: + try: + write_obj[var_name].attrs['axis'] + except KeyError: + write_obj[var_name].attrs['axis'] = 'Z' + + try: + write_obj[var_name].attrs['positive'] + except KeyError: + write_obj[var_name].attrs['positive'] = 'up' + + # Check if the Conventions global attribute lists the CF convention + try: + Conventions = write_obj.attrs['Conventions'] + Conventions = Conventions.split() + cf_listed = False + for ii in Conventions: + if ii.startswith('CF-'): + cf_listed = True + break + if not cf_listed: + Conventions.append(cf_convention) + write_obj.attrs['Conventions'] = ' '.join(Conventions) + + except KeyError: + write_obj.attrs['Conventions'] = str(cf_convention) + + # Reorder global attributes to ensure history is last + try: + global_attrs = write_obj.attrs + history = copy.copy(global_attrs['history']) + del global_attrs['history'] + global_attrs['history'] = history + except KeyError: + pass + write_obj.to_netcdf(encoding=encoding, **kwargs) diff --git a/act/io/csvfiles.py b/act/io/csvfiles.py index 792e057513..44327574be 100644 --- a/act/io/csvfiles.py +++ b/act/io/csvfiles.py @@ -8,6 +8,7 @@ # import standard modules import pandas as pd +import pathlib from .armfiles import check_arm_standards @@ -48,6 +49,14 @@ def read_csv(filename, sep=',', engine='python', column_names=None, act.tests.sample_files.EXAMPLE_CSV_WILDCARD) """ + + # Convert to string if filename is a pathlib + if isinstance(filename, pathlib.PurePath): + filename = str(filename) + + if isinstance(filename, list) and isinstance(filename[0], pathlib.PurePath): + filename = [str(ii) for ii in filename] + # Read data using pandas read_csv arm_ds = pd.read_csv(filename, sep=sep, names=column_names, skipfooter=skipfooter, engine=engine, **kwargs) diff --git a/act/io/noaagml.py b/act/io/noaagml.py new file mode 100644 index 0000000000..db980b3fa0 --- /dev/null +++ b/act/io/noaagml.py @@ -0,0 +1,761 @@ +import act +from pathlib import Path +import numpy as np +from datetime import datetime +import xarray as xr +import re + + +def read_gml(filename, datatype=None, **kwargs): + """ + Function to call or guess what reading NOAA GML daga routine to use. It tries to + guess the correct reading function to call based on filename. It mostly + works, but you may want to specify for best results. + + Parameters + ---------- + filename : str or pathlib.Path + Data file full path name. In theory it should work with a list of + filenames but it is not working as well with that as expected. + datatype : str + Data file type that bypasses the guessing from filename format + and goes directly to the reading routine. Options include + [MET, RADIATION, OZONE, CO2, HALO] + **kwargs : keywords + Keywords to pass through to reading routine. + + Returns + ------- + dataset : Xarray.dataset + Standard ARM Xarray dataset with the data cleaned up to have units, + long_name, correct type and some other stuff. + + """ + + if datatype is not None: + if datatype.upper() == 'MET': + return read_gml_met(filename, **kwargs) +# elif datatype.upper() == 'AEROSOL': +# return None + elif datatype.upper() == 'RADIATION': + return read_gml_radiation(filename, **kwargs) + elif datatype.upper() == 'OZONE': + return read_gml_ozone(filename, **kwargs) + elif datatype.upper() == 'CO2': + return read_gml_co2(filename, **kwargs) + elif datatype.upper() == 'HALO': + return read_gml_halo(filename, **kwargs) + else: + raise ValueError("datatype is unknown") + + else: + test_filename = filename + if isinstance(test_filename, (list, tuple)): + test_filename = filename[0] + + test_filename = str(Path(test_filename).name) + + if test_filename.startswith('met_') and test_filename.endswith('.txt'): + return read_gml_met(filename, **kwargs) + + if test_filename.startswith('co2_') and test_filename.endswith('.txt'): + return read_gml_co2(filename, **kwargs) + +# if test_filename.startswith('ESRL-GMD-AEROSOL_v') and test_filename.endswith('.nc'): +# ds_data = xr.open_dataset(str(filename), group='data') +# ds_absorption = xr.open_dataset(str(filename), group='data/light_absorption') +# ds_concentration = xr.open_dataset(str(filename), group='data/particle_concentration') +# ds_scattering = xr.open_dataset(str(filename), group='data/light_scattering') +# return (ds_data, ds_absorption, ds_concentration, ds_scattering) + + result = re.match(r'([a-z]{3})([\d]{5}).dat', test_filename) + if result is not None: + return read_gml_radiation(filename, **kwargs) + + ozone_pattern = [r'[a-z]{3}_[\d]{2}_[\d]{4}_hour.dat', + r'[a-z]{3}_[\d]{4}_all_minute.dat', + r'[a-z]{3}_[\d]{2}_[\d]{4}_5minute.dat', + r'[a-z]{3}_[\d]{2}_[\d]{4}_min.dat', + r'[a-z]{3}_o3_6m_hour_[\d]{2}_[\d]{4}.dat', + r'[a-z]{3}_ozone_houry__[\d]{4}'] + for pattern in ozone_pattern: + result = re.match(pattern, test_filename) + if result is not None: + return read_gml_ozone(filename, **kwargs) + + ozone_pattern = [r'[a-z]{3}_CCl4_Day.dat', + r'[a-z]{3}_CCl4_All.dat', + r'[a-z]{3}_CCl4_MM.dat', + r'[a-z]{3}_MC_MM.dat'] + for pattern in ozone_pattern: + result = re.match(pattern, test_filename) + if result is not None: + return read_gml_halo(filename, **kwargs) + + +def read_gml_halo(filename, **kwargs): + """ + Function to read Halocarbon data from NOAA GML. + + Parameters + ---------- + filename : str or pathlib.Path + Data file full path name. + + Returns + ------- + dataset : Xarray.dataset + Standard ARM Xarray dataset with the data cleaned up to have units, + long_name, correct type and some other stuff. + **kwargs : keywords + Keywords to pass through to ACT read_csv() routine. + + """ + + ds = None + if filename is None: + return ds + + variables = { + 'CCl4catsBRWm': + {'long_name': 'Carbon Tetrachloride (CCl4) daily median', 'units': 'ppt', + '_FillValue': np.nan, '__type': np.float32, '__rename': 'CCl4'}, + 'CCl4catsBRWmsd': + {'long_name': 'Carbon Tetrachloride (CCl4) standard deviation', 'units': 'ppt', + '_FillValue': np.nan, '__type': np.float32, '__rename': 'CCl4_std_dev'}, + 'CCl4catsBRWsd': + {'long_name': 'Carbon Tetrachloride (CCl4) standard deviation', 'units': 'ppt', + '_FillValue': np.nan, '__type': np.float32, '__rename': 'CCl4_std_dev'}, + 'CCl4catsBRWn': + {'long_name': 'Number of samples', 'units': 'count', + '__type': np.int16, '__rename': 'number_of_samples'}, + 'CCl4catsBRWunc': + {'long_name': 'Carbon Tetrachloride (CCl4) uncertainty', 'units': 'ppt', + '_FillValue': np.nan, + '__type': np.float32, '__rename': 'CCl4_uncertainty'}, + 'MCcatsBRWm': + {'long_name': 'Methyl Chloroform (CH3CCl3)', 'units': 'ppt', + '_FillValue': np.nan, + '__type': np.float32, '__rename': 'methyl_chloroform'}, + 'MCcatsBRWunc': + {'long_name': 'Methyl Chloroform (CH3CCl3) uncertainty', 'units': 'ppt', + '_FillValue': np.nan, + '__type': np.float32, '__rename': 'methyl_chloroform_uncertainty'}, + 'MCcatsBRWsd': + {'long_name': 'Methyl Chloroform (CH3CCl3) standard deviation', 'units': 'ppt', + '_FillValue': np.nan, + '__type': np.float32, '__rename': 'methyl_chloroform_std_dev'}, + 'MCcatsBRWmsd': + {'long_name': 'Methyl Chloroform (CH3CCl3) standard deviation', 'units': 'ppt', + '_FillValue': np.nan, + '__type': np.float32, '__rename': 'methyl_chloroform_std_dev'}, + 'MCcatsBRWn': + {'long_name': 'Number of samples', 'units': 'count', + '__type': np.int16, '__rename': 'number_of_samples'}, + 'MCritsBRWm': + {'long_name': 'Methyl Chloroform (CH3CCl3)', 'units': 'ppt', + '_FillValue': np.nan, + '__type': np.float32, '__rename': 'methyl_chloroform'}, + 'MCritsBRWsd': + {'long_name': 'Methyl Chloroform (CH3CCl3) standard deviation', 'units': 'ppt', + '_FillValue': np.nan, + '__type': np.float32, '__rename': 'methyl_chloroform_std_dev'}, + 'MCritsBRWn': + {'long_name': 'Number of samples', 'units': 'count', + '__type': np.int16, '__rename': 'number_of_samples'}, + } + + test_filename = filename + if isinstance(test_filename, (list, tuple)): + test_filename = test_filename[0] + + with open(test_filename, 'r') as fc: + header = 0 + while True: + line = fc.readline().strip() + if not line.startswith('#'): + break + header += 1 + + ds = act.io.csvfiles.read_csv(filename, sep=r'\s+', header=header, + na_values=['Nan', 'NaN', 'nan', 'NAN']) + var_names = list(ds.data_vars) + year_name, month_name, day_name, hour_name, min_name = None, None, None, None, None + for var_name in var_names: + if var_name.endswith('yr'): + year_name = var_name + elif var_name.endswith('mon'): + month_name = var_name + elif var_name.endswith('day'): + day_name = var_name + elif var_name.endswith('hour'): + hour_name = var_name + elif var_name.endswith('min'): + min_name = var_name + + timestamp = np.full(ds[var_names[0]].size, np.nan, dtype='datetime64[s]') + for ii in range(0, len(timestamp)): + if min_name is not None: + ts = datetime(ds[year_name].values[ii], ds[month_name].values[ii], ds[day_name].values[ii], + ds[hour_name].values[ii], ds[min_name].values[ii]) + elif hour_name is not None: + ts = datetime(ds[year_name].values[ii], ds[month_name].values[ii], ds[day_name].values[ii], + ds[hour_name].values[ii]) + elif day_name is not None: + ts = datetime(ds[year_name].values[ii], ds[month_name].values[ii], ds[day_name].values[ii]) + else: + ts = datetime(ds[year_name].values[ii], ds[month_name].values[ii], 1) + + timestamp[ii] = np.datetime64(ts) + + for var_name in [year_name, month_name, day_name, hour_name, min_name]: + try: + del ds[var_name] + except KeyError: + pass + + ds = ds.rename({'index': 'time'}) + ds = ds.assign_coords(time=timestamp) + ds['time'].attrs['long_name'] = 'Time' + + for var_name, value in variables.items(): + if var_name not in var_names: + continue + + for att_name, att_value in value.items(): + if att_name == '__type': + values = ds[var_name].values + values = values.astype(att_value) + ds[var_name].values = values + elif att_name == '__rename': + ds = ds.rename({var_name: att_value}) + else: + ds[var_name].attrs[att_name] = att_value + + return ds + + +def read_gml_co2(filename=None, convert_missing=True, **kwargs): + """ + Function to read carbon dioxide data from NOAA GML. + + Parameters + ---------- + filename : str or pathlib.Path + Data file full path name. + convert_missing : boolean + Option to convert missing values to NaN. If turned off will + set variable attribute to missing value expected. This works well + to preserve the data type best for writing to a netCDF file. + + Returns + ------- + dataset : Xarray.dataset + Standard ARM Xarray dataset with the data cleaned up to have units, + long_name, correct type and some other stuff. + **kwargs : keywords + Keywords to pass through to ACT read_csv() routine. + + """ + ds = None + if filename is None: + return ds + + variables = {'site_code': None, + 'year': None, + 'month': None, + 'day': None, + 'hour': None, + 'minute': None, + 'second': None, + 'time_decimal': None, + 'value': + {'long_name': 'Carbon monoxide in dry air', 'units': 'ppm', + '_FillValue': -999.99, + 'comment': ('Mole fraction reported in units of micromol mol-1 ' + '(10-6 mol per mol of dry air); abbreviated as ppm (parts per million).'), + '__type': np.float32, '__rename': 'co2'}, + 'value_std_dev': + {'long_name': 'Carbon monoxide in dry air', 'units': 'ppm', + '_FillValue': -99.99, + 'comment': ('This is the standard deviation of the reported mean value ' + 'when nvalue is greater than 1. See provider_comment if available.'), + '__type': np.float32, '__rename': 'co2_std_dev'}, + 'nvalue': + {'long_name': 'Number of measurements contributing to reported value', + 'units': '1', '_FillValue': -9, + '__type': np.int16, '__rename': 'number_of_measurements'}, + 'latitude': + {'long_name': 'Latitude at which air sample was collected', + 'units': 'degrees_north', '_FillValue': -999.999, 'standard_name': "latitude", + '__type': np.float32}, + 'longitude': + {'long_name': 'Latitude at which air sample was collected', + 'units': 'degrees_east', '_FillValue': -999.999, 'standard_name': "longitude", + '__type': np.float32}, + 'altitude': + {'long_name': 'Sample altitude', + 'units': 'm', '_FillValue': -999.999, 'standard_name': "altitude", + 'comment': ('Altitude for this dataset is the sum of surface elevation ' + '(masl) and sample intake height (magl)'), + '__type': np.float32}, + 'intake_height': + {'long_name': 'Sample intake height above ground level', + 'units': 'm', '_FillValue': -999.999, + '__type': np.float32}, + } + + test_filename = filename + if isinstance(test_filename, (list, tuple)): + test_filename = test_filename[0] + + with open(test_filename, 'r') as fc: + skiprows = int(fc.readline().strip().split()[-1]) - 1 + + ds = act.io.csvfiles.read_csv(filename, sep=r'\s+', skiprows=skiprows) + + timestamp = np.full(ds['year'].size, np.nan, dtype='datetime64[s]') + for ii in range(0, len(timestamp)): + ts = datetime(ds['year'].values[ii], ds['month'].values[ii], ds['day'].values[ii], + ds['hour'].values[ii], ds['minute'].values[ii], ds['second'].values[ii]) + timestamp[ii] = np.datetime64(ts) + + ds = ds.rename({'index': 'time'}) + ds = ds.assign_coords(time=timestamp) + ds['time'].attrs['long_name'] = 'Time' + + for var_name, value in variables.items(): + if value is None: + del ds[var_name] + else: + for att_name, att_value in value.items(): + if att_name == '__type': + values = ds[var_name].values + values = values.astype(att_value) + ds[var_name].values = values + elif att_name == '__rename': + ds = ds.rename({var_name: att_value}) + else: + ds[var_name].attrs[att_name] = att_value + + if convert_missing: + try: + var_name = variables[var_name]['__rename'] + except KeyError: + pass + + try: + missing_value = ds[var_name].attrs['_FillValue'] + values = ds[var_name].values.astype(float) + values[np.isclose(missing_value, values)] = np.nan + ds[var_name].values = values + ds[var_name].attrs['_FillValue'] = np.nan + except KeyError: + pass + + values = ds['qcflag'].values + bad_index = [] + suspect_index = [] + for ii, value in enumerate(values): + pts = list(value) + if pts[0] != '.': + bad_index.append(ii) + if pts[1] != '.': + suspect_index.append(ii) + + var_name = 'co2' + qc_var_name = ds.qcfilter.create_qc_variable(var_name) + ds.qcfilter.add_test(var_name, index=bad_index, test_assessment='Bad', + test_meaning='Obvious problems during collection or analysis') + ds.qcfilter.add_test(var_name, index=suspect_index, test_assessment='Indeterminate', + test_meaning=('Likely valid but does not meet selection criteria determined by ' + 'the goals of a particular investigation')) + ds[qc_var_name].attrs['comment'] = 'This quality control flag is provided by the contributing PIs' + del ds['qcflag'] + + return ds + + +def read_gml_ozone(filename=None, **kwargs): + """ + Function to read carbon dioxide data from NOAA GML. + + Parameters + ---------- + filename : str or pathlib.Path + Data file full path name. + **kwargs : keywords + Keywords to pass through to ACT read_csv() routine. + + Returns + ------- + dataset : Xarray.dataset + Standard ARM Xarray dataset with the data cleaned up to have units, + long_name, correct type and some other stuff. + """ + + ds = None + if filename is None: + return ds + + test_filename = filename + if isinstance(test_filename, (list, tuple)): + test_filename = test_filename[0] + + with open(test_filename, 'r') as fc: + skiprows = 0 + while True: + line = fc.readline().strip().split() + try: + if len(line) == 6 and line[0] == 'STN': + break + except IndexError: + pass + skiprows += 1 + + ds = act.io.csvfiles.read_csv(filename, sep=r'\s+', skiprows=skiprows) + ds.attrs['station'] = str(ds['STN'].values[0]).lower() + + timestamp = np.full(ds['YEAR'].size, np.nan, dtype='datetime64[s]') + for ii in range(0, len(timestamp)): + ts = datetime(ds['YEAR'].values[ii], ds['MON'].values[ii], ds['DAY'].values[ii], + ds['HR'].values[ii]) + timestamp[ii] = np.datetime64(ts) + + ds = ds.rename({'index': 'time'}) + ds = ds.assign_coords(time=timestamp) + ds['time'].attrs['long_name'] = 'Time' + + for var_name in ['STN', 'YEAR', 'MON', 'DAY', 'HR']: + del ds[var_name] + + var_name = 'ozone' + ds = ds.rename({'O3(PPB)': var_name}) + ds[var_name].attrs['long_name'] = 'Ozone' + ds[var_name].attrs['units'] = 'ppb' + ds[var_name].attrs['_FillValue'] = np.nan + ds[var_name].values = ds[var_name].values.astype(np.float32) + + return ds + + +def read_gml_radiation(filename=None, convert_missing=True, **kwargs): + """ + Function to read radiation data from NOAA GML. + + Parameters + ---------- + filename : str or pathlib.Path + Data file full path name. + convert_missing : boolean + Option to convert missing values to NaN. If turned off will + set variable attribute to missing value expected. This works well + to preserve the data type best for writing to a netCDF file. + **kwargs : keywords + Keywords to pass through to ACT read_csv() routine. + + Returns + ------- + dataset : Xarray.dataset + Standard ARM Xarray dataset with the data cleaned up to have units, + long_name, correct type and some other stuff. + """ + + ds = None + if filename is None: + return ds + + column_names = {'year': None, + 'jday': None, + 'month': None, + 'day': None, + 'hour': None, + 'minute': None, + 'decimal_time': None, + 'solar_zenith_angle': + {'units': 'degree', + 'long_name': 'Solar zenith angle', + '_FillValue': -9999.9, '__type': np.float32}, + 'downwelling_global_solar': + {'units': 'W/m^2', + 'long_name': 'Downwelling global solar', + '_FillValue': -9999.9, '__type': np.float32}, + 'upwelling_global_solar': + {'units': 'W/m^2', + 'long_name': 'Upwelling global solar', + '_FillValue': -9999.9, '__type': np.float32}, + 'direct_normal_solar': + {'units': 'W/m^2', + 'long_name': 'Direct-normal solar', + '_FillValue': -9999.9, '__type': np.float32}, + 'downwelling_diffuse_solar': + {'units': 'W/m^2', + 'long_name': 'Downwelling diffuse solar', + '_FillValue': -9999.9, '__type': np.float32}, + 'downwelling_thermal_infrared': + {'units': 'W/m^2', + 'long_name': 'Downwelling thermal infrared', + '_FillValue': -9999.9, '__type': np.float32}, + 'downwelling_infrared_case_temp': + {'units': 'degK', + 'long_name': 'Downwelling infrared case temp', + '_FillValue': -9999.9, '__type': np.float32}, + 'downwelling_infrared_dome_temp': + {'units': 'degK', + 'long_name': 'downwelling infrared dome temp', + '_FillValue': -9999.9, '__type': np.float32}, + 'upwelling_thermal_infrared': + {'units': 'W/m^2', + 'long_name': 'Upwelling thermal infrared', + '_FillValue': -9999.9, '__type': np.float32}, + 'upwelling_infrared_case_temp': + {'units': 'degK', + 'long_name': 'Upwelling infrared case temp', + '_FillValue': -9999.9, '__type': np.float32}, + 'upwelling_infrared_dome_temp': + {'units': 'degK', + 'long_name': 'Upwelling infrared dome temp', + '_FillValue': -9999.9, '__type': np.float32}, + 'global_UVB': + {'units': 'mW/m^2', + 'long_name': 'global ultraviolet-B', + '_FillValue': -9999.9, '__type': np.float32}, + 'par': + {'units': 'W/m^2', + 'long_name': 'Photosynthetically active radiation', + '_FillValue': -9999.9, '__type': np.float32}, + 'net_solar': + {'units': 'W/m^2', + 'long_name': 'Net solar (downwelling_global_solar - upwelling_global_solar)', + '_FillValue': -9999.9, '__type': np.float32}, + 'net_infrared': + {'units': 'W/m^2', + 'long_name': ('Net infrared (downwelling_thermal_infrared - ' + 'upwelling_thermal_infrared)'), + '_FillValue': -9999.9, '__type': np.float32}, + 'net_radiation': + {'units': 'W/m^2', + 'long_name': 'Net radiation (net_solar + net_infrared)', + '_FillValue': -9999.9, '__type': np.float32}, + 'air_temperature_10m': + {'units': 'degC', + 'long_name': '10-meter air temperature', + '_FillValue': -9999.9, '__type': np.float32}, + 'relative_humidity': + {'units': '%', + 'long_name': 'Relative humidity', + '_FillValue': -9999.9, '__type': np.float32}, + 'wind_speed': + {'units': 'm/s', + 'long_name': 'Wind speed', + '_FillValue': -9999.9, '__type': np.float32}, + 'wind_direction': + {'units': 'degree', + 'long_name': 'Wind direction (clockwise from north)', + '_FillValue': -9999.9, '__type': np.float32}, + 'station_pressure': + {'units': 'millibar', + 'long_name': 'Station atmospheric pressure', + '_FillValue': -9999.9, '__type': np.float32}, + } + + names = list(column_names.keys()) + skip_vars = ['year', 'jday', 'month', 'day', 'hour', 'minute', 'decimal_time', 'solar_zenith_angle'] + num = 1 + for ii, name in enumerate(column_names.keys()): + if name in skip_vars: + continue + names.insert(ii + num, 'qc_' + name) + num += 1 + + ds = act.io.csvfiles.read_csv(filename, sep=r'\s+', header=0, skiprows=2, column_names=names) + + if ds is not None: + with open(filename, 'r') as fc: + lat = None + lon = None + alt = None + alt_unit = None + station = None + for ii in [0, 1]: + line = fc.readline().strip().split() + if len(line) == 1: + station = line[0] + else: + lat = np.array(line[0], dtype=np.float32) + lon = np.array(line[1], dtype=np.float32) + alt = np.array(line[2], dtype=np.float32) + alt_unit = str(line[3]) + + ds['lat'] = xr.DataArray(lat, attrs={'long_name': 'Latitude', 'units': 'degree_north', + 'standard_name': 'latitude'}) + ds['lon'] = xr.DataArray(lon, attrs={'long_name': 'Longitude', 'units': 'degree_east', + 'standard_name': 'longitude'}) + ds['alt'] = xr.DataArray(alt, attrs={'long_name': 'Latitude', 'units': alt_unit, + 'standard_name': 'altitude'}) + ds.attrs['location'] = station + + timestamp = np.full(ds['year'].size, np.nan, dtype='datetime64[s]') + for ii in range(0, len(timestamp)): + ts = datetime(ds['year'].values[ii], ds['month'].values[ii], ds['day'].values[ii], + ds['hour'].values[ii], ds['minute'].values[ii]) + timestamp[ii] = np.datetime64(ts) + + ds = ds.rename({'index': 'time'}) + ds = ds.assign_coords(time=timestamp) + ds['time'].attrs['long_name'] = 'Time' + for var_name, value in column_names.items(): + if value is None: + ds[var_name] + else: + for att_name, att_value in value.items(): + if att_name == '__type': + values = ds[var_name].values + values = values.astype(att_value) + ds[var_name].values = values + else: + ds[var_name].attrs[att_name] = att_value + + if convert_missing: + try: + missing_value = ds[var_name].attrs['_FillValue'] + values = ds[var_name].values.astype(float) + index = np.isclose(values, missing_value) + values[index] = np.nan + ds[var_name].values = values + ds[var_name].attrs['_FillValue'] = np.nan + except KeyError: + pass + + for var_name in ds.data_vars: + if not var_name.startswith('qc_'): + continue + data_var_name = var_name.replace('qc_', '', 1) + attrs = {'long_name': f"Quality control variable for: {ds[data_var_name].attrs['long_name']}", + 'units': '1', 'standard_name': 'quality_flag', + 'flag_values': [0, 1, 2], + 'flag_meanings': ['Not failing any tests', 'Knowingly bad value', + 'Should be used with scrutiny'], + 'flag_assessments': ['Good', 'Bad', 'Indeterminate']} + ds[var_name].attrs = attrs + ds[data_var_name].attrs['ancillary_variables'] = var_name + + return ds + + +def read_gml_met(filename=None, convert_missing=True, **kwargs): + """ + Function to read meteorlogical data from NOAA GML. + + Parameters + ---------- + filename : str or pathlib.Path + Data file full path name. + convert_missing : boolean + Option to convert missing values to NaN. If turned off will + set variable attribute to missing value expected. This works well + to preserve the data type best for writing to a netCDF file. + **kwargs : keywords + Keywords to pass through to ACT read_csv() routine. + + Returns + ------- + dataset : Xarray.dataset + Standard ARM Xarray dataset with the data cleaned up to have units, + long_name, correct type and some other stuff. + """ + + ds = None + if filename is None: + return ds + + column_names = {'station': None, + 'year': None, + 'month': None, + 'day': None, + 'hour': None, + 'minute': None, + 'wind_direction': + {'units': 'degree', + 'long_name': 'Average wind direction from which the wind is blowing', + '_FillValue': -999, '__type': np.int16}, + 'wind_speed': + {'units': 'm/s', 'long_name': 'Average wind speed', '_FillValue': -999.9, + '__type': np.float32}, + 'wind_steadiness_factor': + {'units': '1', 'long_name': '100 times the ratio of the vector wind speed to the ' + 'average wind speed for the hour', '_FillValue': -9, '__type': np.int16}, + 'barometric_pressure': + {'units': 'hPa', 'long_name': 'Station barometric pressure', '_FillValue': -999.90, + '__type': np.float32}, + 'temperature_2m': + {'units': 'degC', 'long_name': 'Temperature at 2 meters above ground level', + '_FillValue': -999.9, '__type': np.float32}, + 'temperature_10m': + {'units': 'degC', 'long_name': 'Temperature at 10 meters above ground level', + '_FillValue': -999.9, '__type': np.float32}, + 'temperature_tower_top': + {'units': 'degC', 'long_name': 'Temperature at top of instrument tower', + '_FillValue': -999.9, '__type': np.float32}, + 'realitive_humidity': + {'units': 'percent', 'long_name': 'Relative humidity', '_FillValue': -99, + '__type': np.int16}, + 'preciptation_intensity': + {'units': 'mm/hour', 'long_name': 'Amount of precipitation per hour', + '_FillValue': -99, '__type': np.int16, + 'comment': ('The precipitation amount is measured with an unheated ' + 'tipping bucket rain gauge.')} + } + + minutes = True + test_filename = filename + if isinstance(test_filename, (list, tuple)): + test_filename = test_filename[0] + + if '_hour_' in Path(test_filename).name: + minutes = False + del column_names['minute'] + + ds = act.io.csvfiles.read_csv(filename, sep=r'\s+', header=0, column_names=column_names.keys()) + + if ds is not None: + + timestamp = np.full(ds['year'].size, np.nan, dtype='datetime64[s]') + for ii in range(0, len(timestamp)): + if minutes: + ts = datetime(ds['year'].values[ii], ds['month'].values[ii], ds['day'].values[ii], + ds['hour'].values[ii], ds['minute'].values[ii]) + else: + ts = datetime(ds['year'].values[ii], ds['month'].values[ii], ds['day'].values[ii], + ds['hour'].values[ii]) + + timestamp[ii] = np.datetime64(ts) + + ds = ds.rename({'index': 'time'}) + ds = ds.assign_coords(time=timestamp) + ds['time'].attrs['long_name'] = 'Time' + for var_name, value in column_names.items(): + + if value is None: + del ds[var_name] + else: + for att_name, att_value in value.items(): + if att_name == '__type': + values = ds[var_name].values + values = values.astype(att_value) + ds[var_name].values = values + else: + ds[var_name].attrs[att_name] = att_value + + if convert_missing: + try: + missing_value = ds[var_name].attrs['_FillValue'] + values = ds[var_name].values.astype(float) + index = np.isclose(values, missing_value) + values[index] = np.nan + ds[var_name].values = values + ds[var_name].attrs['_FillValue'] = np.nan + except KeyError: + pass + + return ds diff --git a/act/plotting/ContourDisplay.py b/act/plotting/ContourDisplay.py index fd67964749..6392b729da 100644 --- a/act/plotting/ContourDisplay.py +++ b/act/plotting/ContourDisplay.py @@ -24,8 +24,9 @@ def __init__(self, obj, subplot_shape=(1,), ds_name=None, **kwargs): super().__init__(obj, subplot_shape, ds_name, **kwargs) def create_contour(self, fields=None, time=None, function='cubic', - subplot_index=(0,), grid_delta=(0.01, 0.01), - grid_buffer=0.1, **kwargs): + subplot_index=(0,), contour='contourf', + grid_delta=(0.01, 0.01), grid_buffer=0.1, + twod_dim_value=None, **kwargs): """ Extracts, grids, and creates a contour plot. If subplots have not been added yet, an axis will be created assuming that there is only going @@ -42,10 +43,17 @@ def create_contour(self, fields=None, time=None, function='cubic', See scipy.interpolate.Rbf for additional options. subplot_index : 1 or 2D tuple, list, or array The index of the subplot to set the x range of. + contour : str + Whether to use contour or contourf as plotting function. + Default is contourf grid_delta : 1D tuple, list, or array x and y deltas for creating grid. grid_buffer : float Buffer to apply to grid. + twod_dim_value : float + If the field is 2D, which dimension value to pull. + I.e. if dim is depths of [5, 10, 50, 100] specifying 50 + would index the data at 50 **kwargs : keyword arguments The keyword arguments for :func:`plt.contour` @@ -62,10 +70,34 @@ def create_contour(self, fields=None, time=None, function='cubic', z = [] for ds in self._arm: obj = self._arm[ds] + if ds not in fields: + continue field = fields[ds] - x.append(obj[field[0]].sel(time=time).values.tolist()) - y.append(obj[field[1]].sel(time=time).values.tolist()) - z.append(obj[field[2]].sel(time=time).values.tolist()) + if obj[field[2]].sel(time=time).values.size > 1: + dim_values = obj[obj[field[2]].dims[1]].values + if twod_dim_value is None: + dim_index = 0 + else: + dim_index = np.where((dim_values == twod_dim_value)) + if dim_index[0].size == 0: + continue + if np.isnan(obj[field[2]].sel(time=time).values[dim_index]): + continue + z.append(obj[field[2]].sel(time=time).values[dim_index].tolist()) + else: + if np.isnan(obj[field[2]].sel(time=time).values): + continue + z.append(obj[field[2]].sel(time=time).values.tolist()) + + if obj[field[0]].values.size > 1: + x.append(obj[field[0]].sel(time=time).values.tolist()) + else: + x.append(obj[field[0]].values.tolist()) + + if obj[field[1]].values.size > 1: + y.append(obj[field[1]].sel(time=time).values.tolist()) + else: + y.append(obj[field[1]].values.tolist()) # Create a meshgrid for gridding onto xs = np.arange(np.min(x) - grid_buffer, np.max(x) + grid_buffer, grid_delta[0]) @@ -77,7 +109,13 @@ def create_contour(self, fields=None, time=None, function='cubic', zi = rbf(xi, yi) # Create contour plot - self.contourf(xi, yi, zi, subplot_index=subplot_index, **kwargs) + if contour == 'contour': + self.contour(xi, yi, zi, subplot_index=subplot_index, **kwargs) + elif 'contourf': + self.contourf(xi, yi, zi, subplot_index=subplot_index, **kwargs) + else: + raise ValueError( + "Invalid contour plot type. Please choose either 'contourf' or 'contour'") return self.axes[subplot_index] @@ -181,8 +219,15 @@ def plot_vectors_from_spd_dir(self, fields, time=None, subplot_index=(0,), for ds in self._arm: obj = self._arm[ds] field = fields[ds] - x.append(obj[field[0]].sel(time=time).values.tolist()) - y.append(obj[field[1]].sel(time=time).values.tolist()) + if obj[field[0]].values.size > 1: + x.append(obj[field[0]].sel(time=time).values.tolist()) + else: + x.append(obj[field[0]].values.tolist()) + + if obj[field[1]].values.size > 1: + y.append(obj[field[1]].sel(time=time).values.tolist()) + else: + y.append(obj[field[1]].values.tolist()) wspd.append(obj[field[2]].sel(time=time).values.tolist()) wdir.append(obj[field[3]].sel(time=time).values.tolist()) @@ -276,9 +321,15 @@ def plot_station(self, fields, time=None, subplot_index=(0,), field = fields[ds] for i, f in enumerate(field): if i == 0: - x = obj[f].sel(time=time).values.tolist() + if obj[f].values.size > 1: + x = obj[f].sel(time=time).values.tolist() + else: + x = obj[f].values.tolist() elif i == 1: - y = obj[f].sel(time=time).values.tolist() + if obj[f].values.size > 1: + y = obj[f].sel(time=time).values.tolist() + else: + y = obj[f].values.tolist() self.axes[subplot_index].plot(x, y, '*', **kwargs) else: data = obj[f].sel(time=time).values.tolist() diff --git a/act/plotting/HistogramDisplay.py b/act/plotting/HistogramDisplay.py index a22d453d92..ff374936a8 100644 --- a/act/plotting/HistogramDisplay.py +++ b/act/plotting/HistogramDisplay.py @@ -176,6 +176,10 @@ def plot_stacked_bar_graph(self, field, dsname=None, bins=None, label=(str(y_bins[i]) + " to " + str(y_bins[i + 1])), **kwargs) self.axes[subplot_index].legend() else: + if bins is None: + bmin = np.nanmin(xdata) + bmax = np.nanmax(xdata) + bins = np.arange(bmin, bmax, (bmax - bmin) / 10.) my_hist, bins = np.histogram( xdata.values.flatten(), bins=bins, density=density) x_inds = (bins[:-1] + bins[1:]) / 2.0 diff --git a/act/plotting/SkewTDisplay.py b/act/plotting/SkewTDisplay.py index 2416628e8e..4c27636a45 100644 --- a/act/plotting/SkewTDisplay.py +++ b/act/plotting/SkewTDisplay.py @@ -72,7 +72,7 @@ def __init__(self, obj, subplot_shape=(1,), ds_name=None, **kwargs): subplot_kw=dict(projection='skewx'), **new_kwargs) # Make a SkewT object for each subplot - self.add_subplots(subplot_shape) + self.add_subplots(subplot_shape, **kwargs) def add_subplots(self, subplot_shape=(1,), **kwargs): """ @@ -110,8 +110,8 @@ def add_subplots(self, subplot_shape=(1,), **kwargs): subplot_tuple = (subplot_shape[0], subplot_shape[1], i * subplot_shape[1] + j + 1) - self.SkewT[i] = SkewT(fig=self.fig, subplot=subplot_tuple) - self.axes[i] = self.SkewT[i].ax + self.SkewT[i, j] = SkewT(fig=self.fig, subplot=subplot_tuple) + self.axes[i, j] = self.SkewT[i, j].ax else: raise ValueError("Subplot shape must be 1 or 2D!") @@ -130,10 +130,11 @@ def set_xrng(self, xrng, subplot_index=(0,)): if self.axes is None: raise RuntimeError("set_xrng requires the plot to be displayed.") - if not hasattr(self, 'xrng') and len(self.axes.shape) == 2: - self.xrng = np.zeros((self.axes.shape[0], self.axes.shape[1], 2)) - elif not hasattr(self, 'xrng') and len(self.axes.shape) == 1: - self.xrng = np.zeros((self.axes.shape[0], 2)) + if not hasattr(self, 'xrng') or np.all(self.xrng == 0): + if len(self.axes.shape) == 2: + self.xrng = np.zeros((self.axes.shape[0], self.axes.shape[1], 2)) + else: + self.xrng = np.zeros((self.axes.shape[0], 2)) self.axes[subplot_index].set_xlim(xrng) self.xrng[subplot_index, :] = np.array(xrng) @@ -153,6 +154,12 @@ def set_yrng(self, yrng, subplot_index=(0,)): if self.axes is None: raise RuntimeError("set_yrng requires the plot to be displayed.") + if not hasattr(self, 'yrng') or np.all(self.yrng == 0): + if len(self.axes.shape) == 2: + self.yrng = np.zeros((self.axes.shape[0], self.axes.shape[1], 2)) + else: + self.yrng = np.zeros((self.axes.shape[0], 2)) + if not hasattr(self, 'yrng') and len(self.axes.shape) == 2: self.yrng = np.zeros((self.axes.shape[0], self.axes.shape[1], 2)) elif not hasattr(self, 'yrng') and len(self.axes.shape) == 1: @@ -315,6 +322,8 @@ def plot_from_u_and_v(self, u_field, v_field, p_field, u_red[i] = u[index].magnitude * getattr(units, u_units) v_red[i] = v[index].magnitude * getattr(units, v_units) + u_red = u_red.magnitude + v_red = v_red.magnitude self.SkewT[subplot_index].plot(p, T, 'r', **plot_kwargs) self.SkewT[subplot_index].plot(p, Td, 'g', **plot_kwargs) self.SkewT[subplot_index].plot_barbs( @@ -347,17 +356,12 @@ def plot_from_u_and_v(self, u_field, v_field, p_field, self.axes[subplot_index].set_title(set_title) # Set Y Limit - if hasattr(self, 'yrng'): - # Make sure that the yrng is not just the default - if not np.all(self.yrng[subplot_index] == 0): - self.set_yrng(self.yrng[subplot_index], subplot_index) - else: - our_data = p.magnitude - if np.isfinite(our_data).any(): - yrng = [np.nanmax(our_data), np.nanmin(our_data)] - else: - yrng = [1000., 100.] - self.set_yrng(yrng, subplot_index) + our_data = p.magnitude + if np.isfinite(our_data).any(): + yrng = [np.nanmax(our_data), np.nanmin(our_data)] + else: + yrng = [1000., 100.] + self.set_yrng(yrng, subplot_index) # Set X Limit xrng = [np.nanmin(T.magnitude) - 10., np.nanmax(T.magnitude) + 10.] diff --git a/act/plotting/TimeSeriesDisplay.py b/act/plotting/TimeSeriesDisplay.py index 4cb91e8edb..f1d9b82643 100644 --- a/act/plotting/TimeSeriesDisplay.py +++ b/act/plotting/TimeSeriesDisplay.py @@ -11,12 +11,6 @@ import pandas as pd import datetime as dt import warnings -import astral -try: - from astral.sun import sun, elevation, noon - ASTRAL = True -except ImportError: - ASTRAL = False from re import search as re_search from matplotlib import colors as mplcolors @@ -29,6 +23,7 @@ from ..utils.datetime_utils import reduce_time_ranges, determine_time_delta from ..qc.qcfilter import parse_bit from ..utils import data_utils +from ..utils.geo_utils import get_sunrise_sunset_noon from copy import deepcopy from scipy.interpolate import NearestNDInterpolator @@ -144,7 +139,7 @@ def day_night_background(self, dsname=None, subplot_index=(0, )): return try: - if self._arm[dsname].lat.data.size > 1: + if self._arm[dsname][lat_name].data.size > 1: # Look for non-NaN values to use for locaiton. If not found use first value. lat = self._arm[dsname][lat_name].values index = np.where(np.isfinite(lat))[0] @@ -191,58 +186,22 @@ def day_night_background(self, dsname=None, subplot_index=(0, )): rect = ax.patch rect.set_facecolor('0.85') - # Set the the number of degrees the sun must be below the horizon - # for the dawn/dusk calculation. Need to do this so when the calculation - # sends an error it is not going to be an inacurate switch to setting - # the full day. - if ASTRAL: - astral.solar_depression = 0 - else: - a = astral.Astral() - a.solar_depression = 0 - + # Get date ranges to plot + plot_dates = [] for f in all_dates: - # Loop over previous, current and following days to cover all overlaps - # due to local vs UTC times. for ii in [-1, 0, 1]: - new_time = f + dt.timedelta(days=ii) - try: - if ASTRAL: - obs = astral.Observer(latitude=lat, longitude=lon) - s = sun(obs, new_time) - else: - s = a.sun_utc(new_time, lat, lon) + plot_dates.append(f + dt.timedelta(days=ii)) - # add yellow background for specified time period - ax.axvspan( - s['sunrise'], s['sunset'], facecolor='#FFFFCC', zorder=0) + # Get sunrise, sunset and noon times + sunrise, sunset, noon = get_sunrise_sunset_noon(lat, lon, plot_dates) - # add local solar noon line - ax.axvline(x=s['noon'], linestyle='--', color='y', zorder=1) + # Plot daylight + for ii in range(0, len(sunrise)): + ax.axvspan(sunrise[ii], sunset[ii], facecolor='#FFFFCC', zorder=0) - except ValueError: - # Error for all day and all night is the same. Check to see - # if sun is above horizon at solar noon. If so plot. - if ASTRAL: - elev = elevation(obs, new_time) - else: - elev = a.solar_elevation(new_time, lat, lon) - if elev > 0: - # Make whole background yellow for when sun does not reach - # horizon. Use in high latitude locations. - ax.axvspan(dt.datetime(f.year, f.month, f.day, hour=0, - minute=0, second=0), - dt.datetime(f.year, f.month, f.day, hour=23, - minute=59, second=59), - facecolor='#FFFFCC') - - # add local solar noon line - if ASTRAL: - s_noon = noon(obs, f) - else: - s_noon = a.solar_noon_utc(f, lon) - - ax.axvline(x=s_noon, linestyle='--', color='y') + # Plot noon line + for ii in noon: + ax.axvline(x=ii, linestyle='--', color='y', zorder=1) def set_xrng(self, xrng, subplot_index=(0, )): """ @@ -748,6 +707,7 @@ def plot_barbs_from_u_v(self, u_field, v_field, pres_field=None, barb_step_x = round(num_x / num_barbs_x) if barb_step_x == 0: barb_step_x = 1 + if len(dim) > 1 and pres_field is None: if use_var_for_y is None: ydata = self._arm[dsname][dim[1]] @@ -756,12 +716,16 @@ def plot_barbs_from_u_v(self, u_field, v_field, pres_field=None, ydata_dim1 = self._arm[dsname][dim[1]] if np.shape(ydata) != np.shape(ydata_dim1): ydata = ydata_dim1 - ytitle = ''.join(['(', ydata.attrs['units'], ')']) - units = ytitle + if 'units' in ydata.attrs: + units = ydata.attrs['units'] + else: + units = '' + ytitle = ''.join(['(', units, ')']) num_y = ydata.shape[0] barb_step_y = round(num_y / num_barbs_y) if barb_step_y == 0: barb_step_y = 1 + xdata, ydata = np.meshgrid(xdata, ydata, indexing='ij') elif pres_field is not None: # What we will do here is do a nearest-neighbor interpolation @@ -783,7 +747,11 @@ def plot_barbs_from_u_v(self, u_field, v_field, pres_field=None, xdata, ydata = np.meshgrid(x_times, y_levels, indexing='ij') u = u_interp(xdata, ydata) v = v_interp(xdata, ydata) - ytitle = ''.join(['(', pres.attrs['units'], ')']) + if 'units' in pres.attrs: + units = pres.attrs['units'] + else: + units = '' + ytitle = ''.join(['(', units, ')']) else: ydata = None @@ -823,7 +791,6 @@ def plot_barbs_from_u_v(self, u_field, v_field, pres_field=None, map_color = np.sqrt(np.power(u[::barb_step_x, ::barb_step_y], 2) + np.power(v[::barb_step_x, ::barb_step_y], 2)) map_color[np.isnan(map_color)] = 0 - ax = self.axes[subplot_index].barbs( xdata[::barb_step_x, ::barb_step_y], ydata[::barb_step_x, ::barb_step_y], @@ -981,7 +948,7 @@ def plot_time_height_xsection_from_1d_data( self.fig.add_axes(self.axes[0]) mesh = self.axes[subplot_index].pcolormesh( - tdata, ydata, data, **kwargs) + x_times, y_levels, np.transpose(data), **kwargs) if day_night_background is True: self.day_night_background(subplot_index=subplot_index, dsname=dsname) @@ -1123,7 +1090,7 @@ def time_height_scatter( def qc_flag_block_plot( self, data_field=None, dsname=None, - subplot_index=(0, ), time_rng=None, assesment_color=None, + subplot_index=(0, ), time_rng=None, assessment_color=None, edgecolor='face', **kwargs): """ Create a time series plot of embedded quality control values @@ -1143,7 +1110,7 @@ def qc_flag_block_plot( The index of the subplot to set the x range of. time_rng : tuple or list List or tuple with (min, max) values to set the x-axis range limits. - assesment_color : dict + assessment_color : dict Dictionary lookup to override default assessment to color. Make sure assessment work is correctly set with case syntax. **kwargs : keyword arguments @@ -1156,10 +1123,11 @@ def qc_flag_block_plot( 'Indeterminate': 'orange', 'Suspect': 'orange', 'Missing': 'darkgray', - 'Not Failing': 'green'} + 'Not Failing': 'green', + 'Acceptable': 'green'} - if assesment_color is not None: - for asses, color in assesment_color.items(): + if assessment_color is not None: + for asses, color in assessment_color.items(): color_lookup[asses] = color if asses == 'Incorrect': color_lookup['Bad'] = color @@ -1232,6 +1200,8 @@ def qc_flag_block_plot( tick_names.append('Not Failing') for ii, assess in enumerate(cur_assessments): + if assess not in color_lookup: + color_lookup[assess] = list(mplcolors.CSS4_COLORS.keys())[ii] ii += 1 assess_data = self._arm[dsname].qcfilter.get_masked_data(data_field, rm_assessments=assess) @@ -1301,6 +1271,8 @@ def qc_flag_block_plot( test_nums = [] for ii, assess in enumerate(flag_assessments): + if assess not in color_lookup: + color_lookup[assess] = list(mplcolors.CSS4_COLORS.keys())[ii] # Plot green data first. ax.broken_barh(barh_list_green, (ii, ii + 1), facecolors=color_lookup['Not Failing'], edgecolor=edgecolor, **kwargs) diff --git a/act/plotting/WindRoseDisplay.py b/act/plotting/WindRoseDisplay.py index cedf42a8d0..5ea0c76e88 100644 --- a/act/plotting/WindRoseDisplay.py +++ b/act/plotting/WindRoseDisplay.py @@ -52,12 +52,13 @@ def set_thetarng(self, trng=(0., 360.), subplot_index=(0,)): """ if self.axes is not None: - self.axes[subplot_index].set_thetamin(np.deg2rad(trng[0])) - self.axes[subplot_index].set_thetamax(np.deg2rad(trng[1])) + self.axes[subplot_index].set_thetamin(trng[0]) + self.axes[subplot_index].set_thetamax(trng[1]) self.trng = trng else: raise RuntimeError(("Axes must be initialized before" + " changing limits!")) + print(self.trng) def set_rrng(self, rrng, subplot_index=(0,)): """ @@ -72,8 +73,8 @@ def set_rrng(self, rrng, subplot_index=(0,)): """ if self.axes is not None: - self.axes[subplot_index].set_thetamin(rrng[0]) - self.axes[subplot_index].set_thetamax(rrng[1]) + self.axes[subplot_index].set_rmin(rrng[0]) + self.axes[subplot_index].set_rmax(rrng[1]) self.rrng = rrng else: raise RuntimeError(("Axes must be initialized before" + @@ -82,6 +83,7 @@ def set_rrng(self, rrng, subplot_index=(0,)): def plot(self, dir_field, spd_field, dsname=None, subplot_index=(0,), cmap=None, set_title=None, num_dirs=20, spd_bins=None, tick_interval=3, legend_loc=0, legend_bbox=None, legend_title=None, + calm_threshold=1., **kwargs): """ Makes the wind rose plot from the given dataset. @@ -113,6 +115,8 @@ def plot(self, dir_field, spd_field, dsname=None, subplot_index=(0,), Legend bounding box coordinates legend_title : string Legend title + calm_threshold : float + Winds below this threshold are considered to be calm. **kwargs : keyword arguments Additional keyword arguments will be passed into :func:plt.bar @@ -199,6 +203,11 @@ def plot(self, dir_field, spd_field, dsname=None, subplot_index=(0,), self.axes[subplot_index].set_theta_zero_location("N") self.axes[subplot_index].set_theta_direction(-1) + # Add an annulus with text stating % of time calm + pct_calm = np.sum(spd_data <= calm_threshold) / len(spd_data) * 100 + self.axes[subplot_index].set_rorigin(-2.5) + self.axes[subplot_index].annotate("%3.2f%%\n calm" % pct_calm, xy=(0, -2.5), ha='center', va='center') + # Set the ticks to be nice numbers tick_max = tick_interval * round( np.nanmax(np.cumsum(wind_hist, axis=1)) / tick_interval) diff --git a/act/plotting/plot.py b/act/plotting/plot.py index edf1cdf198..b2ca7e934f 100644 --- a/act/plotting/plot.py +++ b/act/plotting/plot.py @@ -148,6 +148,8 @@ def add_subplots(self, subplot_shape=(1,), subplot_kw=None, **kwargs) self.xrng = np.zeros((subplot_shape[0], subplot_shape[1], 2)) self.yrng = np.zeros((subplot_shape[0], subplot_shape[1], 2)) + if subplot_shape[0] == 1: + ax = ax.reshape(1, subplot_shape[1]) elif len(subplot_shape) == 1: fig, ax = plt.subplots( subplot_shape[0], 1, subplot_kw=subplot_kw, **kwargs) diff --git a/act/qc/__init__.py b/act/qc/__init__.py index 25543e2a82..fe520a7581 100644 --- a/act/qc/__init__.py +++ b/act/qc/__init__.py @@ -34,6 +34,7 @@ qcfilter.QCFilter.check_for_ancillary_qc qcfilter.QCFilter.compare_time_series_trends qcfilter.QCFilter.create_qc_variable + qcfilter.QCFilter.datafilter qcfilter.QCFilter.get_qc_test_mask qcfilter.QCFilter.get_masked_data qcfilter.QCFilter.remove_test diff --git a/act/qc/arm.py b/act/qc/arm.py index ead1400ccd..7933200f0d 100644 --- a/act/qc/arm.py +++ b/act/qc/arm.py @@ -106,8 +106,8 @@ def add_dqr_to_qc(obj, variable=None, assessment='incorrect,suspect', if include is not None and line[0] not in include: continue - line[1] = dt.datetime.fromtimestamp(int(line[1])) - line[2] = dt.datetime.fromtimestamp(int(line[2])) + line[1] = dt.datetime.utcfromtimestamp(int(line[1])) + line[2] = dt.datetime.utcfromtimestamp(int(line[2])) ind = np.where((time >= np.datetime64(line[1])) & (time <= np.datetime64(line[2]))) if len(ind[0]) == 0: continue diff --git a/act/qc/clean.py b/act/qc/clean.py index b1c82ed76a..3a5451da1d 100644 --- a/act/qc/clean.py +++ b/act/qc/clean.py @@ -378,7 +378,7 @@ def get_attr_info(self, variable=None, flag=False): if len(flag_masks) > 0 or len(description_bit_num) > 0: return_dict = dict() return_dict['flag_meanings'] = list(np.array(flag_meanings, - dtype=np.str)) + dtype=str)) if len(flag_masks) > 0 and max(flag_masks) > 2**32 - 1: flag_mask_dtype = np.int64 else: @@ -396,11 +396,11 @@ def get_attr_info(self, variable=None, flag=False): dtype=flag_mask_dtype)) return_dict['flag_assessments'] = list(np.array(flag_assessments, - dtype=np.str)) + dtype=str)) return_dict['flag_tests'] = list(np.array(description_bit_num, dtype=dtype)) return_dict['flag_comments'] = list(np.array(flag_comments, - dtype=np.str)) + dtype=str)) return_dict['arm_attributes'] = arm_attributes else: @@ -604,7 +604,10 @@ def clean_arm_qc(self, # Remove replaced attributes if qc_attributes is not None: arm_attributes = qc_attributes['arm_attributes'] - arm_attributes.extend(['description']) # , 'flag_method']) + if 'description' not in arm_attributes: + arm_attributes.append('description') + if 'flag_method' not in arm_attributes: + arm_attributes.append('flag_method') for attr in arm_attributes: try: del self._obj[qc_var].attrs[attr] diff --git a/act/qc/qcfilter.py b/act/qc/qcfilter.py index 6100faff59..2f175b6d65 100644 --- a/act/qc/qcfilter.py +++ b/act/qc/qcfilter.py @@ -468,7 +468,7 @@ def unset_test(self, var_name, index=None, test_number=None, def available_bit(self, qc_var_name, recycle=False): """ - Method to determine next availble bit or flag to use with a QC test. + Method to determine next available bit or flag to use with a QC test. Will check for flag_masks first and if not found will check for flag_values. This will drive how the next value is chosen. @@ -658,7 +658,6 @@ def get_masked_data(self, var_name, rm_assessments=None, fill_value=1e+20, dtype=float32) """ - qc_var_name = self._obj.qcfilter.check_for_ancillary_qc( var_name, add_if_missing=False) @@ -739,11 +738,11 @@ def datafilter(self, variables=None, rm_assessments=None, rm_tests=None, """ Method to apply quality control variables to data variables by changing the data values in the dataset using quality control variables. - The data variable is changed to to a numpy masked array with failing data - masked or, if requested, to numpy array with failing data set to NaN. - This can be used to update the data variable in the xarray dataset for use with - xarray methods to perform analysis on the data since those methods don't - read the quality control variables. + The data variable is changed to to a numpy masked array with failing + data masked or, if requested, to numpy array with failing data set to + NaN. This can be used to update the data variable in the xarray + dataset for use with xarray methods to perform analysis on the data + since those methods don't read the quality control variables. Parameters ---------- @@ -768,10 +767,10 @@ def datafilter(self, variables=None, rm_assessments=None, rm_tests=None, Opttion to delete quality control variable after processing. Since the data values can not be determined after they are set to NaN and xarray method processing would also process the quality control - variables, the default is to remove the quality control data variables. - If numpy masked arrays are used the data are not lost but would need - to be extracted and set to DataArray to return the dataset back to orginal state. - + variables, the default is to remove the quality control data + variables. If numpy masked arrays are used the data are not lost + but would need to be extracted and set to DataArray to return the + dataset back to original state. Examples -------- @@ -794,7 +793,6 @@ def datafilter(self, variables=None, rm_assessments=None, rm_tests=None, All data: 98.86097717285156, Bad Removed: 99.15148162841797 """ - if variables is not None and isinstance(variables, str): variables = [variables] diff --git a/act/qc/qctests.py b/act/qc/qctests.py index a141839eab..f76baeaa06 100644 --- a/act/qc/qctests.py +++ b/act/qc/qctests.py @@ -33,8 +33,7 @@ class declaration. """ def __init__(self, obj, **kwargs): - - print('test') + self._obj = obj def add_missing_value_test(self, var_name, missing_value=None, missing_value_att_name='missing_value', diff --git a/act/qc/radiometer_tests.py b/act/qc/radiometer_tests.py index a0d6ae4f6d..1833aa72d8 100644 --- a/act/qc/radiometer_tests.py +++ b/act/qc/radiometer_tests.py @@ -11,15 +11,10 @@ import pandas as pd import datetime import dask -import astral import warnings -try: - from astral.sun import sun - ASTRAL = True -except ImportError: - ASTRAL = False from act.utils.datetime_utils import determine_time_delta +from act.utils.geo_utils import get_sunrise_sunset_noon, is_sun_visible def fft_shading_test(obj, variable='diffuse_hemisp_narrowband_filter4', @@ -97,6 +92,7 @@ def fft_shading_test(obj, variable='diffuse_hemisp_narrowband_filter4', # Compute the FFT for each point +- window samples task = [] + sun_up = is_sun_visible(latitude=obj['lat'].values, longitude=obj['lon'].values, date_time=time) for t in range(len(time)): sind = t - fft_window eind = t + fft_window @@ -112,14 +108,13 @@ def fft_shading_test(obj, variable='diffuse_hemisp_narrowband_filter4', d = d[index] # Add to task for dask processing - lat = [obj['lat'].values] if not isinstance(obj['lat'].values, list) else obj['lat'].values - lon = [obj['lon'].values] if not isinstance(obj['lon'].values, list) else obj['lon'].values - task.append(dask.delayed(fft_shading_test_process)(time[t], - lat[0], lon[0], d, - shad_freq_lower=shad_freq_lower, - shad_freq_upper=shad_freq_upper, - ratio_thresh=ratio_thresh, - time_interval=dt)) + task.append(dask.delayed(fft_shading_test_process)( + time[t], d, + shad_freq_lower=shad_freq_lower, + shad_freq_upper=shad_freq_upper, + ratio_thresh=ratio_thresh, + time_interval=dt, + is_sunny=sun_up[t])) # Process using dask result = dask.compute(*task) @@ -161,9 +156,9 @@ def fft_shading_test(obj, variable='diffuse_hemisp_narrowband_filter4', return obj -def fft_shading_test_process(time, lat, lon, data, shad_freq_lower=None, +def fft_shading_test_process(time, data, shad_freq_lower=None, shad_freq_upper=None, ratio_thresh=None, - time_interval=None): + time_interval=None, is_sunny=None): """ Processing function to do the FFT calculations/thresholding @@ -171,10 +166,6 @@ def fft_shading_test_process(time, lat, lon, data, shad_freq_lower=None, ---------- time : datetime Center time of calculation used for calculating sunrise/sunset - lat : float - Latitude used for calculating sunrise/sunset - lon : float - Longitude used for calculating sunrise/sunset data : list Data for run through fft processing shad_freq_lower : list @@ -193,40 +184,7 @@ def fft_shading_test_process(time, lat, lon, data, shad_freq_lower=None, """ - # Get sunrise/sunset that are on same days - # This is used to help in the processing and easily exclude - # nighttime data from processing - if ASTRAL: - astral.solar_depression = 0 - obs = astral.Observer(latitude=float(lat), longitude=float(lon)) - s = sun(obs, pd.Timestamp(time)) - else: - a = astral.Astral() - a.solar_depression = 0 - s = a.sun_utc(pd.Timestamp(time), float(lat), float(lon)) - - sr = s['sunrise'].replace(tzinfo=None) - ss = s['sunset'].replace(tzinfo=None) - delta = ss.date() - sr.date() - if delta > datetime.timedelta(days=0): - if ASTRAL: - s = sun(obs, pd.Timestamp(time) - datetime.timedelta(days=1)) - else: - s = a.sun_utc(pd.Timestamp(time), float(lat), float(lon)) - ss = s['sunset'].replace(tzinfo=None) - - # Set if night or not - shading = 0 - night = False - if sr < ss: - if (pd.Timestamp(time) < sr) or (pd.Timestamp(time) > ss): - night = True - if sr > ss: - if (pd.Timestamp(time) < sr) and (pd.Timestamp(time) > ss): - night = True - - # Return shading of 0 if no valid data or it's night - if len(data) == 0 or night is True: + if not is_sunny: return {'shading': 0, 'fft': [np.nan] * len(data), 'freq': [np.nan] * len(data)} # FFT Algorithm @@ -244,6 +202,7 @@ def fft_shading_test_process(time, lat, lon, data, shad_freq_lower=None, # Return if FFT is empty if len(fftv) == 0: return {'shading': 0, 'fft': [np.nan] * len(data), 'freq': [np.nan] * len(data)} + # Commented out as it seems to work better without smoothing # fftv=pd.DataFrame(data=fftv).rolling(min_periods=3,window=3,center=True).mean().values.flatten() diff --git a/act/retrievals/radiation.py b/act/retrievals/radiation.py index 7f70403667..79e34c0ac6 100644 --- a/act/retrievals/radiation.py +++ b/act/retrievals/radiation.py @@ -10,12 +10,7 @@ import xarray as xr from scipy.constants import Stefan_Boltzmann from act.utils.datetime_utils import datetime64_to_datetime -import astral -try: - from astral import Observer - ASTRAL = True -except ImportError: - ASTRAL = False +from act.utils.geo_utils import get_solar_azimuth_elevation def calculate_dsh_from_dsdh_sdn(obj, dsdh='down_short_diffuse_hemisp', @@ -29,7 +24,7 @@ def calculate_dsh_from_dsdh_sdn(obj, dsdh='down_short_diffuse_hemisp', Parameters ---------- - obj : ACT object + obj : Xarray dataset Object where variables for these calculations are stored dsdh : str Name of the downwelling shortwave diffuse hemispheric irradiance field to use. @@ -38,33 +33,22 @@ def calculate_dsh_from_dsdh_sdn(obj, dsdh='down_short_diffuse_hemisp', Name of shortwave direct normal irradiance field to use. Defaults to shortwave_direct_normal_irradiance. lat : str - Name of SIRS lat field to use. Defaults to 'lat'. + Name of latitude field in dataset to use. Defaults to 'lat'. lon : str - Name of SIRS lon field to use. Defaults to 'lon'. + Name of longitued field in dataset to use. Defaults to 'lon'. Returns ------- - obj: ACT Object - Object with calculations included as new variables. + obj: Xarray dataset + ACT Xarray dataset oject with calculations included as new variables. """ # Calculating Derived Down Short Hemisp - - if ASTRAL: - obs = Observer(latitude=obj[lat], longitude=obj[lon]) - else: - a = astral.Astral() - tt = datetime64_to_datetime(obj['time'].values) - solar_zenith = np.full(len(tt), np.nan) - for ii, tm in enumerate(tt): - if ASTRAL: - solar_zenith[ii] = np.cos(np.radians(astral.sun.zenith(obs, tt[ii]))) - else: - solar_zenith[ii] = np.cos(np.radians(a.solar_zenith(tt[ii], obj[lat], obj[lon]))) - + elevation, _, _ = get_solar_azimuth_elevation(obj[lat].values, obj[lon].values, tt) + solar_zenith = np.cos(np.radians(90. - elevation)) dsh = (obj[dsdh].values + (solar_zenith * obj[sdn].values)) # Add data back to object diff --git a/act/tests/__init__.py b/act/tests/__init__.py index 1720597ab9..ba510fae3f 100644 --- a/act/tests/__init__.py +++ b/act/tests/__init__.py @@ -16,9 +16,9 @@ EXAMPLE_SONDE_WILDCARD """ -from .sample_files import (EXAMPLE_SONDE1, EXAMPLE_LCL1, +from .sample_files import (EXAMPLE_SONDE1, EXAMPLE_LCL1, EXAMPLE_MET_CSV, EXAMPLE_SONDE_WILDCARD, EXAMPLE_MET1, - EXAMPLE_METE40, + EXAMPLE_METE40, EXAMPLE_TWP_SONDE_20060121, EXAMPLE_MET_WILDCARD, EXAMPLE_CEIL1, EXAMPLE_CEIL_WILDCARD, EXAMPLE_ANL_CSV, EXAMPLE_MPL_1SAMPLE, EXAMPLE_IRT25m20s, @@ -27,4 +27,6 @@ EXAMPLE_EBBR2, EXAMPLE_BRS, EXAMPLE_AERI, EXAMPLE_MFRSR, EXAMPLE_SURFSPECALB1MLAWER, EXAMPLE_SIGMA_MPLV5, EXAMPLE_RL1, - EXAMPLE_CO2FLX4M, EXAMPLE_SIRS, EXAMPLE_IRTSST) + EXAMPLE_CO2FLX4M, EXAMPLE_SIRS, EXAMPLE_IRTSST, + EXAMPLE_MET_TEST1, EXAMPLE_MET_TEST2, + EXAMPLE_STAMP_WILDCARD) diff --git a/act/tests/baseline/test_2D_timeseries_plot.png b/act/tests/baseline/test_2D_timeseries_plot.png new file mode 100644 index 0000000000..9cd2e5c8dd Binary files /dev/null and b/act/tests/baseline/test_2D_timeseries_plot.png differ diff --git a/act/tests/baseline/test_2d_as_1d.png b/act/tests/baseline/test_2d_as_1d.png index 22fbe8871d..974430a033 100644 Binary files a/act/tests/baseline/test_2d_as_1d.png and b/act/tests/baseline/test_2d_as_1d.png differ diff --git a/act/tests/baseline/test_assessment_overplot.png b/act/tests/baseline/test_assessment_overplot.png index e436dd9861..fb67f04ae6 100644 Binary files a/act/tests/baseline/test_assessment_overplot.png and b/act/tests/baseline/test_assessment_overplot.png differ diff --git a/act/tests/baseline/test_assessment_overplot_multi.png b/act/tests/baseline/test_assessment_overplot_multi.png index 29c3f9bd47..dfcb33c658 100644 Binary files a/act/tests/baseline/test_assessment_overplot_multi.png and b/act/tests/baseline/test_assessment_overplot_multi.png differ diff --git a/act/tests/baseline/test_barb_sounding_plot.png b/act/tests/baseline/test_barb_sounding_plot.png index 9650493077..c79288b69b 100644 Binary files a/act/tests/baseline/test_barb_sounding_plot.png and b/act/tests/baseline/test_barb_sounding_plot.png differ diff --git a/act/tests/baseline/test_contour.png b/act/tests/baseline/test_contour.png index 63bf87cfd2..c2651dee65 100644 Binary files a/act/tests/baseline/test_contour.png and b/act/tests/baseline/test_contour.png differ diff --git a/act/tests/baseline/test_contour2.png b/act/tests/baseline/test_contour2.png new file mode 100644 index 0000000000..c6a9567e89 Binary files /dev/null and b/act/tests/baseline/test_contour2.png differ diff --git a/act/tests/baseline/test_contour_stamp.png b/act/tests/baseline/test_contour_stamp.png new file mode 100644 index 0000000000..25a17ba8f2 Binary files /dev/null and b/act/tests/baseline/test_contour_stamp.png differ diff --git a/act/tests/baseline/test_contourf.png b/act/tests/baseline/test_contourf.png new file mode 100644 index 0000000000..726e7a824b Binary files /dev/null and b/act/tests/baseline/test_contourf.png differ diff --git a/act/tests/baseline/test_contourf2.png b/act/tests/baseline/test_contourf2.png new file mode 100644 index 0000000000..3421f15ad7 Binary files /dev/null and b/act/tests/baseline/test_contourf2.png differ diff --git a/act/tests/baseline/test_fill_between.png b/act/tests/baseline/test_fill_between.png index 38609d88e2..5708721531 100644 Binary files a/act/tests/baseline/test_fill_between.png and b/act/tests/baseline/test_fill_between.png differ diff --git a/act/tests/baseline/test_geoplot.png b/act/tests/baseline/test_geoplot.png index 1b799c36ac..6d42bc13e0 100644 Binary files a/act/tests/baseline/test_geoplot.png and b/act/tests/baseline/test_geoplot.png differ diff --git a/act/tests/baseline/test_heatmap.png b/act/tests/baseline/test_heatmap.png index cede0ad207..f7a7a251aa 100644 Binary files a/act/tests/baseline/test_heatmap.png and b/act/tests/baseline/test_heatmap.png differ diff --git a/act/tests/baseline/test_multi_skewt_plot.png b/act/tests/baseline/test_multi_skewt_plot.png new file mode 100644 index 0000000000..4f222a56b8 Binary files /dev/null and b/act/tests/baseline/test_multi_skewt_plot.png differ diff --git a/act/tests/baseline/test_multidataset_plot_dict.png b/act/tests/baseline/test_multidataset_plot_dict.png index 5c50b0ac6d..f7856a5fa1 100644 Binary files a/act/tests/baseline/test_multidataset_plot_dict.png and b/act/tests/baseline/test_multidataset_plot_dict.png differ diff --git a/act/tests/baseline/test_multidataset_plot_tuple.png b/act/tests/baseline/test_multidataset_plot_tuple.png index 38d61b5d80..0be93f1a15 100644 Binary files a/act/tests/baseline/test_multidataset_plot_tuple.png and b/act/tests/baseline/test_multidataset_plot_tuple.png differ diff --git a/act/tests/baseline/test_plot.png b/act/tests/baseline/test_plot.png index ce1bfc21f6..4ab502ff5a 100644 Binary files a/act/tests/baseline/test_plot.png and b/act/tests/baseline/test_plot.png differ diff --git a/act/tests/baseline/test_plot_barbs_from_u_v.png b/act/tests/baseline/test_plot_barbs_from_u_v.png new file mode 100644 index 0000000000..0ac3c275a4 Binary files /dev/null and b/act/tests/baseline/test_plot_barbs_from_u_v.png differ diff --git a/act/tests/baseline/test_plot_barbs_from_u_v2.png b/act/tests/baseline/test_plot_barbs_from_u_v2.png new file mode 100644 index 0000000000..35feaa0f4a Binary files /dev/null and b/act/tests/baseline/test_plot_barbs_from_u_v2.png differ diff --git a/act/tests/baseline/test_qc_bar_plot.png b/act/tests/baseline/test_qc_bar_plot.png index 717f9a139e..9aa36b8d0a 100644 Binary files a/act/tests/baseline/test_qc_bar_plot.png and b/act/tests/baseline/test_qc_bar_plot.png differ diff --git a/act/tests/baseline/test_qc_flag_block_plot.png b/act/tests/baseline/test_qc_flag_block_plot.png index a01f34b23f..fa3e328303 100644 Binary files a/act/tests/baseline/test_qc_flag_block_plot.png and b/act/tests/baseline/test_qc_flag_block_plot.png differ diff --git a/act/tests/baseline/test_size_distribution.png b/act/tests/baseline/test_size_distribution.png index ecc4dfc61c..2fa19f74d9 100644 Binary files a/act/tests/baseline/test_size_distribution.png and b/act/tests/baseline/test_size_distribution.png differ diff --git a/act/tests/baseline/test_skewt_plot.png b/act/tests/baseline/test_skewt_plot.png index 96e70a7080..3fc8d67b0f 100644 Binary files a/act/tests/baseline/test_skewt_plot.png and b/act/tests/baseline/test_skewt_plot.png differ diff --git a/act/tests/baseline/test_skewt_plot_spd_dir.png b/act/tests/baseline/test_skewt_plot_spd_dir.png index 96e70a7080..3fc8d67b0f 100644 Binary files a/act/tests/baseline/test_skewt_plot_spd_dir.png and b/act/tests/baseline/test_skewt_plot_spd_dir.png differ diff --git a/act/tests/baseline/test_stacked_bar_graph.png b/act/tests/baseline/test_stacked_bar_graph.png index 78614ea9ac..bd52279209 100644 Binary files a/act/tests/baseline/test_stacked_bar_graph.png and b/act/tests/baseline/test_stacked_bar_graph.png differ diff --git a/act/tests/baseline/test_stacked_bar_graph2.png b/act/tests/baseline/test_stacked_bar_graph2.png new file mode 100644 index 0000000000..13e1fc78ef Binary files /dev/null and b/act/tests/baseline/test_stacked_bar_graph2.png differ diff --git a/act/tests/baseline/test_stacked_bar_graph_sorted.png b/act/tests/baseline/test_stacked_bar_graph_sorted.png index 0c031704f3..4607ad3d27 100644 Binary files a/act/tests/baseline/test_stacked_bar_graph_sorted.png and b/act/tests/baseline/test_stacked_bar_graph_sorted.png differ diff --git a/act/tests/baseline/test_stair_graph.png b/act/tests/baseline/test_stair_graph.png index 9ceaaa08bd..8c0e0fb4c6 100644 Binary files a/act/tests/baseline/test_stair_graph.png and b/act/tests/baseline/test_stair_graph.png differ diff --git a/act/tests/baseline/test_stair_graph_sorted.png b/act/tests/baseline/test_stair_graph_sorted.png index 1d671d85e2..5d2f888b0c 100644 Binary files a/act/tests/baseline/test_stair_graph_sorted.png and b/act/tests/baseline/test_stair_graph_sorted.png differ diff --git a/act/tests/baseline/test_wind_rose.png b/act/tests/baseline/test_wind_rose.png index fba043e82b..95e6b4c1ee 100644 Binary files a/act/tests/baseline/test_wind_rose.png and b/act/tests/baseline/test_wind_rose.png differ diff --git a/act/tests/baseline/test_xsection_plot.png b/act/tests/baseline/test_xsection_plot.png index 13edf39bb0..167de71687 100644 Binary files a/act/tests/baseline/test_xsection_plot.png and b/act/tests/baseline/test_xsection_plot.png differ diff --git a/act/tests/baseline/test_xsection_plot_map.png b/act/tests/baseline/test_xsection_plot_map.png index b0104f8af1..76413c84e7 100644 Binary files a/act/tests/baseline/test_xsection_plot_map.png and b/act/tests/baseline/test_xsection_plot_map.png differ diff --git a/act/tests/data/201509021500.bin b/act/tests/data/201509021500.bi similarity index 100% rename from act/tests/data/201509021500.bin rename to act/tests/data/201509021500.bi diff --git a/act/tests/data/brw21001.dat b/act/tests/data/brw21001.dat new file mode 100644 index 0000000000..4ac35acc84 --- /dev/null +++ b/act/tests/data/brw21001.dat @@ -0,0 +1,20 @@ + Barrow + 71.316 -156.600 11 m + 2021 1 1 1 0 0 0.000 95.60 -1.5 0 -0.7 0 0.4 0 0.0 0 139.8 0 246.15 0 245.91 0 209.5 0 246.15 0 246.03 0 -9999.9 2 -9999.9 1 0.0 0 -69.6 0 -69.6 0 -28.2 0 78.6 0 6.8 0 78.1 0 1020.1 0 + 2021 1 1 1 0 1 0.017 95.62 -1.5 0 -0.7 0 0.4 0 0.0 0 139.8 0 246.15 0 245.91 0 209.5 0 246.15 0 246.03 0 -9999.9 2 -9999.9 1 0.0 0 -69.7 0 -69.7 0 -28.2 0 78.6 0 6.5 0 77.8 0 1020.1 0 + 2021 1 1 1 0 2 0.033 95.65 -1.5 0 -0.8 0 0.3 0 0.0 0 139.4 0 246.15 0 245.91 0 209.8 0 246.18 0 246.03 0 -9999.9 2 -9999.9 1 0.0 0 -70.5 0 -70.4 0 -28.3 0 78.8 0 6.8 0 77.2 0 1020.1 0 + 2021 1 1 1 0 3 0.050 95.68 -1.5 0 -0.8 0 0.3 0 0.0 0 139.6 0 246.15 0 245.88 0 210.0 0 246.19 0 246.03 0 -9999.9 2 -9999.9 1 0.0 0 -70.4 0 -70.4 0 -28.3 0 79.2 0 6.2 0 76.1 0 1020.1 0 + 2021 1 1 1 0 4 0.067 95.71 -1.5 0 -0.8 0 0.4 0 0.0 0 139.5 0 246.15 0 245.84 0 210.0 0 246.19 0 246.03 0 -9999.9 2 -9999.9 1 0.0 0 -70.5 0 -70.4 0 -28.4 0 79.2 0 6.1 0 75.3 0 1020.0 0 + 2021 1 1 1 0 5 0.083 95.74 -1.5 0 -0.9 0 0.3 0 0.0 0 139.4 0 246.14 0 245.81 0 210.0 0 246.19 0 246.03 0 -9999.9 2 -9999.9 1 0.0 0 -70.6 0 -70.6 0 -28.4 0 79.1 0 6.3 0 74.7 0 1019.9 0 + 2021 1 1 1 0 6 0.100 95.77 -1.5 0 -0.9 0 0.1 0 -0.1 0 139.3 0 246.11 0 245.77 0 210.3 0 246.21 0 246.03 0 -9999.9 2 -9999.9 1 0.0 0 -71.0 0 -71.0 0 -28.5 0 79.0 0 6.9 0 74.7 0 1019.9 0 + 2021 1 1 1 0 7 0.117 95.80 -1.5 0 -0.9 0 0.2 0 -0.1 0 139.6 0 246.10 0 245.75 0 210.3 0 246.22 0 246.03 0 -9999.9 2 -9999.9 1 0.0 0 -70.7 0 -70.7 0 -28.5 0 79.3 0 7.1 0 74.5 0 1019.8 0 + 2021 1 1 1 0 8 0.133 95.83 -1.5 0 -0.9 0 0.4 0 -0.1 0 142.4 0 246.07 0 245.75 0 210.4 0 246.22 0 246.03 0 -9999.9 2 -9999.9 1 0.0 0 -68.0 0 -68.0 0 -28.4 0 78.4 0 6.6 0 74.3 0 1019.8 0 + 2021 1 1 1 0 9 0.150 95.86 -1.5 0 -0.9 0 0.4 0 -0.1 0 146.6 0 246.07 0 245.75 0 210.5 0 246.22 0 246.03 0 -9999.9 2 -9999.9 1 0.0 0 -63.9 0 -63.9 0 -28.4 0 77.7 0 6.8 0 74.2 0 1019.8 0 + 2021 1 1 1 0 10 0.167 95.89 -1.5 0 -0.9 0 0.4 0 -0.1 0 146.6 0 246.04 0 245.75 0 210.2 0 246.22 0 246.05 0 -9999.9 2 -9999.9 1 0.0 0 -63.6 0 -63.6 0 -28.3 0 77.9 0 6.7 0 74.2 0 1019.8 0 + 2021 1 1 1 0 11 0.183 95.92 -1.5 0 -0.9 0 0.4 0 -0.1 0 146.7 0 246.03 0 245.75 0 210.1 0 246.23 0 246.07 0 -9999.9 2 -9999.9 1 0.0 0 -63.4 0 -63.4 0 -28.3 0 77.7 0 6.8 0 74.1 0 1019.8 0 + 2021 1 1 1 0 12 0.200 95.95 -1.5 0 -0.9 0 0.4 0 0.0 0 148.8 0 246.03 0 245.78 0 210.1 0 246.23 0 246.07 0 -9999.9 2 -9999.9 1 0.0 0 -61.3 0 -61.3 0 -28.2 0 77.8 0 7.2 0 74.7 0 1019.8 0 + 2021 1 1 1 0 13 0.217 95.99 -1.5 0 -0.8 0 0.4 0 0.0 0 149.1 0 246.03 0 245.82 0 209.9 0 246.23 0 246.09 0 -9999.9 2 -9999.9 1 0.0 0 -60.8 0 -60.8 0 -28.1 0 77.2 0 6.5 0 75.0 0 1019.8 0 + 2021 1 1 1 0 14 0.233 96.02 -1.5 0 -0.8 0 1.0 0 0.0 0 151.8 0 246.03 0 245.87 0 209.6 0 246.23 0 246.13 0 -9999.9 2 -9999.9 1 0.0 0 -57.7 0 -57.7 0 -28.0 0 76.7 0 7.5 0 75.4 0 1019.8 0 + 2021 1 1 1 0 15 0.250 96.05 -1.5 0 -0.8 0 1.1 0 0.1 0 155.4 0 246.03 0 245.94 0 209.5 0 246.25 0 246.15 0 -9999.9 2 -9999.9 1 0.0 0 -54.1 0 -54.1 0 -27.9 0 76.6 0 6.7 0 76.3 0 1019.8 0 + 2021 1 1 1 0 16 0.267 96.08 -1.5 0 -0.7 0 1.1 0 0.2 0 157.6 0 246.07 0 246.00 0 209.4 0 246.27 0 246.19 0 -9999.9 2 -9999.9 1 0.0 0 -51.8 0 -51.8 0 -27.8 0 76.6 0 7.2 0 76.6 0 1019.8 0 + 2021 1 1 1 0 17 0.283 96.12 -1.2 0 -0.7 0 1.1 0 0.2 0 161.5 0 246.10 0 246.05 0 209.5 0 246.30 0 246.22 0 -9999.9 2 -9999.9 1 0.0 0 -48.0 0 -48.0 0 -27.7 0 76.7 0 6.6 0 76.9 0 1019.8 0 diff --git a/act/tests/data/brw_12_2020_hour.dat b/act/tests/data/brw_12_2020_hour.dat new file mode 100644 index 0000000000..a56ced2c56 --- /dev/null +++ b/act/tests/data/brw_12_2020_hour.dat @@ -0,0 +1,34 @@ +These data are made available to the scientific community and public with the understanding that authors of publications +and presentations will contact the PI to obtain up-to-date information on the latest version of the data during the early stages of their analysis. +Manuscripts that employ these data should be reviewed by the PI before submission to ensure the quality and limitations of the data are represented. + +Data Disclaimer: Data has undergone extensive quality checks and is available from NOAA-GMD and WDCGG. +However, there still exists the potential for these data to be modified at the discretion of NOAA/ESRL Global Monitoring Division. + +If you have questions regarding the data in this file, please contact: +Irina Petropavlovskikh: Irina.Petro@noaa.gov (303-497-6279) +Audra McClure-Begley: Audra.mcclure@noaa.gov + +Please use the following citation for data use: +McClure-Begley, A., Petropavlovskikh, I., Oltmans, S., (2014) NOAA Global Monitoring Surface Ozone Network. + Station name, Time start-Time end. National Oceanic and Atmospheric Administration, Earth Systems Research Laboratory Global Monitoring Division. + Boulder, CO. DATE ACCESSED http://dx.doi.org/10.7289/V57P8WBF + +STN YEAR MON DAY HR O3(PPB) +BRW 2020 12 01 00 32.15 +BRW 2020 12 01 01 33.07 +BRW 2020 12 01 02 34.06 +BRW 2020 12 01 03 35.18 +BRW 2020 12 01 04 35.41 +BRW 2020 12 01 05 35.40 +BRW 2020 12 01 06 35.23 +BRW 2020 12 01 07 35.38 +BRW 2020 12 01 08 35.57 +BRW 2020 12 01 09 36.04 +BRW 2020 12 01 10 36.13 +BRW 2020 12 01 11 35.37 +BRW 2020 12 01 12 33.13 +BRW 2020 12 01 13 32.62 +BRW 2020 12 01 14 32.53 +BRW 2020 12 01 15 32.51 +BRW 2020 12 01 16 32.98 diff --git a/act/tests/data/brw_CCl4_Day.dat b/act/tests/data/brw_CCl4_Day.dat new file mode 100644 index 0000000000..4fd2044e44 --- /dev/null +++ b/act/tests/data/brw_CCl4_Day.dat @@ -0,0 +1,68 @@ +# file: brw_CCl4_Day.dat +# date: Mon, Mar 8, 2021 2:35:30 PM +# +# Carbon tetrachloride (CCl4) data from hourly in situ samples analyzed on a gas chromatograph located +# at Pt. Barrow (BRW), Alaska (71.3 N, 156.6 W, elevation: 8 m). +# +# Atmospheric Measurements from the NOAA/ESRL Chromatograph for Atmospheric Trace Species +# (CATS) Program. This GC replaces the RITS GCs operated during the period of 1986 through 2000. +# This work was funded in part by the Atmospheric Chemistry Project of NOAA's Climate and +# Global Change Program. +# +# (PIs) Geoffrey S. Dutton, Dr. James W. Elkins, Dr. Bradley D. Hall +# +# National Oceanic and Atmospheric Administration (NOAA) +# Earth System Research Laboratory (ESRL) +# 325 Broadway, R/GML +# Boulder, CO 80305-3328 +# +# Email: Geoff.Dutton@noaa.gov; James.W.Elkins@noaa.gov; Bradley.Hall@noaa.gov +# Phone: (303) 497-6086 (303) 497-6224 (303) 497-7011 +# Fax: (303) 497-6290 +# +# See (http://www.esrl.noaa.gov/gmd/obop/) for station statistics and personnel. +# +# Data use policy: +# If you use this data in a publication or report, we would appreciate that you contact the principal +# investigators (PIs) first and discuss your interests. We are trying to encourage scientific discussion. +# The PIs can discuss the quality of the data and any potential problems with its analysis. Recent data (less +# than 1.5 years) are considered preliminary. Please contact the PIs for more information. +# +# The data should be cited in presentations or publications as the flowing: Carbon tetrachloride data from +# the NOAA/ESRL halocarbons in situ program. +# +# DOI: http://doi.org/10.7289/V5X0659V +# +# Calibration scale used: NOAA 2008 +# +# More information about the calibration scale can be found at: +# http://www.esrl.noaa.gov/gmd/ccl/ +# +# The CATS GCs sample air once an hour. Daily median data are provided in parts-per-trillion, ppt. +# Abbreviations used: +# yr=year, mon=month, m=monthly data; gc=type of sample--in situ GC data, sd=standard deviation of samples, +# n=number of samples, brw=Pt. Barrow, Alaska, sum=Summit, Greenland, nwr= Niwot Ridge, Colorado +# mlo=Mauna Loa, Hawaii, smo = American Samoa, spo=South Pole, Antarctica +# NAN=not a number or no data, and space(s) is the delimiter. +# +# The standard deviation (unc) reported is the estimate of insturmental precision. +# The standard deviation (sd) reported is the estimate of atmospheric variability. +# +# year, month, day, daily median CCl4 in ppt, std. dev. +# +CCl4catsBRWyr CCl4catsBRWmon CCl4catsBRWday CCl4catsBRWm CCl4catsBRWmsd CCl4catsBRWn +1998 6 16 nan nan 0 +1998 6 17 103.76 0.66 21 +1998 6 18 103.44 0.59 24 +1998 6 19 103.27 0.45 17 +1998 6 20 103.37 0.41 20 +1998 6 21 103.20 0.43 19 +1998 6 22 103.11 0.39 18 +1998 6 23 103.20 0.40 24 +1998 6 24 103.11 0.52 22 +1998 6 25 103.07 0.50 17 +1998 6 26 103.35 0.51 15 +1998 6 27 103.09 0.60 24 +1998 6 28 103.09 0.66 14 +1998 6 29 nan nan 0 +1998 6 30 103.59 0.56 12 diff --git a/act/tests/data/co2_brw_surface-insitu_1_ccgg_MonthlyData.txt b/act/tests/data/co2_brw_surface-insitu_1_ccgg_MonthlyData.txt new file mode 100644 index 0000000000..09112800a2 --- /dev/null +++ b/act/tests/data/co2_brw_surface-insitu_1_ccgg_MonthlyData.txt @@ -0,0 +1,164 @@ +# header_lines : 151 +# +# ------------------------------------------------------------->>>> +# DATA SET NAME +# +# dataset_name: co2_brw_surface-insitu_1_ccgg_MonthlyData +# +# ------------------------------------------------------------->>>> +# DESCRIPTION +# +# dataset_description: Atmospheric Carbon Dioxide Dry Air Mole Fractions from quasi-continuous measurements at Barrow, Alaska. +# +# ------------------------------------------------------------->>>> +# CITATION +# +# dataset_citation: K.W. Thoning, A.M. Crotwell, and J.W. Mund (2021), Atmospheric Carbon Dioxide Dry Air Mole Fractions from continuous measurements at Mauna Loa, Hawaii, Barrow, Alaska, American Samoa and South Pole. 1973-2019, Version 2021-02 National Oceanic and Atmospheric Administration (NOAA), Global Monitoring Laboratory (GML), Boulder, Colorado, USA https://doi.org/10.15138/yaf1-bk21 FTP path: ftp://aftp.cmdl.noaa.gov/data/greenhouse_gases/co2/in-situ/surface/ +# +# ------------------------------------------------------------->>>> +# FAIR USE POLICY +# +# dataset_fair_use: These data are made freely available to the public and the scientific community in the belief that their wide dissemination will lead to greater understanding and new scientific insights. The availability of these data does not constitute publication of the data. NOAA relies on the ethics and integrity of the user to ensure that ESRL receives fair credit for their work. If the data are obtained for potential use in a publication or presentation, ESRL should be informed at the outset of the nature of this work. If the ESRL data are essential to the work, or if an important result or conclusion depends on the ESRL data, co-authorship may be appropriate. This should be discussed at an early stage in the work. Manuscripts using the ESRL data should be sent to ESRL for review before they are submitted for publication so we can ensure that the quality and limitations of the data are accurately represented. +# +# ------------------------------------------------------------->>>> +# WARNING +# +# dataset_warning: Every effort is made to produce the most accurate and precise measurements possible. However, we reserve the right to make corrections to the data based on recalibration of standard gases or for other reasons deemed scientifically justified. We are not responsible for results and conclusions based on use of these data without regard to this warning. +# +# ------------------------------------------------------------->>>> +# GLOBAL ATTRIBUTES +# +# site_code : BRW +# site_name : Barrow Atmospheric Baseline Observatory +# site_country : United States +# site_country_flag : http://www.esrl.noaa.gov/gmd/webdata/ccgg/ObsPack/images/flags/UNST0001.GIF +# site_latitude : 71.323 +# site_longitude : -156.6114 +# site_elevation : 11.0 +# site_elevation_unit : masl +# site_position_comment : This is the nominal location of the site. The sampling location at many sites has changed over time, and we report here the most recent nominal location. The actual sampling location for each observation is not necessarily the site location. The sampling locations for each observation are reported in the latitude, longitude, and altitude variables. +# site_utc2lst : -9.0 +# site_utc2lst_comment : Add 'site_utc2lst' hours to convert a time stamp in UTC (Coordinated Universal Time) to LST (Local Standard Time). +# site_url : http://www.esrl.noaa.gov/gmd/obop/brw/index.html +# dataset_creation_date : 2021-03-05T12:36:08.450229 +# dataset_num : 3 +# dataset_name : co2_brw_surface-insitu_1_ccgg_MonthlyData +# dataset_parameter : co2 +# dataset_project : surface-insitu +# dataset_platform : fixed +# dataset_selection : Monthly values derived from hourly data representative of baseline conditions +# dataset_selection_tag : MonthlyData +# dataset_calibration_scale : WMO CO2 X2019 +# dataset_start_date : 1973-01-01T00:00:00Z +# dataset_stop_date : 2019-12-01T00:00:00Z +# dataset_data_frequency : 1 +# dataset_data_frequency_unit : month +# dataset_reciprocity : Use of these data implies an agreement to reciprocate. Laboratories making similar measurements agree to make their own data available to the general public and to the scientific community in an equally complete and easily accessible form. Modelers are encouraged to make available to the community, upon request, their own tools used in the interpretation of the ESRL data, namely well documented model code, transport fields, and additional information necessary for other scientists to repeat the work and to run modified versions. Model availability includes collaborative support for new users of the models. +# dataset_usage_url : https://www.esrl.noaa.gov/gmd/ccgg/obspack/citation.php?product=obspack_co2_1_CCGGObsInsitu_v1.0_2021-03-05 +# dataset_usage_description : Please cite the product's citation when using data from this dataset. Relevant literature references for this dataset are listed below for convenience. +# dataset_reference_total_listed : 2 +# dataset_reference_1_name : Peterson, J.T., W.D. Komhyr, L.S. Waterman, R.H. Gammon, K.W. Thoning, and T.J. Conway, Atmospheric CO2 variations at Barrow, Alaska, 1973 1982, J. Atmos. Chem., 4, 491 510, 1986. +# dataset_reference_2_name : Halter, B. and J. M. Harris (1983), On the variability of atmospheric carbon dioxide concentration at Barrow, Alaska, during winter, Journal of Geophysical Research, 88(C11), 6858-6864. +# dataset_contribution : These data are provided by NOAA. Principal investigators include Kirk Thoning (NOAA) AND Pieter Tans (NOAA). +# lab_total_listed : 1 +# lab_1_number : 1 +# lab_1_abbr : NOAA +# lab_1_name : NOAA Global Monitoring Laboratory +# lab_1_address1 : 325 Broadway +# lab_1_address2 : NOAA R/GML-1 +# lab_1_address3 : Boulder, CO 80305-3328 +# lab_1_country : United States +# lab_1_url : http://www.esrl.noaa.gov/gmd/ccgg/ +# lab_1_parameter : Lab has contributed measurements for: co2 +# lab_1_country_flag : http://www.esrl.noaa.gov/gmd/webdata/ccgg/ObsPack/images/flags/UNST0001.GIF +# lab_1_logo : http://www.esrl.noaa.gov/gmd/webdata/ccgg/ObsPack/images/logos/noaa_medium.png +# lab_1_ongoing_atmospheric_air_comparison : T +# lab_1_comparison_activity : Ongoing comparison with co-located measurements including NOAA surface flask data and independent measurement laboratories. +# provider_total_listed : 2 +# provider_1_name : Kirk Thoning +# provider_1_address1 : NOAA ESRL GML +# provider_1_address2 : 325 Broadway R/GML-1 +# provider_1_address3 : Boulder, CO 80305-3328 +# provider_1_country : United States +# provider_1_affiliation : National Oceanic and Atmospheric Administration +# provider_1_email : kirk.w.thoning@noaa.gov +# provider_1_tel : 303-497-6078 +# provider_1_parameter : Provider has contributed measurements for: co2 +# provider_2_name : Pieter Tans +# provider_2_address1 : NOAA ESRL GML +# provider_2_address2 : 325 Broadway R/GML-1 +# provider_2_address3 : Boulder, CO 80305-3328 +# provider_2_country : United States +# provider_2_affiliation : National Oceanic and Atmospheric Administration +# provider_2_email : pieter.tans@noaa.gov +# provider_2_tel : 303-497-6678 +# provider_2_parameter : Provider has contributed measurements for: co2 +# ------------------------------------------------------------->>>> +# VARIABLE ATTRIBUTES +# +# site_code:long_name : site_name_abbreviation. +# site_code:comment : Site code is an abbreviation for the sampling site name. +# time_components:_FillValue : -9 +# time_components:long_name : integer_components_of_UTC_date/time +# time_components:order : year, month, day, hour, minute, second +# time_components:comment : Calendar time components as integers. Times and dates are UTC. Time-averaged values are reported at the beginning of the averaging interval. +# time_decimal:_FillValue : -999.999 +# time_decimal:long_name : sample_decimal_year_in_UTC +# time_decimal:comment : decimal year in UTC. Time-averaged values are reported at the beginning of the averaging interval. +# value:_FillValue : -999.999 +# value:long_name : measured_mole_fraction_of_trace_gas_in_dry_air +# value:units : micromol mol-1 +# value:comment : Mole fraction reported in units of micromol mol-1 (10-6 mol per mol of dry air); abbreviated as ppm (parts per million). +# value_std_dev:_FillValue : -99.99 +# value_std_dev:long_name : standard_deviation_in_reported_value +# value_std_dev:units : micromol mol-1 +# value_std_dev:comment : This is the standard deviation of the reported mean value when nvalue is greater than 1. See provider_comment if available. +# nvalue:_FillValue : -9 +# nvalue:long_name : number_of_measurements_contributing_to_reported_value +# nvalue:comment : Number of individual measurements used to compute reported value. See provider_comment if available. +# latitude:_FillValue : -999.999 +# latitude:standard_name : latitude +# latitude:long_name : sample_latitude_in_decimal_degrees +# latitude:units : degrees_north +# latitude:comment : Latitude at which air sample was collected. +# longitude:_FillValue : -999.999 +# longitude:standard_name : longitude +# longitude:long_name : sample_longitude_in_decimal_degrees +# longitude:units : degrees_east +# longitude:comment : Longitude at which air sample was collected using a range of -180 degrees to +180 degrees. +# altitude:_FillValue : -999.999 +# altitude:standard_name : altitude +# altitude:long_name : sample_altitude_in_meters_above_sea_level +# altitude:units : m +# altitude:comment : Altitude (in meters above sea level). See provider_comment if available. +# altitude:provider_comment : Altitude for this dataset is the sum of surface elevation (masl) and sample intake height (magl). +# elevation:_FillValue : -999.999 +# elevation:standard_name : elevation +# elevation:long_name : surface_elevation_in_meters_above_sea_level +# elevation:units : m +# elevation:comment : Surface elevation in meters above sea level. See provider_comment if available. +# intake_height:_FillValue : -999.999 +# intake_height:long_name : sample_intake_height_in_meters_above_ground_level +# intake_height:units : m +# intake_height:comment : Sample intake height in meters above ground level (magl). See provider_comment if available. +# qcflag:missing_value : NA +# qcflag:long_name : quality_control_flag +# qcflag:comment : This quality control flag is provided by the contributing PIs. See provider_comment if available. +# qcflag:provider_comment : This is the NOAA 3-character quality control flag. Column 1 is the REJECTION flag. An alphanumeric other than a period (.) in the FIRST column indicates a sample with obvious problems during collection or analysis. This measurement should not be interpreted. Column 2 is the SELECTION flag. An alphanumeric other than a period (.) in the SECOND column indicates a sample that is likely valid but does not meet selection criteria determined by the goals of a particular investigation. Column 3 is the INFORMATION flag. An alphanumeric other than a period (.) in the THIRD column provides additional information about the collection or analysis of the sample. +# +# VARIABLE ORDER +# +site_code year month day hour minute second time_decimal value value_std_dev nvalue latitude longitude altitude elevation intake_height qcflag +BRW 1973 1 1 0 0 0 1973.0 -999.99 -99.99 0 71.323 -156.611 27.0 11.0 16.0 *.. +BRW 1973 2 1 0 0 0 1973.0849315068492 -999.99 -99.99 0 71.323 -156.611 27.0 11.0 16.0 *.. +BRW 1973 3 1 0 0 0 1973.1616438356164 -999.99 -99.99 0 71.323 -156.611 27.0 11.0 16.0 *.. +BRW 1973 4 1 0 0 0 1973.2465753424658 -999.99 -99.99 0 71.323 -156.611 27.0 11.0 16.0 *.. +BRW 1973 5 1 0 0 0 1973.3287671232877 -999.99 -99.99 0 71.323 -156.611 27.0 11.0 16.0 *.. +BRW 1973 6 1 0 0 0 1973.4136986301369 -999.99 -99.99 0 71.323 -156.611 27.0 11.0 16.0 *.. +BRW 1973 7 1 0 0 0 1973.495890410959 324.53 0.39 5 71.323 -156.611 27.0 11.0 16.0 ... +BRW 1973 8 1 0 0 0 1973.5808219178082 322.73 0.84 23 71.323 -156.611 27.0 11.0 16.0 ... +BRW 1973 9 1 0 0 0 1973.6657534246576 324.79 1.81 29 71.323 -156.611 27.0 11.0 16.0 ... +BRW 1973 10 1 0 0 0 1973.7479452054795 329.86 2.78 29 71.323 -156.611 27.0 11.0 16.0 ... +BRW 1973 11 1 0 0 0 1973.8328767123287 333.61 0.77 27 71.323 -156.611 27.0 11.0 16.0 ... +BRW 1973 12 1 0 0 0 1973.9150684931508 334.6 1.6 30 71.323 -156.611 27.0 11.0 16.0 ... +BRW 1974 1 1 0 0 0 1974.0 337.51 1.35 20 71.323 -156.611 27.0 11.0 16.0 ... diff --git a/act/tests/data/met_brw_insitu_1_obop_hour_2020.txt b/act/tests/data/met_brw_insitu_1_obop_hour_2020.txt new file mode 100644 index 0000000000..e6b17550f0 --- /dev/null +++ b/act/tests/data/met_brw_insitu_1_obop_hour_2020.txt @@ -0,0 +1,20 @@ +BRW 2020 01 01 00 5 7.1 100 1004.73 -25.5 -25.0 -999.9 77 -99 +BRW 2020 01 01 01 0 7.0 100 1005.16 -25.8 -25.2 -999.9 77 -99 +BRW 2020 01 01 02 3 7.6 99 1005.48 -25.3 -25.0 -999.9 77 -99 +BRW 2020 01 01 03 4 7.6 99 1006.04 -25.2 -24.7 -999.9 77 -99 +BRW 2020 01 01 04 15 7.5 100 1006.61 -25.1 -24.6 -999.9 77 -99 +BRW 2020 01 01 05 18 7.2 100 1007.19 -24.9 -24.5 -999.9 77 -99 +BRW 2020 01 01 06 17 7.7 100 1007.56 -24.8 -24.5 -999.9 77 -99 +BRW 2020 01 01 07 23 7.0 100 1008.24 -25.4 -25.0 -999.9 77 -99 +BRW 2020 01 01 08 21 7.4 100 1008.89 -25.7 -25.3 -999.9 77 -99 +BRW 2020 01 01 09 16 6.9 100 1009.78 -25.8 -25.2 -999.9 76 -99 +BRW 2020 01 01 10 17 7.0 100 1010.39 -25.3 -24.9 -999.9 77 -99 +BRW 2020 01 01 11 13 7.0 100 1010.57 -25.4 -25.0 -999.9 77 -99 +BRW 2020 01 01 12 28 7.5 100 1010.80 -24.9 -24.7 -999.9 77 -99 +BRW 2020 01 01 13 18 7.5 100 1011.54 -24.7 -24.5 -999.9 77 -99 +BRW 2020 01 01 14 17 7.6 100 1012.06 -24.3 -24.2 -999.9 77 -99 +BRW 2020 01 01 15 16 7.8 100 1012.52 -24.4 -24.1 -999.9 77 -99 +BRW 2020 01 01 16 44 6.9 100 1013.36 -24.7 -24.4 -999.9 77 -99 +BRW 2020 01 01 17 36 7.6 100 1014.02 -24.8 -24.5 -999.9 77 -99 +BRW 2020 01 01 18 28 7.8 100 1014.64 -24.3 -24.1 -999.9 78 -99 +BRW 2020 01 01 19 25 8.4 100 1015.17 -24.4 -24.2 -999.9 78 -99 diff --git a/act/tests/data/sgpmetE13.b1.20210401.000000.csv b/act/tests/data/sgpmetE13.b1.20210401.000000.csv new file mode 100644 index 0000000000..9fc8eff6ab --- /dev/null +++ b/act/tests/data/sgpmetE13.b1.20210401.000000.csv @@ -0,0 +1,4 @@ +time,date_time,time_offset,atmos_pressure,qc_atmos_pressure,temp_mean,qc_temp_mean,temp_std,rh_mean,qc_rh_mean,rh_std,vapor_pressure_mean,qc_vapor_pressure_mean,vapor_pressure_std,wspd_arith_mean,qc_wspd_arith_mean,wspd_vec_mean,qc_wspd_vec_mean,wdir_vec_mean,qc_wdir_vec_mean,wdir_vec_std,tbrg_precip_total,qc_tbrg_precip_total,tbrg_precip_total_corr,qc_tbrg_precip_total_corr,org_precip_rate_mean,qc_org_precip_rate_mean,pwd_err_code,pwd_mean_vis_1min,qc_pwd_mean_vis_1min,pwd_mean_vis_10min,qc_pwd_mean_vis_10min,pwd_pw_code_inst,qc_pwd_pw_code_inst,pwd_pw_code_15min,qc_pwd_pw_code_15min,pwd_pw_code_1hr,qc_pwd_pw_code_1hr,pwd_precip_rate_mean_1min,qc_pwd_precip_rate_mean_1min,pwd_cumul_rain,qc_pwd_cumul_rain,pwd_cumul_snow,qc_pwd_cumul_snow,logger_volt,qc_logger_volt,logger_temp,qc_logger_temp,lat,lon,alt +2021-04-01 00:00:00,2021-04-01 00:00:00,2021-04-01 00:00:00,99.14,0,13.9,0,0.019,36.49,0,0.231,0.579,0,0.003,2.744,0,2.742,0,113.1,0,2.301,0.0,0,0.0,0,0.0,0,0.0,20000.0,0,20000.0,0,0.0,0,0.0,0,0.0,0,0.0,0,26.98,0,457.0,0,12.83,0,20.4,0,36.605,-97.485,318.0 +2021-04-01 00:01:00,2021-04-01 00:00:00,2021-04-01 00:01:00,99.14,0,13.86,0,0.026,37.29,0,0.163,0.59,0,0.002,2.553,0,2.549,0,110.2,0,3.345,0.0,0,0.0,0,0.0,0,0.0,20000.0,0,20000.0,0,0.0,0,0.0,0,0.0,0,0.0,0,26.98,0,457.0,0,12.84,0,20.37,0,36.605,-97.485,318.0 +2021-04-01 00:02:00,2021-04-01 00:00:00,2021-04-01 00:02:00,99.14,0,13.76,0,0.034,38.36,0,0.524,0.603,0,0.007,2.661,0,2.66,0,103.6,0,1.518,0.0,0,0.0,0,0.0,0,0.0,20000.0,0,20000.0,0,0.0,0,0.0,0,0.0,0,0.0,0,26.98,0,457.0,0,12.83,0,20.31,0,36.605,-97.485,318.0 diff --git a/act/tests/data/sgpmet_no_time.nc b/act/tests/data/sgpmet_no_time.nc new file mode 100644 index 0000000000..cf32a36d83 Binary files /dev/null and b/act/tests/data/sgpmet_no_time.nc differ diff --git a/act/tests/data/sgpmet_test_time.nc b/act/tests/data/sgpmet_test_time.nc new file mode 100644 index 0000000000..3f067ea715 Binary files /dev/null and b/act/tests/data/sgpmet_test_time.nc differ diff --git a/act/tests/data/sgpstampE13.b1.20200101.000000.nc b/act/tests/data/sgpstampE13.b1.20200101.000000.nc new file mode 100644 index 0000000000..2bb4d93df3 Binary files /dev/null and b/act/tests/data/sgpstampE13.b1.20200101.000000.nc differ diff --git a/act/tests/data/sgpstampE31.b1.20200101.000000.nc b/act/tests/data/sgpstampE31.b1.20200101.000000.nc new file mode 100644 index 0000000000..863f10def7 Binary files /dev/null and b/act/tests/data/sgpstampE31.b1.20200101.000000.nc differ diff --git a/act/tests/data/sgpstampE32.b1.20200101.000000.nc b/act/tests/data/sgpstampE32.b1.20200101.000000.nc new file mode 100644 index 0000000000..a4619d15ad Binary files /dev/null and b/act/tests/data/sgpstampE32.b1.20200101.000000.nc differ diff --git a/act/tests/data/sgpstampE33.b1.20200101.000000.nc b/act/tests/data/sgpstampE33.b1.20200101.000000.nc new file mode 100644 index 0000000000..3db75441af Binary files /dev/null and b/act/tests/data/sgpstampE33.b1.20200101.000000.nc differ diff --git a/act/tests/data/sgpstampE34.b1.20200101.000000.nc b/act/tests/data/sgpstampE34.b1.20200101.000000.nc new file mode 100644 index 0000000000..84aa0cf4f4 Binary files /dev/null and b/act/tests/data/sgpstampE34.b1.20200101.000000.nc differ diff --git a/act/tests/data/sgpstampE9.b1.20200101.000000.nc b/act/tests/data/sgpstampE9.b1.20200101.000000.nc new file mode 100644 index 0000000000..8180d4031f Binary files /dev/null and b/act/tests/data/sgpstampE9.b1.20200101.000000.nc differ diff --git a/act/tests/sample_files.py b/act/tests/sample_files.py index f8ec87be7c..7432d99411 100644 --- a/act/tests/sample_files.py +++ b/act/tests/sample_files.py @@ -16,8 +16,10 @@ EXAMPLE_SONDE_WILDCARD EXAMPLE_MET_WILDCARD EXAMPLE_MET_CONTOUR + EXAMPLE_MET_CSV EXAMPLE_CEIL_WILDCARD EXAMPLE_TWP_SONDE_WILDCARD + EXAMPLE_TWP_SONDE_20060121 EXAMPLE_ANL_CSV EXAMPLE_VISST EXAMPLE_DLPPI @@ -31,12 +33,21 @@ EXAMPLE_SIGMA_MPLV5 EXAMPLE_CO2FLX4M EXAMPLE_SIRS + EXAMPLE_GML_RADIATION + EXAMPLE_GML_MET + EXAMPLE_GML_OZONE + EXAMPLE_GML_CO2 + EXAMPLE_GML_HALO + EXAMPLE_MET_TEST1 + EXAMPLE_MET_TEST2 + EXAMPLE_STAMP_WILDCARD """ import os DATA_PATH = os.path.join(os.path.dirname(__file__), 'data') EXAMPLE_MET1 = os.path.join(DATA_PATH, 'sgpmetE13.b1.20190101.000000.cdf') +EXAMPLE_MET_CSV = os.path.join(DATA_PATH, 'sgpmetE13.*csv') EXAMPLE_METE40 = os.path.join(DATA_PATH, 'sgpmetE40.b1.20190508.000000.cdf') EXAMPLE_CEIL1 = os.path.join(DATA_PATH, 'sgpceilC1.b1.20190101.000000.nc') EXAMPLE_SONDE1 = os.path.join(DATA_PATH, @@ -47,6 +58,7 @@ EXAMPLE_MET_CONTOUR = os.path.join(DATA_PATH, 'sgpmet*20190508*.cdf') EXAMPLE_CEIL_WILDCARD = os.path.join(DATA_PATH, 'sgpceil*.cdf') EXAMPLE_TWP_SONDE_WILDCARD = os.path.join(DATA_PATH, 'twpsondewnpn*.cdf') +EXAMPLE_TWP_SONDE_20060121 = os.path.join(DATA_PATH, 'twpsondewnpn*20060121*.cdf') EXAMPLE_ANL_CSV = os.path.join(DATA_PATH, 'anltwr_mar19met.data') EXAMPLE_VISST = os.path.join( DATA_PATH, 'twpvisstgridirtemp.c1.20050705.002500.nc') @@ -68,7 +80,15 @@ EXAMPLE_MFRSR = os.path.join(DATA_PATH, 'sgpmfrsr7nchE38.b1.20190514.180000.nc') EXAMPLE_SURFSPECALB1MLAWER = os.path.join( DATA_PATH, 'nsasurfspecalb1mlawerC1.c1.20160609.080000.nc') -EXAMPLE_SIGMA_MPLV5 = os.path.join(DATA_PATH, '201509021500.bin') +EXAMPLE_SIGMA_MPLV5 = os.path.join(DATA_PATH, '201509021500.bi') EXAMPLE_RL1 = os.path.join(DATA_PATH, 'sgprlC1.a0.20160131.000000.nc') EXAMPLE_CO2FLX4M = os.path.join(DATA_PATH, 'sgpco2flx4mC1.b1.20201007.001500.nc') EXAMPLE_SIRS = os.path.join(DATA_PATH, 'sgpsirsE13.b1.20190101.000000.cdf') +EXAMPLE_GML_RADIATION = os.path.join(DATA_PATH, 'brw21001.dat') +EXAMPLE_GML_MET = os.path.join(DATA_PATH, 'met_brw_insitu_1_obop_hour_2020.txt') +EXAMPLE_GML_OZONE = os.path.join(DATA_PATH, 'brw_12_2020_hour.dat') +EXAMPLE_GML_CO2 = os.path.join(DATA_PATH, 'co2_brw_surface-insitu_1_ccgg_MonthlyData.txt') +EXAMPLE_GML_HALO = os.path.join(DATA_PATH, 'brw_CCl4_Day.dat') +EXAMPLE_MET_TEST1 = os.path.join(DATA_PATH, 'sgpmet_no_time.nc') +EXAMPLE_MET_TEST2 = os.path.join(DATA_PATH, 'sgpmet_test_time.nc') +EXAMPLE_STAMP_WILDCARD = os.path.join(DATA_PATH, 'sgpstamp*202001*.nc') diff --git a/act/tests/test_CropScape.py b/act/tests/test_CropScape.py deleted file mode 100644 index e12b00b5f4..0000000000 --- a/act/tests/test_CropScape.py +++ /dev/null @@ -1,18 +0,0 @@ -import act -import json -import requests -from requests.packages.urllib3.exceptions import InsecureRequestWarning -requests.packages.urllib3.disable_warnings(InsecureRequestWarning) - - -def test_cropType(): - year = 2018 - lat = 37.1509 - lon = -98.362 - try: - crop = act.discovery.get_CropScape.croptype(lat, lon, year) - except Exception: - return - - if crop is not None: - assert crop == 'Grass/Pasture' diff --git a/act/tests/test_asos.py b/act/tests/test_asos.py deleted file mode 100644 index d835bba87f..0000000000 --- a/act/tests/test_asos.py +++ /dev/null @@ -1,22 +0,0 @@ -import numpy as np -from act.discovery import get_asos -from datetime import datetime - - -def test_get_ord(): - time_window = [datetime(2020, 2, 4, 2, 0), datetime(2020, 2, 12, 10, 0)] - my_asoses = get_asos(time_window, station="ORD") - assert "ORD" in my_asoses.keys() - assert np.all( - np.equal(my_asoses["ORD"]["sknt"].values[:10], - np.array([13., 11., 11., 11., 9., 10., 10., 11., 11., 11.]))) - - -def test_get_region(): - my_keys = ['MDW', 'IGQ', 'ORD', '06C', 'PWK', 'LOT', 'GYY'] - time_window = [datetime(2020, 2, 4, 2, 0), datetime(2020, 2, 12, 10, 0)] - lat_window = (41.8781 - 0.5, 41.8781 + 0.5) - lon_window = (-87.6298 - 0.5, -87.6298 + 0.5) - my_asoses = get_asos(time_window, lat_range=lat_window, lon_range=lon_window) - asos_keys = [x for x in my_asoses.keys()] - assert asos_keys == my_keys diff --git a/act/tests/test_clean.py b/act/tests/test_clean.py deleted file mode 100644 index e94ffeefee..0000000000 --- a/act/tests/test_clean.py +++ /dev/null @@ -1,64 +0,0 @@ -from act.io.armfiles import read_netcdf -from act.tests import EXAMPLE_CEIL1 - - -def test_clean(): - # Read test data - ceil_ds = read_netcdf([EXAMPLE_CEIL1]) - # Cleanup QC data - ceil_ds.clean.cleanup(clean_arm_state_vars=['detection_status']) - - # Check that global attribures are removed - global_attributes = ['qc_bit_comment', - 'qc_bit_1_description', - 'qc_bit_1_assessment', - 'qc_bit_2_description', - 'qc_bit_2_assessment' - 'qc_bit_3_description', - 'qc_bit_3_assessment' - ] - - for glb_att in global_attributes: - assert glb_att not in ceil_ds.attrs.keys() - - # Check that CF attributes are set including new flag_assessments - var_name = 'qc_first_cbh' - for attr_name in ['flag_masks', 'flag_meanings', 'flag_assessments']: - assert attr_name in ceil_ds[var_name].attrs.keys() - assert isinstance(ceil_ds[var_name].attrs[attr_name], list) - - # Check that the flag_mask values are set correctly - assert ceil_ds['qc_first_cbh'].attrs['flag_masks'] == [1, 2, 4] - - # Check that the flag_meanings values are set correctly - assert (ceil_ds['qc_first_cbh'].attrs['flag_meanings'] == - ['Value is equal to missing_value.', - 'Value is less than the fail_min.', - 'Value is greater than the fail_max.']) - - # Check the value of flag_assessments is as expected - assert ceil_ds['qc_first_cbh'].attrs['flag_assessments'] == ['Bad', 'Bad', 'Bad'] - - # Check that ancillary varibles is being added - assert 'qc_first_cbh' in ceil_ds['first_cbh'].attrs['ancillary_variables'].split() - - # Check that state field is updated to CF - assert 'flag_values' in ceil_ds['detection_status'].attrs.keys() - assert isinstance(ceil_ds['detection_status'].attrs['flag_values'], list) - assert ceil_ds['detection_status'].attrs['flag_values'] == [0, 1, 2, 3, 4, 5] - - assert 'flag_meanings' in ceil_ds['detection_status'].attrs.keys() - assert isinstance(ceil_ds['detection_status'].attrs['flag_meanings'], list) - assert (ceil_ds['detection_status'].attrs['flag_meanings'] == - ['No significant backscatter', - 'One cloud base detected', - 'Two cloud bases detected', - 'Three cloud bases detected', - 'Full obscuration determined but no cloud base detected', - 'Some obscuration detected but determined to be transparent']) - - assert 'flag_0_description' not in ceil_ds['detection_status'].attrs.keys() - assert ('detection_status' in - ceil_ds['first_cbh'].attrs['ancillary_variables'].split()) - - ceil_ds.close() diff --git a/act/tests/test_correct.py b/act/tests/test_correct.py index d7c9b72b82..979203d99d 100644 --- a/act/tests/test_correct.py +++ b/act/tests/test_correct.py @@ -35,7 +35,7 @@ def test_correct_mpl(): 181.9355]) np.testing.assert_allclose( sig_cross_pol, [-0.5823283, -1.6066532, -1.7153032, - -2.520143, -2.275405]) + -2.520143, -2.275405], rtol=4e-07) np.testing.assert_allclose( sig_co_pol, [12.5631485, 11.035495, 11.999875, 11.09393, 11.388968]) diff --git a/act/tests/test_discovery.py b/act/tests/test_discovery.py new file mode 100644 index 0000000000..7f93b55ed1 --- /dev/null +++ b/act/tests/test_discovery.py @@ -0,0 +1,79 @@ +import act +import requests +import numpy as np +import os +import glob +from datetime import datetime +from act.discovery import get_asos +from requests.packages.urllib3.exceptions import InsecureRequestWarning +requests.packages.urllib3.disable_warnings(InsecureRequestWarning) + + +def test_cropType(): + year = 2018 + lat = 37.1509 + lon = -98.362 + try: + crop = act.discovery.get_CropScape.croptype(lat, lon, year) + crop2 = act.discovery.get_CropScape.croptype(lat, lon) + except Exception: + return + + if crop is not None: + assert crop == 'Grass/Pasture' + if crop2 is not None: + assert crop2 == 'Soybeans' + + +def test_get_ord(): + time_window = [datetime(2020, 2, 4, 2, 0), datetime(2020, 2, 12, 10, 0)] + my_asoses = get_asos(time_window, station="ORD") + assert "ORD" in my_asoses.keys() + assert np.all( + np.equal(my_asoses["ORD"]["sknt"].values[:10], + np.array([13., 11., 14., 14., 13., 11., 14., 13., 13., 13.]))) + + +def test_get_region(): + my_keys = ['MDW', 'IGQ', 'ORD', '06C', 'PWK', 'LOT', 'GYY'] + time_window = [datetime(2020, 2, 4, 2, 0), datetime(2020, 2, 12, 10, 0)] + lat_window = (41.8781 - 0.5, 41.8781 + 0.5) + lon_window = (-87.6298 - 0.5, -87.6298 + 0.5) + my_asoses = get_asos(time_window, lat_range=lat_window, lon_range=lon_window) + asos_keys = [x for x in my_asoses.keys()] + assert asos_keys == my_keys + + +def test_get_armfile(): + if not os.path.isdir((os.getcwd() + '/data/')): + os.makedirs((os.getcwd() + '/data/')) + + uname = os.getenv('ARM_USERNAME') + token = os.getenv('ARM_PASSWORD') + + if uname is not None: + datastream = 'sgpmetE13.b1' + startdate = '2020-01-01' + enddate = startdate + outdir = os.getcwd() + '/data/' + + results = act.discovery.get_armfiles.download_data(uname, token, + datastream, + startdate, enddate, + output=outdir) + files = glob.glob(outdir + datastream + '*20200101*cdf') + if len(results) > 0: + assert files is not None + assert 'sgpmetE13' in files[0] + + if files is not None: + if len(files) > 0: + os.remove(files[0]) + + datastream = 'sgpmeetE13.b1' + act.discovery.get_armfiles.download_data(uname, token, + datastream, + startdate, enddate, + output=outdir) + files = glob.glob(outdir + datastream + '*20200101*cdf') + assert len(files) == 0 diff --git a/act/tests/test_io.py b/act/tests/test_io.py index 6000fd7a8c..9555d6fd35 100644 --- a/act/tests/test_io.py +++ b/act/tests/test_io.py @@ -1,15 +1,31 @@ import act +from act.io.noaagml import read_gml import act.tests.sample_files as sample_files from pathlib import Path import tempfile +import numpy as np +import glob def test_io(): - sonde_ds = act.io.armfiles.read_netcdf( - [act.tests.EXAMPLE_MET1]) + sonde_ds = act.io.armfiles.read_netcdf([act.tests.EXAMPLE_MET1]) assert 'temp_mean' in sonde_ds.variables.keys() assert 'rh_mean' in sonde_ds.variables.keys() assert sonde_ds.attrs['_arm_standards_flag'] == (1 << 0) + + with np.testing.assert_raises(OSError): + act.io.armfiles.read_netcdf([]) + + result = act.io.armfiles.read_netcdf([], return_None=True) + assert result is None + result = act.io.armfiles.read_netcdf(['./randomfile.nc'], return_None=True) + assert result is None + + obj = act.io.armfiles.read_netcdf([act.tests.EXAMPLE_MET_TEST1]) + assert 'time' in obj + + obj = act.io.armfiles.read_netcdf([act.tests.EXAMPLE_MET_TEST2]) + assert obj['time'].values[10] == np.datetime64('2019-01-01T00:10:00') sonde_ds.close() @@ -23,7 +39,7 @@ def test_io_mfdataset(): sonde_ds.close() -def test_io_anl_csv(): +def test_io_csv(): headers = ['day', 'month', 'year', 'time', 'pasquill', 'wdir_60m', 'wspd_60m', 'wdir_60m_std', 'temp_60m', 'wdir_10m', 'wspd_10m', @@ -39,15 +55,34 @@ def test_io_anl_csv(): assert anl_ds['temp_60m'].values[10] == -1.7 anl_ds.close() + files = glob.glob(act.tests.EXAMPLE_MET_CSV) + obj = act.io.csvfiles.read_csv(files[0]) + assert 'date_time' in obj + assert '_datastream' in obj.attrs + def test_io_dod(): dims = {'time': 1440, 'drop_diameter': 50} - obj = act.io.armfiles.create_obj_from_arm_dod('vdis.b1', dims, version='1.2', - scalar_fill_dim='time') - assert 'moment1' in obj - assert len(obj['base_time'].values) == 1440 - assert len(obj['drop_diameter'].values) == 50 + try: + obj = act.io.armfiles.create_obj_from_arm_dod('vdis.b1', dims, version='1.2', + scalar_fill_dim='time') + assert 'moment1' in obj + assert len(obj['base_time'].values) == 1440 + assert len(obj['drop_diameter'].values) == 50 + with np.testing.assert_warns(UserWarning): + obj2 = act.io.armfiles.create_obj_from_arm_dod('vdis.b1', dims, + scalar_fill_dim='time') + assert 'moment1' in obj2 + assert len(obj2['base_time'].values) == 1440 + assert len(obj2['drop_diameter'].values) == 50 + with np.testing.assert_raises(ValueError): + obj = act.io.armfiles.create_obj_from_arm_dod('vdis.b1', {}, version='1.2') + + except Exception: + return + obj.close() + obj2.close() def test_io_write(): @@ -61,15 +96,59 @@ def test_io_write(): if var_name not in keep_vars: del sonde_ds[var_name] sonde_ds.write.write_netcdf(path=write_file, FillValue=-9999) - sonde_ds.close() - del sonde_ds - sonde_ds = act.io.armfiles.read_netcdf(str(write_file)) - assert list(sonde_ds.data_vars) == keep_vars - assert isinstance(sonde_ds['qc_tdry'].attrs['flag_meanings'], str) - assert sonde_ds['qc_tdry'].attrs['flag_meanings'].count('__') == 21 + sonde_ds_read = act.io.armfiles.read_netcdf(str(write_file)) + assert list(sonde_ds_read.data_vars) == keep_vars + assert isinstance(sonde_ds_read['qc_tdry'].attrs['flag_meanings'], str) + assert sonde_ds_read['qc_tdry'].attrs['flag_meanings'].count('__') == 21 for attr in ['qc_standards_version', 'qc_method', 'qc_comment']: - assert attr not in list(sonde_ds.attrs) + assert attr not in list(sonde_ds_read.attrs) + sonde_ds_read.close() + del sonde_ds_read + + sonde_ds.close() + + sonde_ds = act.io.armfiles.read_netcdf(sample_files.EXAMPLE_EBBR1) + sonde_ds.clean.cleanup() + assert 'fail_min' in sonde_ds['qc_home_signal_15'].attrs + assert 'standard_name' in sonde_ds['qc_home_signal_15'].attrs + assert 'flag_masks' in sonde_ds['qc_home_signal_15'].attrs + + with tempfile.TemporaryDirectory() as tmpdirname: + cf_convention = 'CF-1.8' + write_file = Path(tmpdirname, Path(sample_files.EXAMPLE_EBBR1).name) + sonde_ds.write.write_netcdf(path=write_file, make_copy=False, join_char='_', + cf_compliant=True, cf_convention=cf_convention) + + sonde_ds_read = act.io.armfiles.read_netcdf(str(write_file)) + + assert cf_convention in sonde_ds_read.attrs['Conventions'].split() + assert sonde_ds_read.attrs['FeatureType'] == 'timeSeries' + global_att_keys = [ii for ii in sonde_ds_read.attrs.keys() if not ii.startswith('_')] + assert global_att_keys[-1] == 'history' + assert sonde_ds_read['alt'].attrs['axis'] == 'Z' + assert sonde_ds_read['alt'].attrs['positive'] == 'up' + + sonde_ds_read.close() + del sonde_ds_read + + sonde_ds.close() + + obj = act.io.armfiles.read_netcdf(sample_files.EXAMPLE_CEIL1) + with tempfile.TemporaryDirectory() as tmpdirname: + cf_convention = 'CF-1.8' + write_file = Path(tmpdirname, Path(sample_files.EXAMPLE_CEIL1).name) + obj.write.write_netcdf(path=write_file, make_copy=False, join_char='_', + cf_compliant=True, cf_convention=cf_convention) + + obj_read = act.io.armfiles.read_netcdf(str(write_file)) + + assert cf_convention in obj_read.attrs['Conventions'].split() + assert obj_read.attrs['FeatureType'] == 'timeSeriesProfile' + assert len(obj_read.dims) > 1 + + obj_read.close() + del obj_read def test_io_mpldataset(): @@ -94,3 +173,97 @@ def test_io_mpldataset(): # Tests attributes assert '_datastream' in mpl_ds.attrs.keys() mpl_ds.close() + + +def test_read_gml(): + # Test Radiation + ds = read_gml(sample_files.EXAMPLE_GML_RADIATION, datatype='RADIATION') + assert np.isclose(np.nansum(ds['solar_zenith_angle']), 1629.68) + assert np.isclose(np.nansum(ds['upwelling_infrared_case_temp']), 4185.73) + assert (ds['upwelling_infrared_case_temp'].attrs['ancillary_variables'] == + 'qc_upwelling_infrared_case_temp') + assert ds['qc_upwelling_infrared_case_temp'].attrs['flag_values'] == [0, 1, 2] + assert (ds['qc_upwelling_infrared_case_temp'].attrs['flag_meanings'] == + ['Not failing any tests', 'Knowingly bad value', 'Should be used with scrutiny']) + assert (ds['qc_upwelling_infrared_case_temp'].attrs['flag_assessments'] == + ['Good', 'Bad', 'Indeterminate']) + assert ds['time'].values[-1] == np.datetime64('2021-01-01T00:17:00') + + ds = read_gml(sample_files.EXAMPLE_GML_RADIATION, convert_missing=False) + assert np.isclose(np.nansum(ds['solar_zenith_angle']), 1629.68) + assert np.isclose(np.nansum(ds['upwelling_infrared_case_temp']), 4185.73) + assert (ds['upwelling_infrared_case_temp'].attrs['ancillary_variables'] == + 'qc_upwelling_infrared_case_temp') + assert ds['qc_upwelling_infrared_case_temp'].attrs['flag_values'] == [0, 1, 2] + assert (ds['qc_upwelling_infrared_case_temp'].attrs['flag_meanings'] == + ['Not failing any tests', 'Knowingly bad value', 'Should be used with scrutiny']) + assert (ds['qc_upwelling_infrared_case_temp'].attrs['flag_assessments'] == + ['Good', 'Bad', 'Indeterminate']) + assert ds['time'].values[-1] == np.datetime64('2021-01-01T00:17:00') + + # Test MET + ds = read_gml(sample_files.EXAMPLE_GML_MET, datatype='MET') + assert np.isclose(np.nansum(ds['wind_speed'].values), 140.999) + assert ds['wind_speed'].attrs['units'] == 'm/s' + assert np.isnan(ds['wind_speed'].attrs['_FillValue']) + assert np.sum(np.isnan(ds['preciptation_intensity'].values)) == 19 + assert ds['preciptation_intensity'].attrs['units'] == 'mm/hour' + assert ds['time'].values[0] == np.datetime64('2020-01-01T01:00:00') + + ds = read_gml(sample_files.EXAMPLE_GML_MET, convert_missing=False) + assert np.isclose(np.nansum(ds['wind_speed'].values), 140.999) + assert ds['wind_speed'].attrs['units'] == 'm/s' + assert np.isclose(ds['wind_speed'].attrs['_FillValue'], -999.9) + assert np.sum(ds['preciptation_intensity'].values) == -1881 + assert ds['preciptation_intensity'].attrs['units'] == 'mm/hour' + assert ds['time'].values[0] == np.datetime64('2020-01-01T01:00:00') + + # Test Ozone + ds = read_gml(sample_files.EXAMPLE_GML_OZONE, datatype='OZONE') + assert np.isclose(np.nansum(ds['ozone'].values), 582.76) + assert ds['ozone'].attrs['long_name'] == 'Ozone' + assert ds['ozone'].attrs['units'] == 'ppb' + assert np.isnan(ds['ozone'].attrs['_FillValue']) + assert ds['time'].values[0] == np.datetime64('2020-12-01T00:00:00') + + ds = read_gml(sample_files.EXAMPLE_GML_OZONE) + assert np.isclose(np.nansum(ds['ozone'].values), 582.76) + assert ds['ozone'].attrs['long_name'] == 'Ozone' + assert ds['ozone'].attrs['units'] == 'ppb' + assert np.isnan(ds['ozone'].attrs['_FillValue']) + assert ds['time'].values[0] == np.datetime64('2020-12-01T00:00:00') + + # Test Carbon Dioxide + ds = read_gml(sample_files.EXAMPLE_GML_CO2, datatype='co2') + assert np.isclose(np.nansum(ds['co2'].values), 2307.630) + assert (ds['qc_co2'].values == + np.array([1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0], dtype=int)).all() + assert ds['co2'].attrs['units'] == 'ppm' + assert np.isnan(ds['co2'].attrs['_FillValue']) + assert ds['qc_co2'].attrs['flag_assessments'] == ['Bad', 'Indeterminate'] + assert ds['latitude'].attrs['standard_name'] == 'latitude' + + ds = read_gml(sample_files.EXAMPLE_GML_CO2, convert_missing=False) + assert np.isclose(np.nansum(ds['co2'].values), -3692.3098) + assert ds['co2'].attrs['_FillValue'] == -999.99 + assert (ds['qc_co2'].values == + np.array([1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0], dtype=int)).all() + assert ds['co2'].attrs['units'] == 'ppm' + assert np.isclose(ds['co2'].attrs['_FillValue'], -999.99) + assert ds['qc_co2'].attrs['flag_assessments'] == ['Bad', 'Indeterminate'] + assert ds['latitude'].attrs['standard_name'] == 'latitude' + + # Test Halocarbon + ds = read_gml(sample_files.EXAMPLE_GML_HALO, datatype='HALO') + assert np.isclose(np.nansum(ds['CCl4'].values), 1342.6499) + assert ds['CCl4'].attrs['units'] == 'ppt' + assert ds['CCl4'].attrs['long_name'] == 'Carbon Tetrachloride (CCl4) daily median' + assert np.isnan(ds['CCl4'].attrs['_FillValue']) + assert ds['time'].values[0] == np.datetime64('1998-06-16T00:00:00') + + ds = read_gml(sample_files.EXAMPLE_GML_HALO) + assert np.isclose(np.nansum(ds['CCl4'].values), 1342.6499) + assert ds['CCl4'].attrs['units'] == 'ppt' + assert ds['CCl4'].attrs['long_name'] == 'Carbon Tetrachloride (CCl4) daily median' + assert np.isnan(ds['CCl4'].attrs['_FillValue']) + assert ds['time'].values[0] == np.datetime64('1998-06-16T00:00:00') diff --git a/act/tests/test_plotting.py b/act/tests/test_plotting.py index 267b912e8f..f1a20f0b4e 100644 --- a/act/tests/test_plotting.py +++ b/act/tests/test_plotting.py @@ -1,21 +1,24 @@ -import act.io.armfiles as arm -import act.tests.sample_files as sample_files -import act.corrections.ceil as ceil import pytest -import os -import boto3 import numpy as np import glob import xarray as xr - +import pandas as pd +import os +import act +import act.io.armfiles as arm +import act.tests.sample_files as sample_files from act.plotting import TimeSeriesDisplay, WindRoseDisplay from act.plotting import SkewTDisplay, XSectionDisplay from act.plotting import GeographicPlotDisplay, HistogramDisplay from act.plotting import ContourDisplay from act.utils.data_utils import accumulate_precip -from botocore.handlers import disable_signing import matplotlib matplotlib.use('Agg') +try: + import metpy.calc as mpcalc + METPY = True +except ImportError: + METPY = False @pytest.mark.mpl_image_compare(tolerance=30) @@ -46,61 +49,195 @@ def test_plot(): return display.fig -@pytest.mark.mpl_image_compare(tolerance=30) -def test_multidataset_plot_tuple(): - conn = boto3.resource('s3') - conn.meta.client.meta.events.register('choose-signer.s3.*', - disable_signing) - bucket = conn.Bucket('act-tests') - if not os.path.isdir((os.getcwd() + '/data/')): - os.makedirs((os.getcwd() + '/data/')) +def test_errors(): + files = sample_files.EXAMPLE_MET_WILDCARD + obj = arm.read_netcdf(files) - for item in bucket.objects.all(): - bucket.download_file(item.key, (os.getcwd() + '/data/' + item.key)) + display = TimeSeriesDisplay(obj) + display.axes = None + with np.testing.assert_raises(RuntimeError): + display.day_night_background() + + display = TimeSeriesDisplay({'met': obj, 'met2': obj}) + with np.testing.assert_raises(ValueError): + display.plot('temp_mean') + with np.testing.assert_raises(ValueError): + display.qc_flag_block_plot('qc_temp_mean') + with np.testing.assert_raises(ValueError): + display.plot_barbs_from_spd_dir('wdir_vec_mean', 'wspd_vec_mean') + with np.testing.assert_raises(ValueError): + display.plot_barbs_from_u_v('wdir_vec_mean', 'wspd_vec_mean') + + del obj.attrs['_file_dates'] + + data = np.empty(len(obj['time'])) * np.nan + lat = obj['lat'].values + lon = obj['lon'].values + obj['lat'].values = data + obj['lon'].values = data - ceil_ds = arm.read_netcdf('data/sgpceilC1.b1*') - sonde_ds = arm.read_netcdf( - sample_files.EXAMPLE_MET_WILDCARD) - # Removing fill value of -9999 as it was causing some warnings - ceil_ds = ceil.correct_ceil(ceil_ds) + display = TimeSeriesDisplay(obj) + display.plot('temp_mean') + display.set_yrng([0, 0]) + with np.testing.assert_warns(RuntimeWarning): + display.day_night_background() + obj['lat'].values = lat + with np.testing.assert_warns(RuntimeWarning): + display.day_night_background() + obj['lon'].values = lon * 100. + with np.testing.assert_warns(RuntimeWarning): + display.day_night_background() + obj['lat'].values = lat * 100. + with np.testing.assert_warns(RuntimeWarning): + display.day_night_background() + + obj.close() + + # Test some of the other errors + obj = arm.read_netcdf(files) + del obj['temp_mean'].attrs['units'] + display = TimeSeriesDisplay(obj) + display.axes = None + with np.testing.assert_raises(RuntimeError): + display.set_yrng([0, 10]) + with np.testing.assert_raises(RuntimeError): + display.set_xrng([0, 10]) + display.fig = None + display.plot('temp_mean', add_nan=True) + + assert display.fig is not None + assert display.axes is not None + + with np.testing.assert_raises(AttributeError): + display = TimeSeriesDisplay([]) + + fig, ax = matplotlib.pyplot.subplots() + display = TimeSeriesDisplay(obj) + display.add_subplots((2, 2), figsize=(15, 10)) + display.assign_to_figure_axis(fig, ax) + assert display.fig is not None + assert display.axes is not None + + obj = arm.read_netcdf(files) + display = TimeSeriesDisplay(obj) + obj.clean.cleanup() + display.axes = None + display.fig = None + display.qc_flag_block_plot('atmos_pressure') + assert display.fig is not None + assert display.axes is not None + + +def test_histogram_errors(): + files = sample_files.EXAMPLE_MET1 + obj = arm.read_netcdf(files) + + histdisplay = HistogramDisplay(obj) + histdisplay.axes = None + with np.testing.assert_raises(RuntimeError): + histdisplay.set_yrng([0, 10]) + with np.testing.assert_raises(RuntimeError): + histdisplay.set_xrng([-40, 40]) + histdisplay.fig = None + histdisplay.plot_stacked_bar_graph('temp_mean', bins=np.arange(-40, 40, 5)) + histdisplay.set_yrng([0, 0]) + assert histdisplay.yrng[0][1] == 1. + assert histdisplay.fig is not None + assert histdisplay.axes is not None + + with np.testing.assert_raises(AttributeError): + HistogramDisplay([]) + + histdisplay.axes = None + histdisplay.fig = None + histdisplay.plot_stairstep_graph('temp_mean', bins=np.arange(-40, 40, 5)) + assert histdisplay.fig is not None + assert histdisplay.axes is not None + + sigma = 10 + mu = 50 + bins = np.linspace(0, 100, 50) + ydata = 1 / (sigma * np.sqrt(2 * np.pi)) * np.exp(-(bins - mu)**2 / (2 * sigma**2)) + y_array = xr.DataArray(ydata, dims={'bins': bins}) + bins = xr.DataArray(bins, dims={'bins': bins}) + my_fake_ds = xr.Dataset({'bins': bins, 'ydata': y_array}) + histdisplay = HistogramDisplay(my_fake_ds) + histdisplay.axes = None + histdisplay.fig = None + histdisplay.plot_size_distribution('ydata', 'bins', set_title='Fake distribution.') + assert histdisplay.fig is not None + assert histdisplay.axes is not None + + sonde_ds = arm.read_netcdf(sample_files.EXAMPLE_SONDE1) + histdisplay = HistogramDisplay({'sgpsondewnpnC1.b1': sonde_ds}) + histdisplay.axes = None + histdisplay.fig = None + histdisplay.plot_heatmap('tdry', 'alt', x_bins=np.arange(-60, 10, 1), + y_bins=np.linspace(0, 10000., 50), cmap='coolwarm') + assert histdisplay.fig is not None + assert histdisplay.axes is not None + + +def test_xsection_errors(): + obj = arm.read_netcdf(sample_files.EXAMPLE_CEIL1) + + display = XSectionDisplay(obj, figsize=(10, 8), subplot_shape=(2,)) + display.axes = None + with np.testing.assert_raises(RuntimeError): + display.set_yrng([0, 10]) + with np.testing.assert_raises(RuntimeError): + display.set_xrng([-40, 40]) + + display = XSectionDisplay(obj, figsize=(10, 8), subplot_shape=(1,)) + with np.testing.assert_raises(RuntimeError): + display.plot_xsection(None, 'backscatter', x='time') + + obj.close() + + +@pytest.mark.mpl_image_compare(tolerance=30) +def test_multidataset_plot_tuple(): + obj = arm.read_netcdf(sample_files.EXAMPLE_MET1) + obj2 = arm.read_netcdf(sample_files.EXAMPLE_SIRS) + obj = obj.rename({'lat': 'fun_time'}) + obj['fun_time'].attrs['standard_name'] = 'latitude' + obj = obj.rename({'lon': 'not_so_fun_time'}) + obj['not_so_fun_time'].attrs['standard_name'] = 'longitude' # You can use tuples if the datasets in the tuple contain a # datastream attribute. This is required in all ARM datasets. display = TimeSeriesDisplay( - (ceil_ds, sonde_ds), subplot_shape=(2,), figsize=(15, 10)) - display.plot('backscatter', 'sgpceilC1.b1', subplot_index=(0,)) + (obj, obj2), subplot_shape=(2,), figsize=(15, 10)) + display.plot('short_direct_normal', 'sgpsirsE13.b1', subplot_index=(0,)) + display.day_night_background('sgpsirsE13.b1', subplot_index=(0,)) display.plot('temp_mean', 'sgpmetE13.b1', subplot_index=(1,)) display.day_night_background('sgpmetE13.b1', subplot_index=(1,)) - ceil_ds.close() - sonde_ds.close() + + ax = act.plotting.common.parse_ax(ax=None) + ax, fig = act.plotting.common.parse_ax_fig(ax=None, fig=None) + obj.close() + obj2.close() return display.fig @pytest.mark.mpl_image_compare(tolerance=30) def test_multidataset_plot_dict(): - conn = boto3.resource('s3') - conn.meta.client.meta.events.register('choose-signer.s3.*', - disable_signing) - bucket = conn.Bucket('act-tests') - if not os.path.isdir((os.getcwd() + '/data/')): - os.makedirs((os.getcwd() + '/data/')) - for item in bucket.objects.all(): - bucket.download_file(item.key, (os.getcwd() + '/data/' + item.key)) - - ceil_ds = arm.read_netcdf('data/sgpceilC1.b1*') - sonde_ds = arm.read_netcdf( - sample_files.EXAMPLE_MET_WILDCARD) - ceil_ds = ceil.correct_ceil(ceil_ds, fill_value=-9999.) + obj = arm.read_netcdf(sample_files.EXAMPLE_MET1) + obj2 = arm.read_netcdf(sample_files.EXAMPLE_SIRS) + # You can use tuples if the datasets in the tuple contain a + # datastream attribute. This is required in all ARM datasets. display = TimeSeriesDisplay( - {'ceiliometer': ceil_ds, 'rawinsonde': sonde_ds}, + {'sirs': obj2, 'met': obj}, subplot_shape=(2,), figsize=(15, 10)) - display.plot('backscatter', 'ceiliometer', subplot_index=(0,)) - display.plot('temp_mean', 'rawinsonde', subplot_index=(1,)) - display.day_night_background('rawinsonde', subplot_index=(1,)) - ceil_ds.close() - sonde_ds.close() + display.plot('short_direct_normal', 'sirs', subplot_index=(0,)) + display.day_night_background('sirs', subplot_index=(0,)) + display.plot('temp_mean', 'met', subplot_index=(1,)) + display.day_night_background('met', subplot_index=(1,)) + obj.close() + obj2.close() + return display.fig @@ -112,7 +249,10 @@ def test_wind_rose(): WindDisplay = WindRoseDisplay(sonde_ds, figsize=(10, 10)) WindDisplay.plot('deg', 'wspd', spd_bins=np.linspace(0, 20, 10), num_dirs=30, - tick_interval=2) + tick_interval=2, cmap='viridis') + WindDisplay.set_thetarng(trng=(0., 360.)) + WindDisplay.set_rrng((0., 14)) + sonde_ds.close() return WindDisplay.fig @@ -137,30 +277,52 @@ def test_skewt_plot(): sonde_ds = arm.read_netcdf( sample_files.EXAMPLE_SONDE1) - try: + if METPY: skewt = SkewTDisplay(sonde_ds) skewt.plot_from_u_and_v('u_wind', 'v_wind', 'pres', 'tdry', 'dp') sonde_ds.close() return skewt.fig - except Exception: - pass @pytest.mark.mpl_image_compare(tolerance=30) def test_skewt_plot_spd_dir(): - sonde_ds = arm.read_netcdf( - sample_files.EXAMPLE_SONDE1) + sonde_ds = arm.read_netcdf(sample_files.EXAMPLE_SONDE1) - try: - skewt = SkewTDisplay(sonde_ds) + if METPY: + skewt = SkewTDisplay(sonde_ds, ds_name='act_datastream') skewt.plot_from_spd_and_dir('wspd', 'deg', 'pres', 'tdry', 'dp') sonde_ds.close() return skewt.fig - except Exception: - pass -@pytest.mark.mpl_image_compare(tolerance=30) +@pytest.mark.mpl_image_compare(tolerance=67) +def test_multi_skewt_plot(): + + files = glob.glob(sample_files.EXAMPLE_TWP_SONDE_20060121) + test = {} + for f in files: + time = f.split('.')[-3] + sonde_ds = arm.read_netcdf(f) + test.update({time: sonde_ds}) + + if METPY: + skewt = SkewTDisplay(test, subplot_shape=(2, 2)) + i = 0 + j = 0 + for f in files: + time = f.split('.')[-3] + skewt.plot_from_spd_and_dir('wspd', 'deg', 'pres', 'tdry', 'dp', + subplot_index=(j, i), dsname=time, + p_levels_to_plot=np.arange(10., 1000., 25)) + if j == 1: + i += 1 + j = 0 + elif j == 0: + j += 1 + return skewt.fig + + +@pytest.mark.mpl_image_compare(tolerance=31) def test_xsection_plot(): visst_ds = arm.read_netcdf( sample_files.EXAMPLE_CEIL1) @@ -176,11 +338,14 @@ def test_xsection_plot(): def test_xsection_plot_map(): radar_ds = arm.read_netcdf( sample_files.EXAMPLE_VISST, combine='nested') - xsection = XSectionDisplay(radar_ds, figsize=(15, 8)) - xsection.plot_xsection_map(None, 'ir_temperature', vmin=220, vmax=300, cmap='Greys', - x='longitude', y='latitude', isel_kwargs={'time': 0}) - radar_ds.close() - return xsection.fig + try: + xsection = XSectionDisplay(radar_ds, figsize=(15, 8)) + xsection.plot_xsection_map(None, 'ir_temperature', vmin=220, vmax=300, cmap='Greys', + x='longitude', y='latitude', isel_kwargs={'time': 0}) + radar_ds.close() + return xsection.fig + except Exception: + pass @pytest.mark.mpl_image_compare(tolerance=30) @@ -189,7 +354,9 @@ def test_geoplot(): sample_files.EXAMPLE_SONDE1) try: geodisplay = GeographicPlotDisplay({'sgpsondewnpnC1.b1': sonde_ds}) - geodisplay.geoplot('tdry', marker='.') + geodisplay.geoplot('tdry', marker='.', cartopy_feature=['STATES', 'LAND', 'OCEAN', 'COASTLINE', + 'BORDERS', 'LAKES', 'RIVERS'], + text={'Ponca City': [-97.0725, 36.7125]}) return geodisplay.fig except Exception: pass @@ -234,6 +401,20 @@ def test_stacked_bar_graph(): return histdisplay.fig +@pytest.mark.mpl_image_compare(tolerance=30) +def test_stacked_bar_graph2(): + sonde_ds = arm.read_netcdf( + sample_files.EXAMPLE_SONDE1) + + histdisplay = HistogramDisplay({'sgpsondewnpnC1.b1': sonde_ds}) + histdisplay.plot_stacked_bar_graph('tdry') + histdisplay.set_yrng([0, 400]) + histdisplay.set_xrng([-70, 0]) + sonde_ds.close() + + return histdisplay.fig + + @pytest.mark.mpl_image_compare(tolerance=30) def test_stacked_bar_graph_sorted(): sonde_ds = arm.read_netcdf( @@ -282,13 +463,114 @@ def test_contour(): time = '2019-05-08T04:00:00.000000000' data = {} fields = {} + wind_fields = {} + station_fields = {} + for f in files: + obj = arm.read_netcdf(f) + data.update({f: obj}) + fields.update({f: ['lon', 'lat', 'temp_mean']}) + wind_fields.update({f: ['lon', 'lat', 'wspd_vec_mean', 'wdir_vec_mean']}) + station_fields.update({f: ['lon', 'lat', 'atmos_pressure']}) + + display = ContourDisplay(data, figsize=(8, 8)) + display.create_contour(fields=fields, time=time, levels=50, + contour='contour', cmap='viridis') + display.plot_vectors_from_spd_dir(fields=wind_fields, time=time, mesh=True, grid_delta=(0.1, 0.1)) + display.plot_station(fields=station_fields, time=time, markersize=7, color='red') + + return display.fig + + +@pytest.mark.mpl_image_compare(tolerance=30) +def test_contour_stamp(): + files = glob.glob(sample_files.EXAMPLE_STAMP_WILDCARD) + test = {} + stamp_fields = {} + time = '2020-01-01T00:00:00.000000000' + for f in files: + ds = f.split('/')[-1] + obj = act.io.armfiles.read_netcdf(f) + test.update({ds: obj}) + stamp_fields.update({ds: ['lon', 'lat', 'plant_water_availability_east']}) + obj.close() + + display = act.plotting.ContourDisplay(test, figsize=(8, 8)) + display.create_contour(fields=stamp_fields, time=time, levels=50, + alpha=0.5, twod_dim_value=5) + + return display.fig + + +@pytest.mark.mpl_image_compare(tolerance=30) +def test_contour2(): + files = glob.glob(sample_files.EXAMPLE_MET_CONTOUR) + time = '2019-05-08T04:00:00.000000000' + data = {} + fields = {} + wind_fields = {} + station_fields = {} + for f in files: + obj = arm.read_netcdf(f) + data.update({f: obj}) + fields.update({f: ['lon', 'lat', 'temp_mean']}) + wind_fields.update({f: ['lon', 'lat', 'wspd_vec_mean', 'wdir_vec_mean']}) + station_fields.update({f: ['lon', 'lat', 'atmos_pressure']}) + + display = ContourDisplay(data, figsize=(8, 8)) + display.create_contour(fields=fields, time=time, levels=50, + contour='contour', cmap='viridis') + display.plot_vectors_from_spd_dir(fields=wind_fields, time=time, mesh=False, grid_delta=(0.1, 0.1)) + display.plot_station(fields=station_fields, time=time, markersize=7, color='pink') + + return display.fig + + +@pytest.mark.mpl_image_compare(tolerance=30) +def test_contourf(): + files = glob.glob(sample_files.EXAMPLE_MET_CONTOUR) + time = '2019-05-08T04:00:00.000000000' + data = {} + fields = {} + wind_fields = {} + station_fields = {} + for f in files: + obj = arm.read_netcdf(f) + data.update({f: obj}) + fields.update({f: ['lon', 'lat', 'temp_mean']}) + wind_fields.update({f: ['lon', 'lat', 'wspd_vec_mean', 'wdir_vec_mean']}) + station_fields.update({f: ['lon', 'lat', 'atmos_pressure', 'temp_mean', 'rh_mean', + 'vapor_pressure_mean', 'temp_std']}) + + display = ContourDisplay(data, figsize=(8, 8)) + display.create_contour(fields=fields, time=time, levels=50, + contour='contourf', cmap='viridis') + display.plot_vectors_from_spd_dir(fields=wind_fields, time=time, mesh=True, grid_delta=(0.1, 0.1)) + display.plot_station(fields=station_fields, time=time, markersize=7, color='red') + + return display.fig + + +@pytest.mark.mpl_image_compare(tolerance=30) +def test_contourf2(): + files = glob.glob(sample_files.EXAMPLE_MET_CONTOUR) + time = '2019-05-08T04:00:00.000000000' + data = {} + fields = {} + wind_fields = {} + station_fields = {} for f in files: obj = arm.read_netcdf(f) data.update({f: obj}) fields.update({f: ['lon', 'lat', 'temp_mean']}) + wind_fields.update({f: ['lon', 'lat', 'wspd_vec_mean', 'wdir_vec_mean']}) + station_fields.update({f: ['lon', 'lat', 'atmos_pressure', 'temp_mean', 'rh_mean', + 'vapor_pressure_mean', 'temp_std']}) display = ContourDisplay(data, figsize=(8, 8)) - display.create_contour(fields=fields, time=time, levels=50) + display.create_contour(fields=fields, time=time, levels=50, + contour='contourf', cmap='viridis') + display.plot_vectors_from_spd_dir(fields=wind_fields, time=time, mesh=False, grid_delta=(0.1, 0.1)) + display.plot_station(fields=station_fields, time=time, markersize=7, color='pink') return display.fig @@ -303,7 +585,6 @@ def test_time_height_scatter(): display.time_height_scatter('tdry', day_night_background=False) sonde_ds.close() - del sonde_ds return display.fig @@ -315,14 +596,22 @@ def test_qc_bar_plot(): var_name = 'temp_mean' ds_object.qcfilter.set_test(var_name, index=range(100, 600), test_number=2) + # Testing out when the assessment is not listed + ds_object.qcfilter.set_test(var_name, index=range(500, 800), test_number=4) + ds_object['qc_' + var_name].attrs['flag_assessments'][3] = 'Wonky' + display = TimeSeriesDisplay({'sgpmetE13.b1': ds_object}, subplot_shape=(2, ), figsize=(7, 4)) display.plot(var_name, subplot_index=(0, ), assessment_overplot=True) display.day_night_background('sgpmetE13.b1', subplot_index=(0, )) - display.qc_flag_block_plot(var_name, subplot_index=(1, )) + color_lookup = {'Bad': 'red', 'Incorrect': 'red', + 'Indeterminate': 'orange', 'Suspect': 'orange', + 'Missing': 'darkgray', 'Not Failing': 'green', + 'Acceptable': 'green'} + display.qc_flag_block_plot(var_name, subplot_index=(1, ), + assessment_color=color_lookup) ds_object.close() - del ds_object return display.fig @@ -413,3 +702,45 @@ def test_assessment_overplot_multi(): ds.close() return display.fig + + +@pytest.mark.mpl_image_compare(tolerance=30) +def test_plot_barbs_from_u_v(): + sonde_ds = arm.read_netcdf( + sample_files.EXAMPLE_TWP_SONDE_WILDCARD) + BarbDisplay = TimeSeriesDisplay({'sonde_darwin': sonde_ds}) + BarbDisplay.plot_barbs_from_u_v('u_wind', 'v_wind', 'pres', + num_barbs_x=20) + sonde_ds.close() + return BarbDisplay.fig + + +@pytest.mark.mpl_image_compare(tolerance=30) +def test_plot_barbs_from_u_v2(): + bins = list(np.linspace(0, 1, 10)) + xbins = list(pd.date_range(pd.to_datetime('2020-01-01'), + pd.to_datetime('2020-01-02'), 12)) + y_data = np.full([len(xbins), len(bins)], 1.) + x_data = np.full([len(xbins), len(bins)], 2.) + y_array = xr.DataArray(y_data, dims={'xbins': xbins, 'ybins': bins}, + attrs={'units': 'm/s'}) + x_array = xr.DataArray(x_data, dims={'xbins': xbins, 'ybins': bins}, + attrs={'units': 'm/s'}) + xbins = xr.DataArray(xbins, dims={'xbins': xbins}) + ybins = xr.DataArray(bins, dims={'ybins': bins}) + fake_obj = xr.Dataset({'xbins': xbins, 'ybins': ybins, + 'ydata': y_array, 'xdata': x_array}) + BarbDisplay = TimeSeriesDisplay(fake_obj) + BarbDisplay.plot_barbs_from_u_v('xdata', 'ydata', None, num_barbs_x=20, + num_barbs_y=20, set_title='test plot', cmap='jet') + fake_obj.close() + return BarbDisplay.fig + + +@pytest.mark.mpl_image_compare(tolerance=30) +def test_2D_timeseries_plot(): + obj = arm.read_netcdf(sample_files.EXAMPLE_CEIL1) + display = TimeSeriesDisplay(obj) + display.plot('backscatter', y_rng=[0, 5000], use_var_for_y='range') + matplotlib.pyplot.show() + return display.fig diff --git a/act/tests/test_qc.py b/act/tests/test_qc.py index 7258742f0f..6ab2ad20c3 100644 --- a/act/tests/test_qc.py +++ b/act/tests/test_qc.py @@ -1,11 +1,12 @@ from act.io.armfiles import read_netcdf -from act.tests import (EXAMPLE_IRT25m20s, EXAMPLE_METE40, +from act.tests import (EXAMPLE_IRT25m20s, EXAMPLE_METE40, EXAMPLE_CEIL1, EXAMPLE_MFRSR, EXAMPLE_MET1, EXAMPLE_CO2FLX4M) from act.qc.arm import add_dqr_to_qc from act.qc.radiometer_tests import fft_shading_test from act.qc.qcfilter import parse_bit, set_bit, unset_bit import numpy as np import pytest +import copy def test_fft_shading_test(): @@ -16,6 +17,17 @@ def test_fft_shading_test(): assert np.nansum(qc_data.values) == 456 +def test_qc_test_errors(): + ds_object = read_netcdf(EXAMPLE_MET1) + var_name = 'temp_mean' + + assert ds_object.qcfilter.add_less_test(var_name, None) is None + assert ds_object.qcfilter.add_greater_test(var_name, None) is None + assert ds_object.qcfilter.add_less_equal_test(var_name, None) is None + assert ds_object.qcfilter.add_equal_to_test(var_name, None) is None + assert ds_object.qcfilter.add_not_equal_to_test(var_name, None) is None + + def test_arm_qc(): # Test DQR Webservice using known DQR variable = 'wspd_vec_mean' @@ -27,6 +39,14 @@ def test_arm_qc(): try: obj = add_dqr_to_qc(obj, variable=variable) ran = True + obj.attrs['_datastream'] = obj.attrs['datastream'] + del obj.attrs['datastream'] + obj2 = add_dqr_to_qc(obj, variable=variable) + with np.testing.assert_raises(ValueError): + del obj.attrs['_datastream'] + add_dqr_to_qc(obj, variable=variable) + + obj4 = add_dqr_to_qc(obj) except ValueError: ran = False @@ -34,10 +54,21 @@ def test_arm_qc(): assert qc_variable in obj dqr = [True for d in obj[qc_variable].attrs['flag_meanings'] if 'D190529.4' in d] assert dqr[0] is True - assert 'Suspect' not in obj[qc_variable].attrs['flag_assessments'] assert 'Incorrect' not in obj[qc_variable].attrs['flag_assessments'] + assert qc_variable in obj2 + dqr = [True for d in obj2[qc_variable].attrs['flag_meanings'] if 'D190529.4' in d] + assert dqr[0] is True + assert 'Suspect' not in obj2[qc_variable].attrs['flag_assessments'] + assert 'Incorrect' not in obj2[qc_variable].attrs['flag_assessments'] + + assert qc_variable in obj4 + dqr = [True for d in obj4[qc_variable].attrs['flag_meanings'] if 'D190529.4' in d] + assert dqr[0] is True + assert 'Suspect' not in obj4[qc_variable].attrs['flag_assessments'] + assert 'Incorrect' not in obj4[qc_variable].attrs['flag_assessments'] + def test_qcfilter(): ds_object = read_netcdf(EXAMPLE_IRT25m20s) @@ -165,7 +196,7 @@ def test_qcfilter(): assert 'flag_masks' not in ds_object[qc_var_name].attrs.keys() del ds_object[qc_var_name] - ds_object.qcfilter.add_missing_value_test(var_name, flag_value=True) + ds_object.qcfilter.add_missing_value_test(var_name, flag_value=True, prepend_text='arm') ds_object.qcfilter.add_test(var_name, index=list(range(0, 20)), test_number=2, test_meaning='Testing flag', flag_value=True, test_assessment='Suspect') @@ -194,79 +225,109 @@ def test_qctests(): # less than min test limit_value = 6.8 - result = ds_object.qcfilter.add_less_test(var_name, limit_value) + result = ds_object.qcfilter.add_less_test(var_name, limit_value, prepend_text='arm') data = ds_object.qcfilter.get_masked_data(var_name, rm_tests=result['test_number']) + assert 'arm' in result['test_meaning'] assert np.ma.count_masked(data) == 54 assert 'fail_min' in ds_object[result['qc_variable_name']].attrs.keys() assert (ds_object[result['qc_variable_name']].attrs['fail_min'].dtype == ds_object[result['variable_name']].values.dtype) assert np.isclose(ds_object[result['qc_variable_name']].attrs['fail_min'], limit_value) + result = ds_object.qcfilter.add_less_test(var_name, limit_value, test_assessment='Suspect') + assert 'warn_min' in ds_object[result['qc_variable_name']].attrs.keys() + # greator than max test limit_value = 12.7 - result = ds_object.qcfilter.add_greater_test(var_name, limit_value) + result = ds_object.qcfilter.add_greater_test(var_name, limit_value, prepend_text='arm') data = ds_object.qcfilter.get_masked_data(var_name, rm_tests=result['test_number']) + assert 'arm' in result['test_meaning'] assert np.ma.count_masked(data) == 61 assert 'fail_max' in ds_object[result['qc_variable_name']].attrs.keys() assert (ds_object[result['qc_variable_name']].attrs['fail_max'].dtype == ds_object[result['variable_name']].values.dtype) assert np.isclose(ds_object[result['qc_variable_name']].attrs['fail_max'], limit_value) + result = ds_object.qcfilter.add_greater_test(var_name, limit_value, test_assessment='Suspect') + assert 'warn_max' in ds_object[result['qc_variable_name']].attrs.keys() + # less than or equal test limit_value = 6.9 result = ds_object.qcfilter.add_less_equal_test(var_name, limit_value, - test_assessment='Suspect') + test_assessment='Suspect', + prepend_text='arm') data = ds_object.qcfilter.get_masked_data(var_name, rm_tests=result['test_number']) + assert 'arm' in result['test_meaning'] assert np.ma.count_masked(data) == 149 assert 'warn_min' in ds_object[result['qc_variable_name']].attrs.keys() assert (ds_object[result['qc_variable_name']].attrs['warn_min'].dtype == ds_object[result['variable_name']].values.dtype) assert np.isclose(ds_object[result['qc_variable_name']].attrs['warn_min'], limit_value) + result = ds_object.qcfilter.add_less_equal_test(var_name, limit_value) + assert 'fail_min' in ds_object[result['qc_variable_name']].attrs.keys() + # greater than or equal test limit_value = 12 result = ds_object.qcfilter.add_greater_equal_test(var_name, limit_value, - test_assessment='Suspect') + test_assessment='Suspect', + prepend_text='arm') data = ds_object.qcfilter.get_masked_data(var_name, rm_tests=result['test_number']) + assert 'arm' in result['test_meaning'] assert np.ma.count_masked(data) == 606 assert 'warn_max' in ds_object[result['qc_variable_name']].attrs.keys() assert (ds_object[result['qc_variable_name']].attrs['warn_max'].dtype == ds_object[result['variable_name']].values.dtype) assert np.isclose(ds_object[result['qc_variable_name']].attrs['warn_max'], limit_value) + result = ds_object.qcfilter.add_greater_equal_test(var_name, limit_value) + assert 'fail_max' in ds_object[result['qc_variable_name']].attrs.keys() + # equal to test limit_value = 7.6705 - result = ds_object.qcfilter.add_equal_to_test(var_name, limit_value) + result = ds_object.qcfilter.add_equal_to_test(var_name, limit_value, prepend_text='arm') data = ds_object.qcfilter.get_masked_data(var_name, rm_tests=result['test_number']) + assert 'arm' in result['test_meaning'] assert np.ma.count_masked(data) == 2 assert 'fail_equal_to' in ds_object[result['qc_variable_name']].attrs.keys() assert (ds_object[result['qc_variable_name']].attrs['fail_equal_to'].dtype == ds_object[result['variable_name']].values.dtype) assert np.isclose(ds_object[result['qc_variable_name']].attrs['fail_equal_to'], limit_value) + result = ds_object.qcfilter.add_equal_to_test(var_name, limit_value, + test_assessment='Indeterminate') + assert 'warn_equal_to' in ds_object[result['qc_variable_name']].attrs.keys() + # not equal to test limit_value = 7.6705 result = ds_object.qcfilter.add_not_equal_to_test(var_name, limit_value, - test_assessment='Indeterminate') + test_assessment='Indeterminate', + prepend_text='arm') data = ds_object.qcfilter.get_masked_data(var_name, rm_tests=result['test_number']) + assert 'arm' in result['test_meaning'] assert np.ma.count_masked(data) == 4318 assert 'warn_not_equal_to' in ds_object[result['qc_variable_name']].attrs.keys() assert (ds_object[result['qc_variable_name']].attrs['warn_not_equal_to'].dtype == ds_object[result['variable_name']].values.dtype) assert np.isclose(ds_object[result['qc_variable_name']].attrs['warn_not_equal_to'], limit_value) + result = ds_object.qcfilter.add_not_equal_to_test(var_name, limit_value) + assert 'fail_not_equal_to' in ds_object[result['qc_variable_name']].attrs.keys() + # outside range test limit_value1 = 6.8 limit_value2 = 12.7 - result = ds_object.qcfilter.add_outside_test(var_name, limit_value1, limit_value2) + result = ds_object.qcfilter.add_outside_test(var_name, limit_value1, limit_value2, + prepend_text='arm') data = ds_object.qcfilter.get_masked_data(var_name, rm_tests=result['test_number']) + assert 'arm' in result['test_meaning'] assert np.ma.count_masked(data) == 115 assert 'fail_lower_range' in ds_object[result['qc_variable_name']].attrs.keys() assert (ds_object[result['qc_variable_name']].attrs['fail_lower_range'].dtype == @@ -277,12 +338,19 @@ def test_qctests(): ds_object[result['variable_name']].values.dtype) assert np.isclose(ds_object[result['qc_variable_name']].attrs['fail_upper_range'], limit_value2) + result = ds_object.qcfilter.add_outside_test(var_name, limit_value1, limit_value2, + test_assessment='Indeterminate') + assert 'warn_lower_range' in ds_object[result['qc_variable_name']].attrs.keys() + assert 'warn_upper_range' in ds_object[result['qc_variable_name']].attrs.keys() + # inside range test limit_value1 = 7 limit_value2 = 8 - result = ds_object.qcfilter.add_inside_test(var_name, limit_value1, limit_value2) + result = ds_object.qcfilter.add_inside_test(var_name, limit_value1, limit_value2, + prepend_text='arm') data = ds_object.qcfilter.get_masked_data(var_name, rm_tests=result['test_number']) + assert 'arm' in result['test_meaning'] assert np.ma.count_masked(data) == 479 assert 'fail_lower_range_inner' in ds_object[result['qc_variable_name']].attrs.keys() assert (ds_object[result['qc_variable_name']].attrs['fail_lower_range_inner'].dtype == @@ -295,11 +363,17 @@ def test_qctests(): assert np.isclose(ds_object[result['qc_variable_name']].attrs['fail_upper_range_inner'], limit_value2) + result = ds_object.qcfilter.add_inside_test(var_name, limit_value1, limit_value2, + test_assessment='Indeterminate') + assert 'warn_lower_range_inner' in ds_object[result['qc_variable_name']].attrs.keys() + assert 'warn_upper_range_inner' in ds_object[result['qc_variable_name']].attrs.keys() + # delta test test_limit = 0.05 - result = ds_object.qcfilter.add_delta_test(var_name, test_limit) + result = ds_object.qcfilter.add_delta_test(var_name, test_limit, prepend_text='arm') data = ds_object.qcfilter.get_masked_data(var_name, rm_tests=result['test_number']) + assert 'arm' in result['test_meaning'] assert np.ma.count_masked(data) == 175 assert 'warn_delta' in ds_object[result['qc_variable_name']].attrs.keys() assert (ds_object[result['qc_variable_name']].attrs['warn_delta'].dtype == @@ -308,20 +382,46 @@ def test_qctests(): data = ds_object.qcfilter.get_masked_data(var_name, rm_assessments=['Suspect', 'Bad']) - assert np.ma.count_masked(data) == 1235 + assert np.ma.count_masked(data) == 4320 + + result = ds_object.qcfilter.add_delta_test(var_name, test_limit, test_assessment='Bad') + assert 'fail_delta' in ds_object[result['qc_variable_name']].attrs.keys() comp_object = read_netcdf(EXAMPLE_IRT25m20s) result = ds_object.qcfilter.add_difference_test( var_name, {comp_object.attrs['datastream']: comp_object}, - var_name, diff_limit=1) + var_name, diff_limit=1, prepend_text='arm') data = ds_object.qcfilter.get_masked_data(var_name, rm_tests=result['test_number']) + assert 'arm' in result['test_meaning'] assert not (data.mask).all() comp_object.close() ds_object.close() +def test_qctests_dos(): + ds_object = read_netcdf(EXAMPLE_IRT25m20s) + var_name = 'inst_up_long_dome_resist' + + # persistence test + data = ds_object[var_name].values + data[1000:2500] = data[1000] + ds_object[var_name].values = data + ds_object.qcfilter.add_persistence_test(var_name) + qc_var_name = ds_object.qcfilter.check_for_ancillary_qc( + var_name, add_if_missing=False, cleanup=False, flag_type=False) + test_meaning = ('Data failing persistence test. Standard Deviation over a ' + 'window of 10 values less than 0.0001.') + assert ds_object[qc_var_name].attrs['flag_meanings'][-1] == test_meaning + assert np.sum(ds_object[qc_var_name].values) == 1500 + + ds_object.qcfilter.add_persistence_test(var_name, window=10000, prepend_text='DQO') + test_meaning = ('DQO: Data failing persistence test. Standard Deviation over a window of ' + '4320 values less than 0.0001.') + assert ds_object[qc_var_name].attrs['flag_meanings'][-1] == test_meaning + + def test_datafilter(): ds = read_netcdf(EXAMPLE_MET1) ds.clean.cleanup() @@ -390,3 +490,110 @@ def test_qc_flag_description(): assert len(ds[qc_var_name].attrs['flag_masks']) == 9 unique_flag_assessments = list(set(['Acceptable', 'Indeterminate', 'Bad'])) assert list(set(ds[qc_var_name].attrs['flag_assessments'])) == unique_flag_assessments + + +def test_clean(): + # Read test data + ceil_ds = read_netcdf([EXAMPLE_CEIL1]) + # Cleanup QC data + ceil_ds.clean.cleanup(clean_arm_state_vars=['detection_status']) + + # Check that global attribures are removed + global_attributes = ['qc_bit_comment', + 'qc_bit_1_description', + 'qc_bit_1_assessment', + 'qc_bit_2_description', + 'qc_bit_2_assessment' + 'qc_bit_3_description', + 'qc_bit_3_assessment' + ] + + for glb_att in global_attributes: + assert glb_att not in ceil_ds.attrs.keys() + + # Check that CF attributes are set including new flag_assessments + var_name = 'qc_first_cbh' + for attr_name in ['flag_masks', 'flag_meanings', 'flag_assessments']: + assert attr_name in ceil_ds[var_name].attrs.keys() + assert isinstance(ceil_ds[var_name].attrs[attr_name], list) + + # Check that the flag_mask values are set correctly + assert ceil_ds['qc_first_cbh'].attrs['flag_masks'] == [1, 2, 4] + + # Check that the flag_meanings values are set correctly + assert (ceil_ds['qc_first_cbh'].attrs['flag_meanings'] == + ['Value is equal to missing_value.', + 'Value is less than the fail_min.', + 'Value is greater than the fail_max.']) + + # Check the value of flag_assessments is as expected + assert ceil_ds['qc_first_cbh'].attrs['flag_assessments'] == ['Bad', 'Bad', 'Bad'] + + # Check that ancillary varibles is being added + assert 'qc_first_cbh' in ceil_ds['first_cbh'].attrs['ancillary_variables'].split() + + # Check that state field is updated to CF + assert 'flag_values' in ceil_ds['detection_status'].attrs.keys() + assert isinstance(ceil_ds['detection_status'].attrs['flag_values'], list) + assert ceil_ds['detection_status'].attrs['flag_values'] == [0, 1, 2, 3, 4, 5] + + assert 'flag_meanings' in ceil_ds['detection_status'].attrs.keys() + assert isinstance(ceil_ds['detection_status'].attrs['flag_meanings'], list) + assert (ceil_ds['detection_status'].attrs['flag_meanings'] == + ['No significant backscatter', + 'One cloud base detected', + 'Two cloud bases detected', + 'Three cloud bases detected', + 'Full obscuration determined but no cloud base detected', + 'Some obscuration detected but determined to be transparent']) + + assert 'flag_0_description' not in ceil_ds['detection_status'].attrs.keys() + assert ('detection_status' in + ceil_ds['first_cbh'].attrs['ancillary_variables'].split()) + + ceil_ds.close() + + +def test_compare_time_series_trends(): + + drop_vars = ['base_time', 'time_offset', 'atmos_pressure', 'qc_atmos_pressure', + 'temp_std', 'rh_mean', 'qc_rh_mean', 'rh_std', 'vapor_pressure_mean', + 'qc_vapor_pressure_mean', 'vapor_pressure_std', 'wspd_arith_mean', + 'qc_wspd_arith_mean', 'wspd_vec_mean', 'qc_wspd_vec_mean', 'wdir_vec_mean', + 'qc_wdir_vec_mean', 'wdir_vec_std', 'tbrg_precip_total', 'qc_tbrg_precip_total', + 'tbrg_precip_total_corr', 'qc_tbrg_precip_total_corr', 'org_precip_rate_mean', + 'qc_org_precip_rate_mean', 'pwd_err_code', 'pwd_mean_vis_1min', 'qc_pwd_mean_vis_1min', + 'pwd_mean_vis_10min', 'qc_pwd_mean_vis_10min', 'pwd_pw_code_inst', + 'qc_pwd_pw_code_inst', 'pwd_pw_code_15min', 'qc_pwd_pw_code_15min', + 'pwd_pw_code_1hr', 'qc_pwd_pw_code_1hr', 'pwd_precip_rate_mean_1min', + 'qc_pwd_precip_rate_mean_1min', 'pwd_cumul_rain', 'qc_pwd_cumul_rain', + 'pwd_cumul_snow', 'qc_pwd_cumul_snow', 'logger_volt', 'qc_logger_volt', + 'logger_temp', 'qc_logger_temp', 'lat', 'lon', 'alt'] + ds = read_netcdf(EXAMPLE_MET1, drop_variables=drop_vars) + ds.clean.cleanup() + ds2 = copy.deepcopy(ds) + + var_name = 'temp_mean' + qc_var_name = ds.qcfilter.check_for_ancillary_qc(var_name, add_if_missing=False, + cleanup=False, flag_type=False) + ds.qcfilter.compare_time_series_trends(var_name=var_name, time_shift=60, + comp_var_name=var_name, comp_dataset=ds2, + time_qc_threshold=60 * 10) + + test_description = ('Time shift detected with Minimum Difference test. Comparison of ' + 'temp_mean with temp_mean off by 0 seconds exceeding absolute ' + 'threshold of 600 seconds.') + assert ds[qc_var_name].attrs['flag_meanings'][-1] == test_description + + time = ds2['time'].values + np.timedelta64(1, 'h') + time_attrs = ds2['time'].attrs + ds2 = ds2.assign_coords({'time': time}) + ds2['time'].attrs = time_attrs + + ds.qcfilter.compare_time_series_trends(var_name=var_name, comp_dataset=ds2, time_step=60, + time_match_threshhold=50) + + test_description = ('Time shift detected with Minimum Difference test. Comparison of ' + 'temp_mean with temp_mean off by 3600 seconds exceeding absolute ' + 'threshold of 900 seconds.') + assert ds[qc_var_name].attrs['flag_meanings'][-1] == test_description diff --git a/act/tests/test_retrievals.py b/act/tests/test_retrievals.py index 78f9a6ac04..a7b0748083 100644 --- a/act/tests/test_retrievals.py +++ b/act/tests/test_retrievals.py @@ -20,10 +20,10 @@ def test_get_stability_indices(): np.testing.assert_allclose( sonde_ds["parcel_temperature"].values[0:5], [269.85000005, 269.74530704, 269.67805708, - 269.62251119, 269.57241322]) + 269.62251119, 269.57241322], rtol=1e-5) assert sonde_ds["parcel_temperature"].attrs["units"] == "kelvin" np.testing.assert_almost_equal( - sonde_ds["surface_based_cape"], 1.628, decimal=3) + sonde_ds["surface_based_cape"], 1.62, decimal=2) assert sonde_ds["surface_based_cape"].attrs["units"] == "J/kg" assert sonde_ds[ "surface_based_cape"].attrs["long_name"] == "Surface-based CAPE" @@ -44,7 +44,7 @@ def test_get_stability_indices(): "most_unstable_cin"].attrs["long_name"] == "Most unstable CIN" np.testing.assert_almost_equal( - sonde_ds["lifted_index"], 28.4639, decimal=3) + sonde_ds["lifted_index"], 28.4, decimal=1) assert sonde_ds["lifted_index"].attrs["units"] == "kelvin" assert sonde_ds["lifted_index"].attrs["long_name"] == "Lifted index" np.testing.assert_equal( @@ -56,7 +56,7 @@ def test_get_stability_indices(): "long_name"] == "Level of free convection" np.testing.assert_almost_equal( sonde_ds["lifted_condensation_level_temperature"], - -8.079, decimal=3) + -8.07, decimal=2) assert sonde_ds[ "lifted_condensation_level_temperature"].attrs[ "units"] == "degree_Celsius" @@ -64,7 +64,7 @@ def test_get_stability_indices(): "lifted_condensation_level_temperature"].attrs[ "long_name"] == "Lifted condensation level temperature" np.testing.assert_almost_equal( - sonde_ds["lifted_condensation_level_pressure"], 927.143, decimal=3) + sonde_ds["lifted_condensation_level_pressure"], 927.1, decimal=1) assert sonde_ds[ "lifted_condensation_level_pressure"].attrs[ "units"] == "hectopascal" @@ -159,16 +159,13 @@ def test_calculate_sirs_variable(): act.tests.sample_files.EXAMPLE_MET1) obj = act.retrievals.radiation.calculate_dsh_from_dsdh_sdn(sirs_object) - assert 61159 <= np.ceil( - np.nansum(obj['derived_down_short_hemisp'].values)) <= 61160 + assert np.isclose(np.nansum(obj['derived_down_short_hemisp'].values), 60350, atol=1) obj = act.retrievals.radiation.calculate_irradiance_stats( obj, variable='derived_down_short_hemisp', variable2='down_short_hemisp', threshold=60) - assert 1336 <= np.ceil( - np.nansum(obj['diff_derived_down_short_hemisp'].values)) <= 1337 - assert np.ceil( - np.nansum(obj['ratio_derived_down_short_hemisp'].values)) == 401 + assert np.isclose(np.nansum(obj['diff_derived_down_short_hemisp'].values), 527, atol=1) + assert np.isclose(np.nansum(obj['ratio_derived_down_short_hemisp'].values), 392, atol=1) obj = act.retrievals.radiation.calculate_net_radiation(obj, smooth=30) assert np.ceil(np.nansum(obj['net_radiation'].values)) == 21915 diff --git a/act/tests/test_utils.py b/act/tests/test_utils.py index 32f5c59fb9..3c24f4757c 100644 --- a/act/tests/test_utils.py +++ b/act/tests/test_utils.py @@ -1,27 +1,20 @@ """ Unit tests for ACT utils module. """ from datetime import datetime +import pytz import act import numpy as np import pandas as pd import pytest import xarray as xr - - -def test_astral_optional_dependency(): - try: - from astral import foo - except ImportError: - act.utils.geo_utils.ASTRAL = False - assert act.utils.geo_utils.ASTRAL is False - act.utils.geo_utils.ASTRAL = True +import tempfile +from pathlib import Path def test_dates_between(): start_date = '20190101' end_date = '20190110' - date_list = act.utils.dates_between(start_date, end_date) answer = [datetime(2019, 1, 1), datetime(2019, 1, 2), @@ -66,13 +59,22 @@ def test_add_in_nan(): assert(time_filled[8].values == time_answer[8]) assert(time_filled[5].values == time_answer[5]) + time_filled, data_filled = act.utils.add_in_nan( + time_list[0], data[0]) + assert time_filled.values == time_list[0] + assert data_filled.values == data[0] + def test_get_missing_value(): obj = act.io.armfiles.read_netcdf(act.tests.sample_files.EXAMPLE_EBBR1) missing = act.utils.data_utils.get_missing_value( - obj, 'lv_e', use_FillValue=True) + obj, 'latent_heat_flux', use_FillValue=True, add_if_missing_in_obj=True) assert missing == -9999 + obj['latent_heat_flux'].attrs['missing_value'] = -9998 + missing = act.utils.data_utils.get_missing_value(obj, 'latent_heat_flux') + assert missing == -9998 + def test_convert_units(): obj = act.io.armfiles.read_netcdf(act.tests.sample_files.EXAMPLE_EBBR1) @@ -84,6 +86,47 @@ def test_convert_units(): data = act.utils.data_utils.convert_units(r_data, 'K', 'C') assert np.ceil(data[0]) == 12 + try: + obj.utils.change_units() + except ValueError as error: + assert str(error) == "Need to provide 'desired_unit' keyword for .change_units() method" + + desired_unit = 'degF' + skip_vars = [ii for ii in obj.data_vars if ii.startswith('qc_')] + obj.utils.change_units(variables=None, desired_unit=desired_unit, + skip_variables=skip_vars, skip_standard=True) + units = [] + for var_name in obj.data_vars: + try: + units.append(obj[var_name].attrs['units']) + except KeyError: + pass + indices = [i for i, x in enumerate(units) if x == desired_unit] + assert indices == [0, 2, 4, 6, 8, 32, 34, 36, 38, 40] + + var_name = 'home_signal_15' + desired_unit = 'V' + obj.utils.change_units(var_name, desired_unit, skip_variables='lat') + assert obj[var_name].attrs['units'] == desired_unit + + var_names = ['home_signal_15', 'home_signal_30'] + obj.utils.change_units(var_names, desired_unit) + for var_name in var_names: + assert obj[var_name].attrs['units'] == desired_unit + + obj.close() + del obj + + obj = act.io.armfiles.read_netcdf(act.tests.sample_files.EXAMPLE_CEIL1) + var_name = 'range' + desired_unit = 'km' + obj = obj.utils.change_units(var_name, desired_unit) + assert obj[var_name].attrs['units'] == desired_unit + assert np.isclose(np.sum(obj[var_name].values), 952.56, atol=0.01) + + obj.close() + del obj + def test_ts_weighted_average(): obj = act.io.armfiles.read_netcdf(act.tests.sample_files.EXAMPLE_MET_WILDCARD) @@ -100,11 +143,18 @@ def test_accum_precip(): act.tests.sample_files.EXAMPLE_MET_WILDCARD) obj = act.utils.accumulate_precip(obj, 'tbrg_precip_total') - dmax = round(np.nanmax(obj['tbrg_precip_total_accumulated'])) + assert dmax == 13.0 + obj = act.utils.accumulate_precip(obj, 'tbrg_precip_total', time_delta=60) + dmax = round(np.nanmax(obj['tbrg_precip_total_accumulated'])) assert dmax == 13.0 + obj['tbrg_precip_total'].attrs['units'] = 'mm/hr' + obj = act.utils.accumulate_precip(obj, 'tbrg_precip_total') + dmax = np.round(np.nanmax(obj['tbrg_precip_total_accumulated']), 2) + assert dmax == 0.22 + def test_calc_cog_sog(): obj = act.io.armfiles.read_netcdf( @@ -154,6 +204,14 @@ def test_calculate_dqr_times(): assert ebbr2_result == [('2019-11-30 00:00:00', '2019-11-30 11:00:00')] assert brs_result == [('2019-07-05 01:57:00', '2019-07-05 11:07:00')] assert ebbr3_result is None + with tempfile.TemporaryDirectory() as tmpdirname: + write_file = Path(tmpdirname) + brs_result = act.utils.calculate_dqr_times(brs_ds, variable='down_short_hemisp_min', + qc_bit=2, threshold=30, txt_path=str(write_file)) + + brs_result = act.utils.calculate_dqr_times(brs_ds, variable='down_short_hemisp_min', + qc_bit=2, threshold=30, return_missing=False) + assert len(brs_result[0]) == 2 ebbr1_ds.close() ebbr2_ds.close() @@ -247,16 +305,16 @@ def test_add_solar_variable(): assert 'sun_variable' in list(new_obj.keys()) assert new_obj['sun_variable'].values[10] == 1 - assert np.sum(new_obj['sun_variable'].values) >= 519 + assert np.sum(new_obj['sun_variable'].values) >= 598 new_obj = act.utils.geo_utils.add_solar_variable(obj, dawn_dusk=True) assert 'sun_variable' in list(new_obj.keys()) assert new_obj['sun_variable'].values[10] == 1 - assert np.sum(new_obj['sun_variable'].values) >= 519 + assert np.sum(new_obj['sun_variable'].values) >= 1234 obj = act.io.armfiles.read_netcdf(act.tests.EXAMPLE_MET1) new_obj = act.utils.geo_utils.add_solar_variable(obj, dawn_dusk=True) - assert np.sum(new_obj['sun_variable'].values) >= 582 + assert np.sum(new_obj['sun_variable'].values) >= 1046 obj = act.io.armfiles.read_netcdf(act.tests.EXAMPLE_IRTSST) obj = obj.fillna(0) @@ -297,3 +355,99 @@ def test_planck_converter(): result = act.utils.radiance_utils.planck_converter(wnum=wnum, radiance=radiance) assert np.ceil(result) == temp np.testing.assert_raises(ValueError, act.utils.radiance_utils.planck_converter) + + +def test_solar_azimuth_elevation(): + + obj = act.io.armfiles.read_netcdf(act.tests.EXAMPLE_NAV) + + elevation, azimuth, distance = act.utils.geo_utils.get_solar_azimuth_elevation( + latitude=obj['lat'].values[0], longitude=obj['lon'].values[0], time=obj['time'].values, + library='skyfield', temperature_C='standard', pressure_mbar='standard') + assert np.isclose(np.nanmean(elevation), 26.408, atol=0.001) + assert np.isclose(np.nanmean(azimuth), 179.732, atol=0.001) + assert np.isclose(np.nanmean(distance), 0.985, atol=0.001) + + +def test_get_sunrise_sunset_noon(): + + obj = act.io.armfiles.read_netcdf(act.tests.EXAMPLE_NAV) + + sunrise, sunset, noon = act.utils.geo_utils.get_sunrise_sunset_noon( + latitude=obj['lat'].values[0], longitude=obj['lon'].values[0], + date=obj['time'].values[0], library='skyfield') + assert sunrise[0].replace(microsecond=0) == datetime(2018, 1, 31, 22, 36, 32) + assert sunset[0].replace(microsecond=0) == datetime(2018, 2, 1, 17, 24, 4) + assert noon[0].replace(microsecond=0) == datetime(2018, 2, 1, 8, 2, 10) + + sunrise, sunset, noon = act.utils.geo_utils.get_sunrise_sunset_noon( + latitude=obj['lat'].values[0], longitude=obj['lon'].values[0], + date=obj['time'].values[0], library='skyfield', timezone=True) + assert sunrise[0].replace(microsecond=0) == datetime(2018, 1, 31, 22, 36, 32, tzinfo=pytz.UTC) + assert sunset[0].replace(microsecond=0) == datetime(2018, 2, 1, 17, 24, 4, tzinfo=pytz.UTC) + assert noon[0].replace(microsecond=0) == datetime(2018, 2, 1, 8, 2, 10, tzinfo=pytz.UTC) + + sunrise, sunset, noon = act.utils.geo_utils.get_sunrise_sunset_noon( + latitude=obj['lat'].values[0], longitude=obj['lon'].values[0], + date='20180201', library='skyfield') + assert sunrise[0].replace(microsecond=0) == datetime(2018, 1, 31, 22, 36, 32) + assert sunset[0].replace(microsecond=0) == datetime(2018, 2, 1, 17, 24, 4) + assert noon[0].replace(microsecond=0) == datetime(2018, 2, 1, 8, 2, 10) + + sunrise, sunset, noon = act.utils.geo_utils.get_sunrise_sunset_noon( + latitude=obj['lat'].values[0], longitude=obj['lon'].values[0], + date=['20180201'], library='skyfield') + assert sunrise[0].replace(microsecond=0) == datetime(2018, 1, 31, 22, 36, 32) + assert sunset[0].replace(microsecond=0) == datetime(2018, 2, 1, 17, 24, 4) + assert noon[0].replace(microsecond=0) == datetime(2018, 2, 1, 8, 2, 10) + + sunrise, sunset, noon = act.utils.geo_utils.get_sunrise_sunset_noon( + latitude=obj['lat'].values[0], longitude=obj['lon'].values[0], + date=datetime(2018, 2, 1), library='skyfield') + assert sunrise[0].replace(microsecond=0) == datetime(2018, 1, 31, 22, 36, 32) + assert sunset[0].replace(microsecond=0) == datetime(2018, 2, 1, 17, 24, 4) + assert noon[0].replace(microsecond=0) == datetime(2018, 2, 1, 8, 2, 10) + + sunrise, sunset, noon = act.utils.geo_utils.get_sunrise_sunset_noon( + latitude=obj['lat'].values[0], longitude=obj['lon'].values[0], + date=datetime(2018, 2, 1, tzinfo=pytz.UTC), library='skyfield') + assert sunrise[0].replace(microsecond=0) == datetime(2018, 1, 31, 22, 36, 32) + assert sunset[0].replace(microsecond=0) == datetime(2018, 2, 1, 17, 24, 4) + assert noon[0].replace(microsecond=0) == datetime(2018, 2, 1, 8, 2, 10) + + sunrise, sunset, noon = act.utils.geo_utils.get_sunrise_sunset_noon( + latitude=obj['lat'].values[0], longitude=obj['lon'].values[0], + date=[datetime(2018, 2, 1)], library='skyfield') + assert sunrise[0].replace(microsecond=0) == datetime(2018, 1, 31, 22, 36, 32) + assert sunset[0].replace(microsecond=0) == datetime(2018, 2, 1, 17, 24, 4) + assert noon[0].replace(microsecond=0) == datetime(2018, 2, 1, 8, 2, 10) + + sunrise, sunset, noon = act.utils.geo_utils.get_sunrise_sunset_noon( + latitude=85.0, longitude=-140., date=[datetime(2018, 6, 1)], library='skyfield') + assert sunrise[0].replace(microsecond=0) == datetime(2018, 3, 30, 10, 48, 48) + assert sunset[0].replace(microsecond=0) == datetime(2018, 9, 12, 8, 50, 14) + assert noon[0].replace(microsecond=0) == datetime(2018, 6, 1, 21, 17, 52) + + +def test_is_sun_visible(): + obj = act.io.armfiles.read_netcdf(act.tests.sample_files.EXAMPLE_EBBR1) + result = act.utils.geo_utils.is_sun_visible( + latitude=obj['lat'].values, longitude=obj['lon'].values, + date_time=obj['time'].values) + assert len(result) == 48 + assert sum(result) == 20 + + result = act.utils.geo_utils.is_sun_visible( + latitude=obj['lat'].values, longitude=obj['lon'].values, + date_time=obj['time'].values[0]) + assert result == [False] + + result = act.utils.geo_utils.is_sun_visible( + latitude=obj['lat'].values, longitude=obj['lon'].values, + date_time=[datetime(2019, 11, 25, 13, 30, 00)]) + assert result == [True] + + result = act.utils.geo_utils.is_sun_visible( + latitude=obj['lat'].values, longitude=obj['lon'].values, + date_time=datetime(2019, 11, 25, 13, 30, 00)) + assert result == [True] diff --git a/act/utils/__init__.py b/act/utils/__init__.py index a33b9dd42b..02a0f20e8e 100644 --- a/act/utils/__init__.py +++ b/act/utils/__init__.py @@ -49,3 +49,4 @@ from .geo_utils import add_solar_variable from .inst_utils import decode_present_weather from .radiance_utils import planck_converter +from .data_utils import ChangeUnits diff --git a/act/utils/conf/de421.bsp b/act/utils/conf/de421.bsp new file mode 100644 index 0000000000..17eece18e5 Binary files /dev/null and b/act/utils/conf/de421.bsp differ diff --git a/act/utils/data_utils.py b/act/utils/data_utils.py index a8284343bc..1eb457f48e 100644 --- a/act/utils/data_utils.py +++ b/act/utils/data_utils.py @@ -20,6 +20,79 @@ PYART_AVAILABLE = False +@xr.register_dataset_accessor('utils') +class ChangeUnits(object): + """ + Class for updating units in the object. Data values and units attribute + are updated in place. Coordinate variables can not be updated in place. Must + use new returned dataset when updating coordinage varibles. + """ + def __init__(self, xarray_obj): + self._obj = xarray_obj + + def change_units(self, variables=None, desired_unit=None, + skip_variables=None, skip_standard=True): + """ + Parameters + ---------- + variables : None, str or list of str + Variable names to attempt to change units + desired_unit : str + Desired udunits unit string + skip_variables : None, str or list of str + Varible names to skip. Works well when not providing a variables keyword + skip_standard : boolean + Flag indicating the QC variables that will not need changing are skipped. + Makes the processing faster when processing all variables in dataset. + + Returns + ------- + dataset : xarray.dataset + A new dataset if the coordinate variables are updated. Required to use + returned dataset if coordinage variabels are updated, otherwise the + dataset is updated in place. + """ + + if variables is not None and isinstance(variables, str): + variables = [variables] + + if skip_variables is not None and isinstance(skip_variables, str): + skip_variables = [skip_variables] + + if desired_unit is None: + raise ValueError("Need to provide 'desired_unit' keyword for .change_units() method") + + if variables is None: + variables = list(self._obj.data_vars) + + if skip_variables is not None: + variables = list(set(variables) - set(skip_variables)) + + for var_name in variables: + try: + if self._obj[var_name].attrs['standard_name'] == 'quality_flag': + continue + except KeyError: + pass + + try: + data = convert_units(self._obj[var_name].values, + self._obj[var_name].attrs['units'], desired_unit) + try: + self._obj[var_name].values = data + self._obj[var_name].attrs['units'] = desired_unit + except ValueError: + attrs = self._obj[var_name].attrs + self._obj = self._obj.assign_coords({var_name: data}) + attrs['units'] = desired_unit + self._obj[var_name].attrs = attrs + except (KeyError, pint.errors.DimensionalityError, pint.errors.UndefinedUnitError, + np.core._exceptions.UFuncTypeError): + continue + + return self._obj + + def assign_coordinates(ds, coord_list): """ This procedure will create a new ACT dataset whose coordinates are @@ -274,7 +347,8 @@ def convert_units(data, in_units, out_units): data_type_kind = data.dtype.kind # Do the conversion magic - data = np.asarray((data * ureg(in_units)).to(out_units)) + data = (data * ureg(in_units)).to(out_units) + data = data.magnitude # The data type may be changed by pint. This is a side effect # of pint changing the datatype to float. Check if the converted values @@ -395,10 +469,6 @@ def accumulate_precip(act_obj, variable, time_delta=None): accum = np.nancumsum(data.values) - # Add time as a variable if not already a variable - if 'time' not in act_obj: - act_obj['time'] = xr.DataArray(time, coords=act_obj[variable].coords) - # Add accumulated variable back to ACT object long_name = 'Accumulated precipitation' attrs = {'long_name': long_name, 'units': 'mm'} diff --git a/act/utils/datetime_utils.py b/act/utils/datetime_utils.py index c19c05db44..7c56c8fc93 100644 --- a/act/utils/datetime_utils.py +++ b/act/utils/datetime_utils.py @@ -39,7 +39,7 @@ def dates_between(sdate, edate): return all_dates -def numpy_to_arm_date(_date): +def numpy_to_arm_date(_date, returnTime=False): """ Given a numpy datetime64, return an ARM standard date (yyyymmdd). @@ -47,6 +47,8 @@ def numpy_to_arm_date(_date): ---------- date : numpy.datetime64 Numpy datetime64 date. + returnTime : boolean + If set to true, returns time instead of date Returns ------- @@ -55,7 +57,10 @@ def numpy_to_arm_date(_date): """ date = pd.to_datetime(str(_date)) - date = date.strftime('%Y%m%d') + if returnTime is False: + date = date.strftime('%Y%m%d') + else: + date = date.strftime('%H%M%S') return date diff --git a/act/utils/geo_utils.py b/act/utils/geo_utils.py index 7067a3b92f..33ebf0bf77 100644 --- a/act/utils/geo_utils.py +++ b/act/utils/geo_utils.py @@ -9,16 +9,20 @@ import numpy as np import pandas as pd -import astral -try: - from astral.sun import sunrise, sunset, dusk, dawn - from astral import Observer - ASTRAL = True -except ImportError: - ASTRAL = False +from datetime import datetime, timezone, timedelta +from skyfield.api import wgs84, N, W, load_file, load +from skyfield import almanac +import re +import dateutil.parser +import pytz +from pathlib import Path +from act.utils.datetime_utils import datetime64_to_datetime +from act.utils.data_utils import convert_units +skyfield_bsp_file = str(Path(Path(__file__).parent, "conf", "de421.bsp")) -def destination_azimuth_distance(lat, lon, az, dist): + +def destination_azimuth_distance(lat, lon, az, dist, dist_units='m'): """ This procedure will calculate a destination lat/lon from an initial lat/lon and azimuth and distance. @@ -32,24 +36,26 @@ def destination_azimuth_distance(lat, lon, az, dist): az : float Azimuth in degrees. dist : float - Distance in meters. + Distance + dist_units : str + Units for dist Returns ------- lat2 : float - Latitude of new point. + Latitude of new point in degrees lon2 : float - Longitude of new point. + Longitude of new point in degrees """ - # Volumetric Mean Radius of Earth - R = 6371. + # Volumetric Mean Radius of Earth in km + R = 6378. # Convert az to radian brng = np.radians(az) - # Assuming meters as input - d = dist / 1000. + # Convert distance to km + d = convert_units(dist, dist_units, 'km') # Convert lat/lon to radians lat = np.radians(lat) @@ -66,21 +72,18 @@ def destination_azimuth_distance(lat, lon, az, dist): def add_solar_variable(obj, latitude=None, longitude=None, solar_angle=0., dawn_dusk=False): """ - Calculate solar times depending on location on earth. - - Astral 2.2 is recommended for best performance and for the dawn/dusk feature as it - seems like the dawn calculations are wrong with earlier versions. + Add variable to the object to denote night (0) or sun (1). If dawk_dusk is True + will also return dawn (2) and dusk (3). If at a high latitude and there's sun, will + label twilight as dawn; if dark{2}, will label twilight as dusk(3). Parameters ---------- - obj : act object + obj : xarray dataset ACT object latitude : str - Latitude variable, default will look for matching variables - in object + Latitude variable name, default will look for matching variables in object longitude : str - Longitude variable, default will look for matching variables - in object + Longitude variable name, default will look for matching variables in object solar_angle : float Number of degress to use for dawn/dusk calculations dawn_dusk : boolean @@ -88,8 +91,8 @@ def add_solar_variable(obj, latitude=None, longitude=None, solar_angle=0., dawn_ Returns ------- - obj : act object - ACT object + obj : xarray dataset + Xarray object """ variables = list(obj.keys()) @@ -113,115 +116,331 @@ def add_solar_variable(obj, latitude=None, longitude=None, solar_angle=0., dawn_ lat = obj[latitude[0]].values lon = obj[longitude[0]].values - # Set the the number of degrees the sun must be below the horizon - # for the dawn/dusk calculation. Need to do this so when the calculation - # sends an error it is not going to be an inacurate switch to setting - # the full day. - if ASTRAL: - astral.solar_depression = solar_angle - else: - a = astral.Astral() - a.solar_depression = 0. - - # If only one lat/lon value then set up the observer location - # for Astral. If more than one, it will get set up in the loop - if lat.size == 1 and ASTRAL: - loc = Observer(latitude=lat, longitude=lon) - # Loop through each time to ensure that the sunrise/set calcuations # are correct for each time and lat/lon if multiple - results = [] - time = obj['time'].values - for i in range(len(time)): - # Set up an observer if multiple lat/lon - if lat.size > 1: - if ASTRAL: - loc = Observer(latitude=lat[i], longitude=lon[i]) - else: - s = a.sun_utc(pd.to_datetime(time[i]), lat[i], lon[i]) - elif ASTRAL is False: - s = a.sun_utc(pd.to_datetime(time[i]), float(lat), float(lon)) - - # Get sunrise and sunset - if ASTRAL: - sr = sunrise(loc, pd.to_datetime(time[i])) - ss = sunset(loc, pd.to_datetime(time[i])) + results = is_sun_visible(latitude=lat, longitude=lon, date_time=obj['time'].values, + dawn_dusk=dawn_dusk) + + # Set longname + longname = 'Daylight indicator; 0-Night; 1-Sun' + + if dawn_dusk is False: + results = results * 1 + else: + # If dawn_dusk is True, add 2 more indicators + longname += '; 2-Dawn; 3-Dusk; 4-Twilight' + dark_ind = np.where((results == 0))[0] + twil_ind = np.where((results > 0) & (results < 4))[0] + sun_ind = np.where((results == 4))[0] + + if len(sun_ind) == 0: + results[twil_ind] = 3 + elif len(dark_ind) == 0: + results[twil_ind] = 2 + results[sun_ind] = 1 else: - sr = s['sunrise'] - ss = s['sunset'] - - # Set longname - longname = 'Daylight indicator; 0-Night; 1-Sun' - - # Check to see if dawn/dusk calculations can be performed before preceeding - if dawn_dusk: - try: - if ASTRAL: - dwn = dawn(loc, pd.to_datetime(time[i])) - dsk = dusk(loc, pd.to_datetime(time[i])) - else: - if lat.size > 1: - dsk = a.dusk_utc(pd.to_datetime(time[i]), lat[i], lon[i]) - dwn = a.dawn_utc(pd.to_datetime(time[i]), lat[i], lon[i]) - else: - dsk = a.dusk_utc(pd.to_datetime(time[i]), float(lat), float(lon)) - dwn = a.dawn_utc(pd.to_datetime(time[i]), float(lat), float(lon)) - except ValueError: - print('Dawn/Dusk calculations are not available at this location') - dawn_dusk = False - - if dawn_dusk and ASTRAL: - # If dawn_dusk is True, add 2 more indicators - longname += '; 2-Dawn; 3-Dusk' - # Need to ensure if the sunset if off a day to grab the previous - # days value to catch the early UTC times - if ss.day > sr.day: - if ASTRAL: - ss = sunset(loc, pd.to_datetime(time[i] - np.timedelta64(1, 'D'))) - dsk = dusk(loc, pd.to_datetime(time[i] - np.timedelta64(1, 'D'))) - else: - if lat.size > 1: - dsk = a.dusk_utc(pd.to_datetime(time[i]) - np.timedelta64(1, 'D'), - lat[i], lon[i]) - s = a.sun_utc(pd.to_datetime(time[i]) - np.timedelta64(1, 'D'), - lat[i], lon[i]) - else: - dsk = a.dusk_utc(pd.to_datetime(time[i]) - np.timedelta64(1, 'D'), - float(lat), float(lon)) - s = a.sun_utc(pd.to_datetime(time[i]) - np.timedelta64(1, 'D'), - float(lat), float(lon)) - ss = s['sunset'] - - if dwn <= pd.to_datetime(time[i], utc=True) < sr: - results.append(2) - elif ss <= pd.to_datetime(time[i], utc=True) < dsk: - results.append(3) - elif not(dsk <= pd.to_datetime(time[i], utc=True) < dwn): - results.append(1) - else: - results.append(0) + # Set Dawn between dark and sun + if dark_ind[-1] < sun_ind[0]: + dawn_ind = list(range(dark_ind[-1], sun_ind[0])) else: - if dwn <= pd.to_datetime(time[i], utc=True) < sr: - results.append(2) - elif sr <= pd.to_datetime(time[i], utc=True) < ss: - results.append(1) - elif ss <= pd.to_datetime(time[i], utc=True) < dsk: - results.append(3) - else: - results.append(0) - else: - if ss.day > sr.day: - if ASTRAL: - ss = sunset(loc, pd.to_datetime(time[i] - np.timedelta64(1, 'D'))) - else: - s = a.sun_utc(pd.to_datetime(time[i]) - np.timedelta64(1, 'D'), lat, lon) - ss = s['sunset'] - results.append(int(not(ss < pd.to_datetime(time[i], utc=True) < sr))) + dawn_ind = list(range(dark_ind[-1], len(results))) + list(range(0, sun_ind[0])) + results[dawn_ind] = 2 + + # Set Dusk between sun and dark + if sun_ind[-1] < dark_ind[0]: + dusk_ind = list(range(sun_ind[-1], dark_ind[0])) else: - results.append(int(sr < pd.to_datetime(time[i], utc=True) < ss)) + dusk_ind = list(range(sun_ind[-1], len(results))) + list(range(0, dark_ind[0])) + results[dusk_ind] = 3 + results[sun_ind] = 1 # Add results to object and return obj['sun_variable'] = ('time', np.array(results), {'long_name': longname, 'units': ' '}) return obj + + +def get_solar_azimuth_elevation(latitude=None, longitude=None, time=None, library='skyfield', + temperature_C='standard', pressure_mbar='standard'): + """ + Calculate solar azimuth, elevation and solar distance. + + + Parameters + ---------- + latitude : int, float + Latitude in degrees north positive. Must be a scalar. + longitude : int, float + Longitude in degrees east positive. Must be a scalar. + time : datetime.datetime, numpy.datetime64, list, numpy.array + Time in UTC. May be a scalar or vector. datetime must be timezone aware. + library : str + Library to use for making calculations. Options include ['skyfield'] + temperature_C : string or list of float + If library is 'skyfield' the temperature in degrees C at the surface for + atmospheric compensation of the positon of the sun. Set to None for no + compensation or 'standard' for standard model with a standard temperature. + pressure_mbar : string or list of float + If library is 'skyfield' the pressure in milibars at the surface for + atmospheric compensation of the positon of the sun. Set to None for no + compensation or 'standard' for standard model with a standard pressure. + + Returns + ------- + result : tuple of float + Values returned are a tuple of elevation, azimuth and distance. Elevation and + azimuth are in degrees, with distance in Astronomical Units. + + """ + + result = {'elevation': None, 'azimuth': None, 'distance': None} + + if library == 'skyfield': + planets = load_file(skyfield_bsp_file) + earth, sun = planets['earth'], planets['sun'] + + if isinstance(time, datetime) and time.tzinfo is None: + time = time.replace(tzinfo=pytz.UTC) + + if isinstance(time, (list, tuple)) and time[0].tzinfo is None: + time = [ii.replace(tzinfo=pytz.UTC) for ii in time] + + if type(time).__module__ == np.__name__ and np.issubdtype(time.dtype, np.datetime64): + time = time.astype('datetime64[s]').astype(int) + if time.size > 1: + time = [datetime.fromtimestamp(tm, timezone.utc) for tm in time] + else: + time = [datetime.fromtimestamp(time, timezone.utc)] + + if not isinstance(time, (list, tuple, np.ndarray)): + time = [time] + + ts = load.timescale() + t = ts.from_datetimes(time) + location = earth + wgs84.latlon(latitude * N, longitude * W) + astrometric = location.at(t).observe(sun) + alt, az, distance = astrometric.apparent().altaz(temperature_C=temperature_C, + pressure_mbar=pressure_mbar) + result = (alt.degrees, az.degrees, distance.au) + planets.close() + + return result + + +def get_sunrise_sunset_noon(latitude=None, longitude=None, date=None, library='skyfield', + timezone=False): + """ + Calculate sunrise, sunset and local solar noon times. + + Parameters + ---------- + latitude : int, float + Latitude in degrees north positive. Must be a scalar. + longitude : int, float + Longitude in degrees east positive. Must be a scalar. + date : (datetime.datetime, numpy.datetime64, list of datetime.datetime, + numpy.array of numpy.datetime64, string, list of string) + Date(s) to return sunrise, sunset and noon times spaning the first date to last + date if more than one provided. May be a scalar or vector. If entered as a string must follow + YYYYMMDD format. + library : str + Library to use for making calculations. Options include ['skyfield'] + timezone : boolean + Have timezone with datetime. + + Returns + ------- + result : tuple of three numpy.array + Tuple of three values sunrise, sunset, noon. Values will be a list. + If no values can be calculated will return empty list. If the date is within + polar night will return empty lists. If spans the transition to polar day + will return previous sunrise or next sunset outside of date range provided. + """ + + sunrise, sunset, noon = np.array([]), np.array([]), np.array([]) + + if library == 'skyfield': + ts = load.timescale() + eph = load_file(skyfield_bsp_file) + sf_dates = [] + + # Parse datetime object + if isinstance(date, datetime): + if date.tzinfo is None: + sf_dates = [date.replace(tzinfo=pytz.UTC)] + else: + sf_dates = [date] + + if isinstance(date, (list, tuple)) and isinstance(date[0], datetime): + if date[0].tzinfo is not None: + sf_dates = date + else: + sf_dates = [ii.replace(tzinfo=pytz.UTC) for ii in date] + + # Parse string date + if isinstance(date, str): + sf_dates = [datetime.strptime(date, '%Y%m%d').replace(tzinfo=pytz.UTC)] + + # Parse list of string dates + if isinstance(date, (list, tuple)) and isinstance(date[0], str): + sf_dates = [datetime.strptime(dt, '%Y%m%d').replace(tzinfo=pytz.UTC) for dt in date] + + # Convert datetime64 to datetime + if type(date).__module__ == np.__name__ and np.issubdtype(date.dtype, np.datetime64): + sf_dates = datetime64_to_datetime(date) + sf_dates = [ii.replace(tzinfo=pytz.UTC) for ii in sf_dates] + + # Function for calculating solar noon + # Convert location into skyfield location object + location = wgs84.latlon(latitude, longitude) + # Set up function to indicate calculating locatin of Sun from Earth + f = almanac.meridian_transits(eph, eph['Sun'], location) + # Set up dates to be start of day and end of day so have a range + t0 = sf_dates[0] + t0 = t0.replace(hour=0, minute=0, second=0) + t1 = sf_dates[-1] + t1 = t1.replace(hour=23, minute=59, second=59) + # Convert times from datetime to skyfild times + t0 = ts.from_datetime(t0) + t1 = ts.from_datetime(t1) + # Calculate Meridian Transits. n contains times and x contains 1 and 0's + # indicating when transit time is above or below location. + n, x = almanac.find_discrete(t0, t1, f) + + # Determine if time is during daylight + f = almanac.sunrise_sunset(eph, location) + sun_up = f(n) + + # Filter out times when sun is below location or in polar night + n = n[(x == 1) & sun_up] + noon = n.utc_datetime() + if noon.size == 0: + return sunrise, sunset, noon + + # Calcuate sunrise and sunset times. Calcuate over range 12 less than minimum + # noon time and 12 hours greater than maximum noon time. + t0 = min(noon) - timedelta(hours=12) + t1 = max(noon) + timedelta(hours=12) + t0 = ts.from_datetime(t0) + t1 = ts.from_datetime(t1) + f = almanac.sunrise_sunset(eph, location) + t, y = almanac.find_discrete(t0, t1, f) + times = t.utc_datetime() + sunrise = times[y == 1] + sunset = times[y == 0] + + # Fill in sunrise and sunset if asked to during polar day + if len(noon) > 0 and (y.size == 0 or len(sunrise) != len(sunset)): + days = 200 + t0 = min(noon) - timedelta(days=days) + t1 = max(noon) + timedelta(days=days) + t0 = ts.from_datetime(t0) + t1 = ts.from_datetime(t1) + t, yy = almanac.find_discrete(t0, t1, f) + times = t.utc_datetime() + + # If first time is sunset and/or last time is sunrise filter + # from times + if yy[0] == 0: + yy = yy[1:] + times = times[1:] + if yy[-1] == 1: + yy = yy[:-1] + times = times[:-1] + + # Extract sunrise times + temp_sunrise = times[yy == 1] + # Extract sunset times + temp_sunset = times[yy == 0] + # Look for index closest to first noon time to get the time of last sunrise + # since we are in polar day. + diff = temp_sunrise - min(noon) + sunrise_index = np.max(np.where(diff < timedelta(seconds=1))) + # Look for index closest to last noon time to get the time of first sunset + # since we are in polar day. + diff = max(noon) - temp_sunset + sunset_index = np.min(np.where(diff < timedelta(seconds=1))) + 1 + sunrise = temp_sunrise[sunrise_index: sunset_index] + sunset = temp_sunset[sunrise_index: sunset_index] + + eph.close() + + if timezone is False: + for ii in range(0, sunset.size): + sunrise[ii] = sunrise[ii].replace(tzinfo=None) + sunset[ii] = sunset[ii].replace(tzinfo=None) + for ii in range(0, noon.size): + noon[ii] = noon[ii].replace(tzinfo=None) + + return sunrise, sunset, noon + + +def is_sun_visible(latitude=None, longitude=None, date_time=None, dawn_dusk=False): + """ + Determine if sun is above horizon at for a list of times. + + Parameters + ---------- + latitude : int, float + Latitude in degrees north positive. Must be a scalar. + longitude : int, float + Longitude in degrees east positive. Must be a scalar. + date_time : datetime.datetime, numpy.array.datetime64, list of datetime.datetime + Datetime with timezone, datetime with no timezone in UTC, or numpy.datetime64 + format in UTC. Can be a single datetime object or list of datetime objects. + dawn_dusk : boolean + If set to True, will use skyfields dark_twilight_day function to calculate sun up + Returns a list of int's instead of boolean. + 0 - Dark of Night + 1 - Astronomical Twilight + 2 - Nautical Twilight + 3 - Civil Twilight + 4 - Sun Is Up + + Returns + ------- + result : list + List matching size of date_time containing True/False if sun is above horizon. + """ + + sf_dates = None + + # Check if datetime object is scalar and if has no timezone. + if isinstance(date_time, datetime): + if date_time.tzinfo is None: + sf_dates = [date_time.replace(tzinfo=pytz.UTC)] + else: + sf_dates = [date_time] + + # Check if datetime objects in list have timezone. If not add. + if isinstance(date_time, (list, tuple)) and isinstance(date_time[0], datetime): + if isinstance(date_time[0], datetime) and date_time[0].tzinfo is not None: + sf_dates = date_time + else: + sf_dates = [ii.replace(tzinfo=pytz.UTC) for ii in date_time] + + # Convert datetime64 to datetime with timezone. + if type(date_time).__module__ == np.__name__ and np.issubdtype(date_time.dtype, np.datetime64): + sf_dates = datetime64_to_datetime(date_time) + sf_dates = [ii.replace(tzinfo=pytz.UTC) for ii in sf_dates] + + if sf_dates is None: + raise ValueError('The date_time values entered into is_sun_visible() ' + 'do not match input types.') + + ts = load.timescale() + eph = load_file(skyfield_bsp_file) + + t0 = ts.from_datetimes(sf_dates) + location = wgs84.latlon(latitude, longitude) + if dawn_dusk: + f = almanac.dark_twilight_day(eph, location) + else: + f = almanac.sunrise_sunset(eph, location) + + sun_up = f(t0) + + eph.close() + + return sun_up diff --git a/continuous_integration/environment-3.7.yml b/continuous_integration/environment-3.7.yml index 20872216e2..710779a3a8 100644 --- a/continuous_integration/environment-3.7.yml +++ b/continuous_integration/environment-3.7.yml @@ -8,7 +8,7 @@ dependencies: - pyproj - numpy - scipy - - matplotlib<=3.2.2 + - matplotlib - netcdf4 - pytest - pytest-mpl @@ -18,8 +18,8 @@ dependencies: - flake8 - ipython - pint=0.8.1 + - skyfield - pip - pip: - mpl2nc - arm-pyart - - astral>=2.2 diff --git a/continuous_integration/environment-3.6.yml b/continuous_integration/environment-3.8.yml similarity index 82% rename from continuous_integration/environment-3.6.yml rename to continuous_integration/environment-3.8.yml index 039224c338..aa85455918 100644 --- a/continuous_integration/environment-3.6.yml +++ b/continuous_integration/environment-3.8.yml @@ -3,12 +3,12 @@ channels: - conda-forge - defaults dependencies: - - python=3.6 + - python=3.8 - cartopy - pyproj - numpy - scipy - - matplotlib<=3.2.2 + - matplotlib - netcdf4 - pytest - pytest-mpl @@ -18,8 +18,8 @@ dependencies: - flake8 - ipython - pint=0.8.1 + - skyfield - pip - pip: - mpl2nc - arm-pyart - - astral>=2.2 diff --git a/continuous_integration/install.sh b/continuous_integration/install.sh index a1597c5d4b..4fc61a1768 100644 --- a/continuous_integration/install.sh +++ b/continuous_integration/install.sh @@ -20,6 +20,5 @@ conda info -a ## Create a testenv with the correct Python version conda env create -f continuous_integration/environment-$PYTHON_VERSION.yml source activate testenv -conda install boto3 pip install -e . pip install sphinx sphinx_rtd_theme diff --git a/docs/source/act_plots.png b/docs/source/act_plots.png new file mode 100644 index 0000000000..c0d6dbbf62 Binary files /dev/null and b/docs/source/act_plots.png differ diff --git a/docs/source/index.rst b/docs/source/index.rst index 0120c962ff..f621d7ecaa 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -38,37 +38,35 @@ Atmospheric data Community Toolkit Documentation Atmospheric Community Toolkit (ACT) =================================== +The Atmospheric data Community Toolkit (ACT) is an open source Python toolkit for working with atmospheric time-series datasets of varying dimensions. The toolkit is meant to have functions for every part of the scientific process; discovery, IO, quality control, corrections, retrievals, visualization, and analysis. It is meant to be a community platform for sharing code with the goal of reducing duplication of effort and better connecting the science community with programs such as the `Atmospheric Radiation Measurement (ARM) User Facility `_. Overarching development goals will be updated on a regular basis as part of the `Roadmap `_ . -The Atmospheric Community Toolkit (ACT) is a Python toolkit for working with -atmospheric time-series datasets of varying dimensions. The toolkit is meant -to have functions for every part of the scientific process; discovery, IO, -quality control, corrections, retrievals, visualization, and analysis. Initial -efforts were heavily focused on the static visualization aspect of the process, -but future efforts will look to build up the other areas of interest include -discovery, corrections, retirevals, and interactive plots. +|act| -* Free software: 3-clause BSD license - -|ceil| - -.. |ceil| image:: plot_ceil_example.png +.. |act| image:: act_plots.png Dependencies ============ +* `xarray `_ * `NumPy `_ * `SciPy `_ * `matplotlib `_ -* `xarray `_ -* `astral `_ +* `skyfield `_ * `pandas `_ * `dask `_ * `Pint `_ -* `Cartopy `_ -* `Boto3 `_ * `PyProj `_ +* `Proj `_ +* `Six `_ * `Requests `_ +Optional Dependencies +===================== + +* `MPL2NC `_ Reading binary MPL data. +* `Cartopy `_ Mapping and geoplots +* `MetPy `_ >= V1.0 Skew-T plotting and some stabilities indices calculations + Contributing ============ @@ -93,7 +91,7 @@ Testing After installation, you can launch the test suite from outside the source directory (you will need to have pytest installed):: - $ pytest --pyargs act + $ pytest --mpl --pyargs act In-place installs can be tested using the `pytest` command from within the source directory. diff --git a/environment.yml b/environment.yml index c85d5fd63a..29256d8537 100644 --- a/environment.yml +++ b/environment.yml @@ -14,5 +14,6 @@ dependencies: - pint - requests - pyproj>=2.0.0 - - pip - - astral + - cartopy + - skyfield + diff --git a/examples/change_units.py b/examples/change_units.py new file mode 100644 index 0000000000..6844c9eb35 --- /dev/null +++ b/examples/change_units.py @@ -0,0 +1,61 @@ +""" +===================================== +Example for changing units in dataset +===================================== + +This is an example of how to change +units in the xarray dataset. + +""" + +import act +import numpy as np + + +def print_summary(obj, variables): + for var_name in variables: + print(f"{var_name}: mean={np.nanmean(obj[var_name].values)} " + f"units={obj[var_name].attrs['units']}") + print() + + +variables = ['first_cbh', 'second_cbh', 'alt'] + +# Read in some example data +obj = act.io.armfiles.read_netcdf(act.tests.sample_files.EXAMPLE_CEIL1) + +# Print the variable name, mean of values and units +print('Variables in read data') +print_summary(obj, variables) + +# Change units of one varible from m to km +obj.utils.change_units(variables='first_cbh', desired_unit='km') +print('Variables with one changed to km') +print_summary(obj, variables) + +# Change units of more than one varible from to km +obj.utils.change_units(variables=variables, desired_unit='km') +print('Variables with both changed to km') +print_summary(obj, variables) + +# Can change all data variables in the dataset that are units of length by not providing +# a list of variables. Here we are changing back to orginal meters. +# Also, because it needs to loop over all variables and try to convert, will take +# longer if we keep the QC variables. Faseter if we exclude them. +# The method will return a dataset. In this case the dataset returned is the same +# dataset. +skip_variables = [ii for ii in obj.data_vars if ii.startswith('qc_')] +new_obj = obj.utils.change_units(variables=None, desired_unit='m', skip_variables=skip_variables) +print('Variables changed back to m by looping over all variables in dataset') +print('Orginal dataset is same as retured dataset:', obj is new_obj) +print_summary(new_obj, variables) + +# For coordinate variables need to explicitly give coordinage variable name and use +# the returned dataset. The xarray method used to change values on coordinate +# values requries returning a new updated dataet. +var_name = 'range' +variables.append(var_name) +new_obj = obj.utils.change_units(variables=variables, desired_unit='km') +print('Variables and coordinate variable values changed to km') +print('Orginal dataset is same as retured dataset:', obj is new_obj) +print_summary(new_obj, variables) diff --git a/requirements.txt b/requirements.txt index 693f8bfdc1..3530f5e642 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,15 +1,15 @@ # List required packages in this file, one per line. -cartopy pyproj numpy pandas matplotlib scipy xarray -astral dask distributed pint -boto3 requests six +skyfield +cftime +netcdf4 diff --git a/setup.py b/setup.py index 9c997783d9..e340c380e3 100644 --- a/setup.py +++ b/setup.py @@ -43,7 +43,7 @@ author="Adam Theisen", author_email='atheisen@anl.gov', url='https://github.com/ARM-DOE/ACT', - packages=find_packages(exclude=['docs', 'tests']), + packages=find_packages(exclude=['docs']), entry_points={'console_scripts': []}, include_package_data=True, package_data={'act': []},