diff --git a/.github/workflows/build-docs.yml b/.github/workflows/build-docs.yml new file mode 100644 index 0000000..dfcbce6 --- /dev/null +++ b/.github/workflows/build-docs.yml @@ -0,0 +1,46 @@ +# Build xlma-python documentation. + +name: Build Documentation + +on: [push, pull_request] + +permissions: + id-token: write + contents: read + pages: write + +jobs: + build: + runs-on: ubuntu-latest + steps: + - name: Checkout repository + uses: actions/checkout@v4 + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: 3.13 + - name: Install environment with micromamba + uses: mamba-org/setup-micromamba@v1 + with: + environment-file: docs/docs-env.yml + init-shell: bash + - name: Build with mkdocs + run: | + eval "$(micromamba shell hook --shell bash)" + micromamba activate pyxlma-docs + mkdocs build + - name: Upload static site artifact + id: upload + uses: actions/upload-pages-artifact@v3 + with: + path: ./site/ + deploy: + environment: + name: github-pages + url: ${{ steps.upload.outputs.page_url }} + runs-on: ubuntu-latest + needs: build + if: github.event_name == 'push' && github.ref == 'refs/heads/master' + steps: + - name: Deploy artifact to GH Pages + uses: actions/deploy-pages@v4 diff --git a/.github/workflows/python-package.yml b/.github/workflows/run-tests.yml similarity index 93% rename from .github/workflows/python-package.yml rename to .github/workflows/run-tests.yml index 6c1374e..da105a7 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/run-tests.yml @@ -1,13 +1,9 @@ # This workflow will install Python dependencies, run tests and lint with a variety of Python versions # For more information see: https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-python -name: Run pytests +name: Run tests -on: - push: - branches: [ "master" ] - pull_request: - branches: [ "master" ] +on: [push, pull_request] jobs: test: diff --git a/README.md b/README.md deleted file mode 100644 index e03ddf7..0000000 --- a/README.md +++ /dev/null @@ -1,103 +0,0 @@ -# xlma-python -A future, Python-based version of xlma? - -XLMA is a venerable IDL GUI program that diplays VHF Lightning Mapping Array data. It was originally written by New Mexico Tech in the late 1990s. The purpose of this repository is to collect ideas, sketches of code, and a community that could facilitate rewriting xlma in Python. - -Please use the issues tracker to discuss ideas and pull requests to contribute examples. - -# Installation -Clone this repostiory install with pip. - -``` -git clone https://github.com/deeplycloudy/xlma-python.git -cd xlma-python -pip install -e . -``` - -Then, copy the `XLMA_plots.ipynb` notebook to wherever you'd like and start changing files, dates and times to show data from your case of interest. There also a notebook showing how to do flash sorting and save a new NetCDF file with those data. - -# Dependencies -Required: - -- xarray (I/O requires the netcdf4 backend) -- pandas -- numpy - -Flash clustering: - -- scikit-learn -- scipy -- pyproj - -Plotting: - -- matplotlib -- cartopy -- metpy (optionally, for US county lines) - -GLM Plotting: -- glmtools (https://github.com/deeplycloudy/glmtools) - -Interactive: - -- jupyterlab (or, notebook) -- ipywidgets -- ipympl - -Building: - -- setuptools -- pytest-cov -- pytest-mpl -- lmatools (https://github.com/deeplycloudy/lmatools) -- ...and all of the above - -# Technical architecture - -We envision a two-part package that keeps a clean separation between the core data model, analysis, and display. XLMA utilized a large, global `state` structure that stored all data, as well as the current data selection corresponding to the view in the GUI. Analysis then operated on whatever data was in the current selection. - -## Data model and subsetting - -`xarray` is the obvious choice for the core data structure, because it supports multidimensional data with metadata for each variable, subsetting of all varaibles sharing a dimension, and fast indexing. Data can be easily saved and read from the NetCDF format, and converted from ASCII to `Dataset` using standard tools. - -## Analysis - -Some core features of LMA data analysis will be built in, TBD after surveying capabilities in XLMA. - -## Display - -Keeping the core data structure and selection operations separate from dislpay is good programming practice. It is doubly important in Python, where there is not one obvious solution for high performance *and* publication-quality graphics as in IDL. - -### Plotting library - -There are many options, so we want a design that: -1. Permits a GUI to provide the bounds of the current view (or a polygon lasso) to the data model, changing the subset -2. Allows many GUIs to read from the data model so that we maintain loose coupling as the GUI and plotting landscape evolves. - -- Matplotlib: publication quality figures, mature, and can plot everything, including weather symbols. Perhaps too slow for interactive, colormapped scatterplots of O(1e6) points, as is common for LMA data. -- Vispy: OpenGL, so can display plots quickly, but will require more low-level coding. See the SatPy/PyTroll/SIFT efforts. -- ipyvolume: originally written for volume rendering, but has grown to display points, etc. Browser/notebook based. -- Bokeh/holoviews/geoviews: browser-based; Javascript and WebGL front end with Python backend. Works well with Binder. -- Datashader might be useful as a method of data reduction prior to visualization even if we don't use Bokeh. -- Yt - written by the astronomy community in Python … is it fast enough? - -### GUI - -There is no obvious choice here, either. - -- Jupyter notebooks -- JupyterLab with custom plugin -- PyQT -- [Glue](https://github.com/glue-viz/glue/wiki/SciPy-2019-Tutorial-on-Multi-dimensional-Linked-Data-Exploration-with-Glue) - seems about 60% there out of the box?! - -# Prior art - -- [`lmatools`](https://github.com/deeplycloudy/lmatools/) - - Includes readers for LMA and NLDN data (using older methods from 2010) - - Flash sorting and gridding - - Has code for all necessary coordinate transforms. -- [`brawl4d`](https://github.com/deeplycloudy/brawl4d/) A working version of the basic GUI functionality of xlma. - - Based on matplotlib; plots can be drag-repositioned. Slow for large numbers of data points. - - Includes charge analyis that auto-saves to disk - - At one point, could display radar data underneath LMA data - - Built around a data pipeline, with a pool of data at the start, subsetting, projection, and finally display. diff --git a/README.md b/README.md new file mode 120000 index 0000000..ed9e8db --- /dev/null +++ b/README.md @@ -0,0 +1 @@ +./docs/index.md \ No newline at end of file diff --git a/docs/contributing.md b/docs/contributing.md new file mode 100644 index 0000000..401a395 --- /dev/null +++ b/docs/contributing.md @@ -0,0 +1,21 @@ +# Contributor's Guide +--- + +Discussion of contributions takes place in [GitHub's issues](https://github.com/deeplycloudy/xlma-python/issues) and [pull requests](https://github.com/deeplycloudy/xlma-python/pulls) sections of the xlma-python repository. + +To create a "developer install" of xlma-python: + +```sh +git clone https://github.com/deeplycloudy/xlma-python.git +cd xlma-python +pip install -e . +``` + +This will allow you to edit files in the xlma-python directory on your local machine, and the `pip install -e .` command will allow those changes to be applied to your python environment immediately. + + +xlma-python is built using [setuptools](https://setuptools.pypa.io/en/latest/). The library is entirely pure python, so no additional compliation steps are required for installation. + +Automated testing is included to prevent future code contributions from breaking compatibility with previous versions of the library. Tests are stored in the `tests/` directory. pytest is used for the basic test architecture, with the `pytest-mpl` plugin providing image differences. See the `tests/truth/images/` directory for baseline images. `pytest-cov` is used for coverage checking, results are uploaded to [CodeCov](https://codecov.io). + +Documentation is built by [mkdocstring-python](https://mkdocstrings.github.io/python/). See the `mkdocs.yml` file for configuration and the `docs/reference/` directory for page layout. \ No newline at end of file diff --git a/docs/docs-env.yml b/docs/docs-env.yml new file mode 100644 index 0000000..5a22375 --- /dev/null +++ b/docs/docs-env.yml @@ -0,0 +1,19 @@ +name: pyxlma-docs +channels: +- conda-forge +dependencies: +- metpy +- mkdocs-material +- mkdocstrings-python +- netcdf4 +- numpy +- pandas +- pip +- pyproj +- scikit-learn +- scipy +- xarray +- pip: + - git+https://github.com/deeplycloudy/lmatools + - git+https://github.com/deeplycloudy/glmtools@energy_threshold_util + - -e ../ \ No newline at end of file diff --git a/docs/index.md b/docs/index.md new file mode 100644 index 0000000..84c623c --- /dev/null +++ b/docs/index.md @@ -0,0 +1,56 @@ +logo + +# xlma-python +--- +*A future, Python-based version of XLMA?* + +## Motivation + +xlma-python is used for analysis and visualization of [Lightning Mapping Array](https://doi.org/10.1029/1999GL010856) (LMA) data. + +XLMA is a venerable IDL GUI program that diplays VHF Lightning Mapping Array data. It was originally written by New Mexico Tech in the late 1990s. This repository represents progress on community-contributed ideas and code to facilitate reimplementing workflows and features from XLMA in Python. + +## Features + +### Data model and subsetting + +xlma-python includes file IO support for the ASCII .dat format used by XLMA and `lma_analysis`. Datasets are read in to the widely-used [xarray.Dataset](https://docs.xarray.dev/en/stable/generated/xarray.Dataset.html#xarray.Dataset) object. VHF sources and network information are added as multidimensional data, with metadata for each variable. Data can be easily saved to and read from the NetCDF format, allowing for conversion from ASCII to `xarray.Dataset` to a NetCDF file. + +### Analysis + +xlma-python uses a DBSCAN clustering algorithm from [scikit-learn](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.DBSCAN.html) to assign groups of VHF sources recorded by the LMA into discrete lightning flashes, as well as tools to compute properties of the flashes. + +Subsetting of all varaibles sharing a dimension and fast indexing are supported. VHF Sources and clustered flashes can be filtered by their properties, and individual flashes/events can be interrogated easily. + +### Display + +Publication quality plots in the style of XLMA (with a plan view, latitude-altitude and longitude-altitude cross sections, altitude-histogram, and altitude-timeseries plots) are supported. These plots are generated with high level functions, using the popular [matplotlib](https://matplotlib.org/stable/index.html) plotting library. + +### Interactive + +There are several ongoing, parallel efforts to create an interactive analysis/display tool in python, using different GUI tools, each with their own pros/cons: + +- This repository contains an "interactive" module which utilizes matplotlib/jupyter/ipywidgets for interactivity. This was easy to implement with the already existing matplotlib plotting tools, and allows for "good enough" performance in jupyter notebooks, which are widely used in geosciences. + +- An older example using [Glue](https://github.com/glue-viz/glue) is also included in the examples directory. + +- Other options (for external projects): + - Vispy: OpenGL, so can display plots quickly, but will require more low-level coding. See the SatPy/PyTroll/SIFT efforts. + - ipyvolume: originally written for volume rendering, but has grown to display points, etc. Browser/notebook based. + - HoloViz Suite (holoviews/geoviews plotting libraries with bokeh backend, panel for dashboard layout): browser-based; Javascript and WebGL front end with Python backend. + - HoloViz's Datashader might be useful as a method of data reduction prior to visualization even if we don't use HoloViews/Bokeh. + - Yt: written by the astronomy community in Python ... is it fast enough? + - lonboard: python binding for deck.gl, very fast vector rendering of many points, but limited interactive tools and ONLY allows rendering geospatial data (ie, not cross sections or histograms) + - pydeck: lower-level deck.gl python binding which removes the geospatial data requirement but requires more work on the developer's side. + +## Prior art + +- [`lmatools`](https://github.com/deeplycloudy/lmatools/) + - Includes readers for LMA and NLDN data (using older methods from 2010) + - Flash sorting and gridding + - Has code for all necessary coordinate transforms. +- [`brawl4d`](https://github.com/deeplycloudy/brawl4d/) A working version of the basic GUI functionality of xlma. + - Based on matplotlib; plots can be drag-repositioned. Slow for large numbers of data points. + - Includes charge analyis that auto-saves to disk + - At one point, could display radar data underneath LMA data + - Built around a data pipeline, with a pool of data at the start, subsetting, projection, and finally display. diff --git a/docs/install.md b/docs/install.md new file mode 100644 index 0000000..34c3dae --- /dev/null +++ b/docs/install.md @@ -0,0 +1,46 @@ +## Installation + +To install the latest "bleeding edge" commit of xlma-python, you can pip install the git repo directly: +```sh +pip install git+https://github.com/deeplycloudy/xlma-python` +``` +Note that this does **NOT** automatically update in a conda environment when running `conda update --all`, you must fully reinstall the library to obtain new changes. + +For a development install see [contributing](../contributing/). + +## Dependencies +Required: + +- xarray (I/O requires the netcdf4 backend) +- pandas +- numpy + +Flash clustering: + +- scikit-learn +- scipy +- pyproj + +Plotting: + +- matplotlib +- cartopy +- metpy (optionally, for US county lines) + +GLM Plotting: + +- glmtools (https://github.com/deeplycloudy/glmtools) + +Interactive: + +- jupyterlab (or, notebook) +- ipywidgets +- ipympl + +Building: + +- setuptools +- pytest-cov +- pytest-mpl +- lmatools (https://github.com/deeplycloudy/lmatools) +- ...and all of the above \ No newline at end of file diff --git a/docs/reference/coords/index.md b/docs/reference/coords/index.md new file mode 100644 index 0000000..60e5a6e --- /dev/null +++ b/docs/reference/coords/index.md @@ -0,0 +1,13 @@ +# Coordinates +--- +The coords module handles conversions between coordinate systems. The `CoordinateSystem` class represents a generic transform. Subclasses of `CoordinateSystem` are specific coordinate systems that have implemented transformations. These subclasses are documented on the "Transforms" page, see the sidebar. + +There are a few useful tools related to coordinate transforms included as well. + +::: pyxlma.coords + options: + members: + - centers_to_edges + - centers_to_edges_2d + - semiaxes_to_invflattening + - CoordinateSystem \ No newline at end of file diff --git a/docs/reference/coords/transforms.md b/docs/reference/coords/transforms.md new file mode 100644 index 0000000..934d0b5 --- /dev/null +++ b/docs/reference/coords/transforms.md @@ -0,0 +1,12 @@ +# Transforms +--- +Classes representing implemented coordinate transformations. All classes use the [Earth-Centered, Earth-Fixed (ECEF)](https://en.wikipedia.org/wiki/Earth-centered,_Earth-fixed_coordinate_system) system as a common intermediary system, allowing for conversions between systems. ![Figure demonstrating Earth-Centered, Earth-Fixed](https://upload.wikimedia.org/wikipedia/commons/6/62/Ecef_coordinates.svg) + + +::: pyxlma.coords + options: + filters: + - "!^centers_to_edges$" + - "!^centers_to_edges_2d$" + - "!^semiaxes_to_invflattening$" + - "!^CoordinateSystem$" \ No newline at end of file diff --git a/docs/reference/lmalib/cf_netcdf.md b/docs/reference/lmalib/cf_netcdf.md new file mode 100644 index 0000000..bf9d07c --- /dev/null +++ b/docs/reference/lmalib/cf_netcdf.md @@ -0,0 +1,6 @@ +# Climate/Forecasting NetCDF + +This module is designed to provide a standardized way to name and attribute lightning data in NetCDF4 dataset, +akin to the [Climate and Forecasting](https://cfconventions.org) specification. These tools are primarily used internally for I/O operations, but are provided to the user for convenience. + +::: pyxlma.lmalib.io.cf_netcdf diff --git a/docs/reference/lmalib/flash.md b/docs/reference/lmalib/flash.md new file mode 100644 index 0000000..8fef466 --- /dev/null +++ b/docs/reference/lmalib/flash.md @@ -0,0 +1,7 @@ +# Flash processing + +::: pyxlma.lmalib.flash.cluster + +::: pyxlma.lmalib.flash.properties + + diff --git a/docs/reference/lmalib/grid.md b/docs/reference/lmalib/grid.md new file mode 100644 index 0000000..784113a --- /dev/null +++ b/docs/reference/lmalib/grid.md @@ -0,0 +1,3 @@ +# Gridding + +::: pyxlma.lmalib.grid \ No newline at end of file diff --git a/docs/reference/lmalib/lma_intercept_rhi.md b/docs/reference/lmalib/lma_intercept_rhi.md new file mode 100644 index 0000000..6fc39b4 --- /dev/null +++ b/docs/reference/lmalib/lma_intercept_rhi.md @@ -0,0 +1,3 @@ +# LMA with a radar RHI + +::: pyxlma.lmalib.lma_intercept_rhi \ No newline at end of file diff --git a/docs/reference/lmalib/read.md b/docs/reference/lmalib/read.md new file mode 100644 index 0000000..d9a6df6 --- /dev/null +++ b/docs/reference/lmalib/read.md @@ -0,0 +1,7 @@ +# Reading data + +::: pyxlma.lmalib.io.read + options: + filters: + - "!^to_dataset$" + - "!^lmafile$" \ No newline at end of file diff --git a/docs/reference/lmalib/traversal.md b/docs/reference/lmalib/traversal.md new file mode 100644 index 0000000..c16df93 --- /dev/null +++ b/docs/reference/lmalib/traversal.md @@ -0,0 +1,3 @@ +# Traversal + +::: pyxlma.lmalib.traversal \ No newline at end of file diff --git a/docs/reference/plot/base_plot.md b/docs/reference/plot/base_plot.md new file mode 100644 index 0000000..efb6ae1 --- /dev/null +++ b/docs/reference/plot/base_plot.md @@ -0,0 +1,3 @@ +# Base Plot + +::: pyxlma.plot.xlma_base_plot \ No newline at end of file diff --git a/docs/reference/plot/interactive.md b/docs/reference/plot/interactive.md new file mode 100644 index 0000000..02d7157 --- /dev/null +++ b/docs/reference/plot/interactive.md @@ -0,0 +1,3 @@ +# Jupyter Interactive + +::: pyxlma.plot.interactive \ No newline at end of file diff --git a/docs/reference/plot/plot_feature.md b/docs/reference/plot/plot_feature.md new file mode 100644 index 0000000..8386b39 --- /dev/null +++ b/docs/reference/plot/plot_feature.md @@ -0,0 +1,3 @@ +# LMA Feature Plotting + +::: pyxlma.plot.xlma_plot_feature \ No newline at end of file diff --git a/docs/reference/xarray/index.md b/docs/reference/xarray/index.md new file mode 100644 index 0000000..03c99df --- /dev/null +++ b/docs/reference/xarray/index.md @@ -0,0 +1,5 @@ +# [xarray](https://docs.xarray.dev/) utilities + +Some useful tools for handling xarray [datasets](https://docs.xarray.dev/en/stable/generated/xarray.Dataset.html), the primary data store used by pyxlma to store LMA data. + +::: pyxlma.xarray_util diff --git a/docs/xlma_logo.svg b/docs/xlma_logo.svg new file mode 100644 index 0000000..9ab0705 --- /dev/null +++ b/docs/xlma_logo.svg @@ -0,0 +1,332 @@ + + + +xlma diff --git a/mkdocs.yml b/mkdocs.yml new file mode 100644 index 0000000..0aba0e4 --- /dev/null +++ b/mkdocs.yml @@ -0,0 +1,76 @@ +site_name: XLMA Python +repo_name: deeplycloudy/xlma-python +repo_url: https://github.com/deeplycloudy/xlma-python + + +theme: + name: material + logo: xlma_logo.svg + favicon: xlma_logo.svg + icon: + repo: fontawesome/brands/github + palette: + + # Palette toggle for light mode + - scheme: default + toggle: + icon: material/brightness-7 + name: Switch to dark mode + + # Palette toggle for dark mode + - scheme: slate + toggle: + icon: material/brightness-4 + name: Switch to light mode + +plugins: + - search + - mkdocstrings: + enable_inventory: true + handlers: + python: + options: + docstring_style: numpy + show_root_heading: true + show_symbol_type_heading: true + show_symbol_type_toc: true + show_source: false + docstring_section_style: table + annotations_path: source + import: + - https://matplotlib.org/stable/objects.inv + - https://numpy.org/doc/stable/objects.inv + - https://scitools.org.uk/cartopy/docs/latest/objects.inv + - https://ipywidgets.readthedocs.io/en/stable/objects.inv + - https://docs.xarray.dev/en/stable/objects.inv + - https://pyproj4.github.io/pyproj/stable/objects.inv + - https://pandas.pydata.org/pandas-docs/stable/objects.inv + - https://docs.python.org/3/objects.inv + +markdown_extensions: + - pymdownx.highlight: + anchor_linenums: true + - pymdownx.superfences + +nav: + - Home: index.md + - Installation: install.md + - API Reference: + - coords: + - reference/coords/index.md + - reference/coords/transforms.md + - xarray: + - reference/xarray/index.md + - lmalib: + - reference/lmalib/flash.md + - reference/lmalib/cf_netcdf.md + - reference/lmalib/read.md + - reference/lmalib/lma_intercept_rhi.md + - reference/lmalib/grid.md + - reference/lmalib/traversal.md + - plot: + - reference/plot/base_plot.md + - reference/plot/plot_feature.md + - reference/plot/interactive.md + - Contributing: contributing.md + diff --git a/pyxlma/coords.py b/pyxlma/coords.py index d3f0b20..3efddd8 100644 --- a/pyxlma/coords.py +++ b/pyxlma/coords.py @@ -13,57 +13,91 @@ class CoordinateSystem(object): - """The abstract coordinate system handling provided here works as follows. + """Superclass representing a generic coordinate system. Subclasses represent specific coordinate systems. + + Each subclass coordinate system must be able to convert data to a common coordinate system, + which is chosen to be Earth-Centered, Earth-Fixed (ECEF) cartesian. - Each coordinate system must be able to convert data to a common coordinate system, which is chosen to be ECEF cartesian. - data -> common system - common system -> dislpay coordinates This is implemented by the fromECEF and toECEF methods in each coordinate system object. + User code is responsible for taking data in its native coord system, - transforming it using to/fromECEF using the a coord system appropriate to the data, and then - transforming that data to the final coordinate system using another coord system. - - Subclasses should maintain an attribute ERSxyz that can be used in - transformations to/from an ECEF cartesian system, e.g. - >>> self.ERSxyz = proj4.Proj(proj='geocent', ellps='WGS84', datum='WGS84') - >>> self.ERSlla = proj4.Proj(proj='latlong', ellps='WGS84', datum='WGS84') - >>> projectedData = proj4.Transformer.from_crs(self.ERSlla.crs, self.ERSxyz.crs).transform(lon, lat, alt) - The ECEF system has its origin at the center of the earth, with the +Z toward the north pole, - +X toward (lat=0, lon=0), and +Y right-handed orthogonal to +X, +Z - - Depends on pyproj, http://code.google.com/p/pyproj/ to handle the ugly details of + transforming it using to/fromECEF using the a coord system appropriate to the data, and then + transforming that data to the final coordinate system using another coord system. + + Subclasses should maintain an attribute ERSxyz that can be used in transformations to/from an ECEF cartesian system, e.g. + + ```py + import pyproj as proj4 + self.ERSxyz = proj4.Proj(proj='geocent', ellps='WGS84', datum='WGS84') + self.ERSlla = proj4.Proj(proj='latlong', ellps='WGS84', datum='WGS84') + projectedData = proj4.Transformer.from_crs(self.ERSlla.crs, self.ERSxyz.crs).transform(lon, lat, alt) + ``` + + The ECEF system has its origin at the center of the earth, with the +Z toward the north pole, +X toward (lat=0, lon=0), and +Y right-handed orthogonal to +X, +Z + + Depends on [pyproj](https://github.com/pyproj4/pyproj) to handle the ugly details of various map projections, geodetic transforms, etc. - "You can think of a coordinate system as being something like character encodings, - but messier, and without an obvious winner like UTF-8." - Django OSCON tutorial, 2007 + *"You can think of a coordinate system as being something like character encodings, + but messier, and without an obvious winner like UTF-8." - Django OSCON tutorial, 2007* http://toys.jacobian.org/presentations/2007/oscon/tutorial/ + + Notes + ----- + This class is not intended to be instantiated directly. Instead, use one of the subclasses documented on the "Transforms" page. """ # WGS84xyz = proj4.Proj(proj='geocent', ellps='WGS84', datum='WGS84') - def coordinates(): + def coordinates(self): """Return a tuple of standarized coordinate names""" - raise NotImplemented + raise NotImplementedError() def fromECEF(self, x, y, z): """Take ECEF x, y, z values and return x, y, z in the coordinate system defined by the object subclass""" - raise NotImplemented + raise NotImplementedError() def toECEF(self, x, y, z): """Take x, y, z in the coordinate system defined by the object subclass and return ECEF x, y, z""" - raise NotImplemented + raise NotImplementedError() class GeographicSystem(CoordinateSystem): - """ - Coordinate system defined on the surface of the earth using latitude, - longitude, and altitude, referenced by default to the WGS84 ellipse. + """Coordinate system defined using latitude, longitude, and altitude. + An ellipsoid is used to define the shape of the earth. Latitude and longitude represent the + location of a point on the ellipsoid, and altitude is the height above the ellipsoid. - Alternately, specify the ellipse shape using an ellipse known - to pyproj, or [NOT IMPLEMENTED] specify r_equator and r_pole directly. + Attributes + ---------- + ERSlla : pyproj.Proj + A Proj object representing the geographic coordinate system. + ERSxyz : pyproj.Proj + A Proj object representing the Earth-Centered, Earth-Fixed (ECEF) cartesian coordinate system. """ def __init__(self, ellipse='WGS84', datum='WGS84', r_equator=None, r_pole=None): + """Initialize a GeographicSystem object. + + Parameters + ---------- + + ellipse : str, default='WGS84' + Ellipse name recognized by pyproj. + + *Ignored if r_equator or r_pole are provided.* + datum : str, default='WGS84' + Datum name recognized by pyproj. + + *Ignored if r_equator or r_pole are provided.* + r_equator : float, optional + Semi-major axis of the ellipse in meters. + + *If only one of r_equator or r_pole is provided, the resulting ellipse is assumed to be spherical.* + r_pole : float, optional + Semi-minor axis of the ellipse in meters. + + *If only one of r_equator or r_pole is provided, the resulting ellipse is assumed to be spherical.* + """ if (r_equator is not None) | (r_pole is not None): if r_pole is None: r_pole=r_equator @@ -76,6 +110,26 @@ def __init__(self, ellipse='WGS84', datum='WGS84', self.ERSlla = proj4.Proj(proj='latlong', ellps=ellipse, datum=datum) self.ERSxyz = proj4.Proj(proj='geocent', ellps=ellipse, datum=datum) def toECEF(self, lon, lat, alt): + """Converts longitude, latitude, and altitude to Earth-Centered, Earth-Fixed (ECEF) X, Y, Z coordinates. + + Parameters + ---------- + lon : float or array_like + Longitude in decimal degrees East of the Prime Meridian. + lat : float or array_like + Latitude in decimal degrees North of the equator. + alt: float or array_like + Altitude in meters above the ellipsoid. + + Returns + ------- + X : float or array_like + ECEF X in meters from the center of the Earth. + Y : float or array_like + ECEF Y in meters from the center of the Earth. + Z : float or array_like + ECEF Z in meters from the center of the Earth. + """ lat = atleast_1d(lat) # proj doesn't like scalars lon = atleast_1d(lon) alt = atleast_1d(alt) @@ -87,6 +141,26 @@ def toECEF(self, lon, lat, alt): return projectedData[0,:], projectedData[1,:], projectedData[2,:] def fromECEF(self, x, y, z): + """Converts Earth-Centered, Earth-Fixed (ECEF) X, Y, Z to longitude, latitude, and altitude coordinates. + + Parameters + ---------- + x : float or array_like + ECEF X in meters from the center of the Earth. + y : float or array_like + ECEF Y in meters from the center of the Earth. + z : float or array_like + ECEF Z in meters from the center of the Earth. + + Returns + ------- + lon : float or array_like + Longitude in decimal degrees East of the Prime Meridian. + lat : float or array_like + Latitude in decimal degrees North of the equator. + alt: float or array_like + Altitude in meters above the ellipsoid. + """ x = atleast_1d(x) # proj doesn't like scalars y = atleast_1d(y) z = atleast_1d(z) @@ -99,11 +173,49 @@ def fromECEF(self, x, y, z): class MapProjection(CoordinateSystem): - """Map projection coordinate system. Wraps proj4, and uses its projecion names. Defaults to - equidistant cylindrical projection + """Coordinate system defined using meters x, y, z in a specified map projection. + Wraps pyproj, and uses its projecion names. Converts location in any map projection to ECEF, and vice versa. + + Attributes + ---------- + ERSxyz : pyproj.Proj + A Proj object representing the Earth-Centered, Earth-Fixed (ECEF) cartesian coordinate system. + projection : pyproj.Proj + A Proj object representing the map projection. + ctrLat : float + Latitude of the center of the map projection in decimal degrees North of the equator, if required for the projection. + ctrLon : float + Longitude of the center of the map projection in decimal degrees East of the Prime Meridian, if required for the projection. + ctrAlt : float + Altitude of the center of the map projection in meters. + geoCS : GeographicSystem + GeographicSystem object used to convert the map projection's center to ECEF. + cx : float + X coordinate of the map projection's center in meters. + cy : float + Y coordinate of the map projection's center in meters. + cz : float + Z coordinate of the map projection's center in meters. """ def __init__(self, projection='eqc', ctrLat=None, ctrLon=None, ellipse='WGS84', datum='WGS84', **kwargs): + """Initialize a MapProjection object. + + Parameters + ---------- + projection : str, default='eqc' + Projection name recognized by pyproj. Defaults to 'eqc' (equidistant cylindrical). + ctrLat : float, optional + Latitude of the center of the map projection in decimal degrees North of the equator, if required for the projection. + ctrLon : float, optional + Longitude of the center of the map projection in decimal degrees East of the Prime Meridian, if required for the projection. + ellipse : str, default='WGS84' + Ellipse name recognized by pyproj. + datum : str, default='WGS84' + Datum name recognized by pyproj. + **kwargs + Additional keyword arguments passed to pyproj.Proj() + """ self.ERSxyz = proj4.Proj(proj='geocent', ellps=ellipse, datum=datum) self.projection = proj4.Proj(proj=projection, ellps=ellipse, datum=datum, **kwargs) self.ctrLat=ctrLat @@ -114,6 +226,17 @@ def __init__(self, projection='eqc', ctrLat=None, ctrLon=None, ellipse='WGS84', self.cx, self.cy, self.cz = self.ctrPosition() def ctrPosition(self): + """Get the map projection's center position as projected in the specified map projection. + + Returns + ------- + cx : float + X coordinate of the map projection's center in meters. + cy : float + Y coordinate of the map projection's center in meters. + cz : float + Z coordinate of the map projection's center in meters. + """ if (self.ctrLat != None) & (self.ctrLon != None): ex, ey, ez = self.geoCS.toECEF(self.ctrLon, self.ctrLat, self.ctrAlt) cx, cy, cz = self.fromECEF(ex, ey, ez) @@ -122,6 +245,26 @@ def ctrPosition(self): return cx, cy, cz def toECEF(self, x, y, z): + """Converts x, y, z meters in the map projection to Earth-Centered, Earth-Fixed (ECEF) X, Y, Z coordinates. + + Parameters + ---------- + x : float or array_like + x position in the map projection in meters. + y : float or array_like + y position in the map projection in meters. + z: float or array_like + z position in the map projection in meters. + + Returns + ------- + X : float or array_like + ECEF X in meters from the center of the Earth. + Y : float or array_like + ECEF Y in meters from the center of the Earth. + Z : float or array_like + ECEF Z in meters from the center of the Earth. + """ x += self.cx y += self.cy z += self.cz @@ -133,6 +276,26 @@ def toECEF(self, x, y, z): return px, py, pz def fromECEF(self, x, y, z): + """Converts x, y, z meters in the map projection to Earth-Centered, Earth-Fixed (ECEF) X, Y, Z coordinates. + + Parameters + ---------- + x : float or array_like + x position in the map projection in meters. + y : float or array_like + y position in the map projection in meters. + z: float or array_like + z position in the map projection in meters. + + Returns + ------- + X : float or array_like + ECEF X in meters from the center of the Earth. + Y : float or array_like + ECEF Y in meters from the center of the Earth. + Z : float or array_like + ECEF Z in meters from the center of the Earth. + """ projectedData = array(proj4.Transformer.from_crs(self.ERSxyz.crs, self.projection.crs).transform(x, y, z)) if len(projectedData.shape) == 1: px, py, pz = projectedData[0], projectedData[1], projectedData[2] @@ -141,35 +304,60 @@ def fromECEF(self, x, y, z): return px-self.cx, py-self.cy, pz-self.cz class PixelGrid(CoordinateSystem): + """Coordinate system defined using arbitrary pixel coordinates in a 2D pixel array. + + Attributes + ---------- + geosys : GeographicSystem + GeographicSystem object used to convert pixel coordinates to ECEF. + lookup : object + Lookup object used to find the nearest pixel to a given lat/lon. See `__init__` for more details. + x : list or array_like + 1D integer array of pixel row IDs. + y : list or array_like + 1D integer array of pixel column IDs. + lons : array_like + 2D array of longitudes of pixel centers. + lats : array_like + 2D array of latitudes of pixel centers. + alts : array_like + 2D array of altitudes of pixel centers. + """ def __init__(self, lons, lats, lookup, x, y, alts=None, geosys=None): - """ - Coordinate system for arbitrary pixel coordinates in a 2D pixel array. - Arguments: - lons: 2D array of longitudes of pixel centers - lats: 2D array of longitudes of pixel centers - alts: 2D array of longitudes of pixel centers. If None, zeros are assumed. - Each array is of shape (nx, ny) with pixel coordinate (x=0, y=0) - corresponding to grid index [0, 0] - - lookup is an object with a method 'query' that accepts a single argument, - a (N,2) array of lats, lons and returns pixel IDs that can be used to - index lons and lats, as well as the distances between the pixel centers - and the queried locations. X and Y flattened arrays of pixel coordinates - that align with indices of the flattened lon and lat arrays used to - create the lookup table. - >>> test_events = np.vstack([(-101.5, 33.5), (-102.8, 32.5), (-102.81,32.5)]) - >>> distances, idx = lookup.query(test_events) - >>> loni, lati = lons[X[idx], Y[idx]], lats[X[idx], Y[idx]] - An instance of sklearn.neighbors.KDTree is one such lookup. - - If geosys is provided, it should be an instance of GeographicSystem; - otherwise a GeographicSystem instance with default arguments is created. + """Initialize a PixelGrid object. + Parameters + ---------- + lons : array_like + 2D array of longitudes of pixel centers. + lats : array_like + 2D array of latitudes of pixel centers. + lookup : object + Object with instance method `query` that accepts a single argument, a (N,2) array of lats, lons and + returns a tuple of distances to the nearest pixel centers and the pixel ID of the nearest pixel to each + requested lat/lon. A [`sklearn.neighbors.KDTree`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.cKDTree.html) is the intended use of this argument, but any class with a `query` method is accepted. + + Example of a valid lookup object: + + test_events = np.vstack([(-101.5, 33.5), (-102.8, 32.5), (-102.81,32.5)]) + distances, idx = lookup.query(test_events) + loni, lati = lons[X[idx], Y[idx]], lats[X[idx], Y[idx]] + + x : list or array_like + 1D integer array of pixel row IDs + y : list or array_like + 1D integer array of pixel column IDs + alts : array_like, optional + 2D array of altitudes of pixel centers. If None, zeros are assumed. + geosys : GeographicSystem, optional + GeographicSystem object used to convert pixel coordinates to ECEF. If None, a GeographicSystem instance with default arguments is created. + + Notes + ----- When converting toECEF, which accepts pixel coordinates, the z pixel coordinate is ignored, as it has no meaning. When converting fromECEF, zeros in the shape of x are returned as the z coordinate. - """ if geosys is None: self.geosys = GeographicSystem() @@ -184,7 +372,27 @@ def __init__(self, lons, lats, lookup, x, y, alts=None, geosys=None): alts = zeros_like(lons) self.alts = alts - def toECEF(self, x, y, z): + def toECEF(self, x, y, z=None): + """Converts x, y pixel IDs to Earth-Centered, Earth-Fixed (ECEF) X, Y, Z coordinates. + + Parameters + ---------- + x : float or array_like + row ID of the pixel in the pixel grid. + y : float or array_like + column ID of the pixel in the pixel grid. + z : object, optional + unused. If provided, it is ignored. + + Returns + ------- + X : float or array_like + ECEF X in meters from the center of the Earth. + Y : float or array_like + ECEF Y in meters from the center of the Earth. + Z : float or array_like + ECEF Z in meters from the center of the Earth. + """ x = x.astype('int64') y = y.astype('int64') lons = self.lons[x, y] @@ -193,6 +401,26 @@ def toECEF(self, x, y, z): return self.geosys.toECEF(lons, lats, alts) def fromECEF(self, x, y, z): + """Converts Earth-Centered, Earth-Fixed (ECEF) X, Y, Z to x, y pixel ID coordinates. + + Parameters + ---------- + x : float or array_like + ECEF X in meters from the center of the Earth. + y : float or array_like + ECEF Y in meters from the center of the Earth. + z : float or array_like + ECEF Z in meters from the center of the Earth. + + Returns + ------- + x : array_like + row ID of the pixel in the pixel grid. + y : array_like + column ID of the pixel in the pixel grid. + z : array_like + Zeros array in the shape of x. + """ lons, lats, alts = self.geosys.fromECEF(x, y, z) locs = vstack((lons.flatten(), lats.flatten())).T if locs.shape[0] > 0: @@ -204,28 +432,38 @@ def fromECEF(self, x, y, z): return x, y, zeros_like(x) class GeostationaryFixedGridSystem(CoordinateSystem): + """Coordinate system defined using scan angles from the perspective of a geostationary satellite. + + The pixel locations are a 2D grid of scan angles (in radians) from the perspective of a + geostationary satellite above an arbitrary ellipsoid. + Attributes + ---------- + ECEFxyz : pyproj.Proj + A Proj object representing the Earth-Centered, Earth-Fixed (ECEF) cartesian coordinate system. + fixedgrid : pyproj.Proj + A Proj object representing the geostationary fixed grid coordinate system. + h : float + Height of the satellite in meters above the specified ellipsoid. + """ def __init__(self, subsat_lon=0.0, subsat_lat=0.0, sweep_axis='y', sat_ecef_height=35785831.0, ellipse='WGS84'): - """ - Coordinate system representing a grid of scan angles from the perspective of a - geostationary satellite above an arbitrary ellipsoid. Fixed grid coordinates are - in radians. + """Initialize a GeostationaryFixedGridSystem object. Parameters ---------- - subsat_lon : float + subsat_lon : float, default=0.0 Longitude of the subsatellite point in degrees. - subsat_lat : float + subsat_lat : float, default=0.0 Latitude of the subsatellite point in degrees. - sweep_axis : str + sweep_axis : str, default='y' Axis along which the satellite sweeps. 'x' or 'y'. Use 'x' for GOES - and 'y' (default) for EUMETSAT. - sat_ecef_height : float - Height of the satellite in meters above the specified ellipsoid. - ellipse : str or list - A string representing a known ellipse to pyproj, or a list of [a, b] (semi-major + and 'y' for EUMETSAT. + sat_ecef_height : float, default=35785831.0 + Height of the satellite in meters above the specified ellipsoid. Defaults to the height of the GOES satellite. + ellipse : str or iterable, default='WGS84' + A string representing a known ellipse to pyproj, or iterable of [a, b] (semi-major and semi-minor axes) of the ellipse. Default is 'WGS84'. """ if type(ellipse) == str: @@ -234,7 +472,7 @@ def __init__(self, subsat_lon=0.0, subsat_lat=0.0, sweep_axis='y', rf = semiaxes_to_invflattening(ellipse[0], ellipse[1]) ellipse_args = {'a': ellipse[0], 'rf': rf} else: - raise ValueError("Ellipse must be a string or list of [a, b].") + raise ValueError("Ellipse must be a string or iterable of [a, b].") self.ECEFxyz = proj4.Proj(proj='geocent', **ellipse_args) self.fixedgrid = proj4.Proj(proj='geos', lon_0=subsat_lon, lat_0=subsat_lat, h=sat_ecef_height, x_0=0.0, y_0=0.0, @@ -242,10 +480,50 @@ def __init__(self, subsat_lon=0.0, subsat_lat=0.0, sweep_axis='y', self.h=sat_ecef_height def toECEF(self, x, y, z): + """Converts x, y, z satellite scan angles to Earth-Centered, Earth-Fixed (ECEF) X, Y, Z coordinates. + + Parameters + ---------- + x : float or array_like + horizontal scan angle in radians from the perspective of the satellite. + y : float or array_like + vertical scan angle in radians from the perspective of the satellite. + z : float or array_like + altitude above the ellipsoid expressed as a fraction of the satellite's height above the ellipsoid. + + Returns + ------- + X : float or array_like + ECEF X in meters from the center of the Earth. + Y : float or array_like + ECEF Y in meters from the center of the Earth. + Z : float or array_like + ECEF Z in meters from the center of the Earth. + """ X, Y, Z = x*self.h, y*self.h, z*self.h return proj4.Transformer.from_crs(self.fixedgrid.crs, self.ECEFxyz.crs).transform(X, Y, Z) def fromECEF(self, x, y, z): + """Converts Earth-Centered, Earth-Fixed (ECEF) X, Y, Z to longitude, latitude, and altitude coordinates. + + Parameters + ---------- + x : float or array_like + ECEF X in meters from the center of the Earth. + y : float or array_like + ECEF Y in meters from the center of the Earth. + z : float or array_like + ECEF Z in meters from the center of the Earth. + + Returns + ------- + x : float or array_like + horizontal scan angle in radians from the perspective of the satellite. + y : float or array_like + vertical scan angle in radians from the perspective of the satellite. + z : float or array_like + altitude above the ellipsoid expressed as a fraction of the satellite's height above the ellipsoid. + """ X, Y, Z = proj4.Transformer.from_crs(self.ECEFxyz.crs, self.fixedgrid.crs).transform(x, y, z) return X/self.h, Y/self.h, Z/self.h @@ -259,13 +537,57 @@ def fromECEF(self, x, y, z): # return px, py, z class RadarCoordinateSystem(CoordinateSystem): - """ - Converts spherical (range, az, el) radar coordinates to lat/lon/alt, and then to ECEF. + """Coordinate system defined using the range, azimuth, and elevation angles from a radar. + + Locations are defined using the latitude, longitude, and altitude of the radar and the azimuth and elevation angles of the radar beam. - An earth's effective radius of 4/3 is assumed to correct for atmospheric refraction. + Attributes + ---------- + ctrLat : float + Latitude of the radar in decimal degrees North of the equator. + ctrLon : float + Longitude of the radar in decimal degrees East of the Prime Meridian. + ctrAlt : float + Altitude of the radar in meters above sea level. + datum : str + Datum name recognized by pyproj. + ellps : str + Ellipse name recognized by pyproj. + lla : pyproj.Proj + A Proj object representing the geographic coordinate system. + xyz : pyproj.Proj + A Proj object representing the Earth-Centered, Earth-Fixed (ECEF) cartesian coordinate system. + Requator : float + Equatorial radius of the earth in meters. + Rpolar : float + Polar radius of the earth in meters. + flattening : float + Ellipsoid flattening parameter. + eccen : float + First eccentricity squared. + effectiveRadiusMultiplier : float + Multiplier to scale the earth's radius to account for the beam bending due to atmospheric refraction. """ - def __init__(self, ctrLat, ctrLon, ctrAlt, datum='WGS84', ellps='WGS84', effectiveRadiusMultiplier=4./3.): + def __init__(self, ctrLat, ctrLon, ctrAlt, datum='WGS84', ellps='WGS84', effectiveRadiusMultiplier=4/3): + """Initialize a RadarCoordinateSystem object. + + Parameters + ---------- + + ctrLat : float + Latitude of the radar in decimal degrees North of the equator. + ctrLon : float + Longitude of the radar in decimal degrees East of the Prime Meridian. + ctrAlt : float + Altitude of the radar in meters above sea level. + datum : str, default='WGS84' + Datum name recognized by pyproj. + ellps : str, default='WGS84' + Ellipse name recognized by pyproj. + effectiveRadiusMultiplier : float, default=4/3 + Multiplier to scale the earth's radius to account for the beam bending due to atmospheric refraction. + """ self.ctrLat = float(ctrLat) self.ctrLon = float(ctrLon) self.ctrAlt = float(ctrAlt) @@ -283,9 +605,29 @@ def __init__(self, ctrLat, ctrLon, ctrAlt, datum='WGS84', ellps='WGS84', effecti self.effectiveRadiusMultiplier = effectiveRadiusMultiplier def getGroundRangeHeight(self, r, elevationAngle): - """Convert slant range (along the beam) and elevation angle into - ground range (great circle distance) and height above the earth's surface - Follows Doviak and Zrnic 1993, eq. 2.28.""" + """Convert slant range (along the beam) and elevation angle into ground range and height. + Ground range given in great circle distance and height above the surface of the ellipsoid. + Follows [Doviak and Zrnic 1993, eq. 2.28.](https://doi.org/10.1016/C2009-0-22358-0) + + Parameters + ---------- + r : float or array_like + slant range in meters. + elevationAngle : float or array_like + elevation angle in degrees above the horizon. + + Returns + ------- + s : float or array_like + Ground range (great circle distance) in meters. + z : float or array_like + Height above the surface of the ellipsoid in meters. + + Warning + ------- + By default, an imaginary earth radius of 4/3 the actual earth radius is assumed to correct for atmospheric refraction. + This is a common assumption in radar meteorology, but may not be accurate in all cases. + """ #Double precison arithmetic is crucial to proper operation. lat = self.ctrLat * pi / 180.0 @@ -311,9 +653,23 @@ def getGroundRangeHeight(self, r, elevationAngle): return s, h def getSlantRangeElevation(self, groundRange, z): - """Convert ground range (great circle distance) and height above - the earth's surface to slant range (along the beam) and elevation angle. - Follows Doviak and Zrnic 1993, eq. 2.28""" + """Convert ground range (great circle distance) and height above the earth's surface to slant range (along the beam) and elevation angle. + Follows [Doviak and Zrnic 1993, eq. 2.28.](https://doi.org/10.1016/C2009-0-22358-0) + + Parameters + ---------- + groundRange : float or array_like + Ground range (great circle distance) in meters. + z : float or array_like + Height above the surface of the ellipsoid in meters. + + Returns + ------- + r : float or array_like + Slant range in meters. + el : float or array_like + Elevation angle in degrees above the horizon. + """ lat = self.ctrLat * pi / 180.0 @@ -344,8 +700,26 @@ def getSlantRangeElevation(self, groundRange, z): return r, el def toLonLatAlt(self, r, az, el): - """Convert slant range r, azimuth az, and elevation el to ECEF system""" - geoSys = GeographicSystem() + """Convert slant range r, azimuth az, and elevation el to latitude, longitude, altitude coordiantes. + + Parameters + ---------- + r : float or array_like + Slant range in meters. + az : float or array_like + Azimuth angle in degrees clockwise from North. + el : float or array_like + Elevation angle in degrees above the horizon. + + Returns + ------- + lon : float or array_like + Longitude in decimal degrees East of the Prime Meridian. + lat : float or array_like + Latitude in decimal degrees North of the equator. + z : float or array_like + Altitude in meters above the surface of the ellipsoid. + """ geodetic = proj4.Geod(ellps=self.ellps) try: @@ -358,12 +732,51 @@ def toLonLatAlt(self, r, az, el): return lon, lat, z def toECEF(self, r, az, el): + """Converts range, azimuth, and elevation to Earth-Centered, Earth-Fixed (ECEF) X, Y, Z coordinates. + + Parameters + ---------- + r : float or array_like + Distance in meters along the radar beam from the target to the radar. + az : float or array_like + Azimuth angle of the target in degrees clockwise from North. + el: float or array_like + Elevation angle of the target in degrees above the horizon. + + Returns + ------- + X : float or array_like + ECEF X in meters from the center of the Earth. + Y : float or array_like + ECEF Y in meters from the center of the Earth. + Z : float or array_like + ECEF Z in meters from the center of the Earth. + """ geoSys = GeographicSystem() lon, lat, z = self.toLonLatAlt(r, az, el) return geoSys.toECEF(lon, lat, z.ravel()) def fromECEF(self, x, y, z): - """Convert ECEF system to slant range r, azimuth az, and elevation el""" + """Converts Earth-Centered, Earth-Fixed (ECEF) X, Y, Z to longitude, latitude, and altitude coordinates. + + Parameters + ---------- + x : float or array_like + ECEF X in meters from the center of the Earth. + y : float or array_like + ECEF Y in meters from the center of the Earth. + z : float or array_like + ECEF Z in meters from the center of the Earth. + + Returns + ------- + r : float or array_like + Distance in meters along the radar beam from the target to the radar. + az : float or array_like + Azimuth angle of the target in degrees clockwise from North. + el: float or array_like + Elevation angle of the target in degrees above the horizon. + """ # x = np.atleast1d(x) geoSys = GeographicSystem() geodetic = proj4.Geod(ellps=self.ellps) @@ -384,12 +797,35 @@ def fromECEF(self, x, y, z): return r, az, el -class TangentPlaneCartesianSystem(object): - """ TODO: This function needs to be updated to inherit from CoordinateSystem +class TangentPlaneCartesianSystem(CoordinateSystem): + """Coordinate system defined by a meters relative to plane tangent to the earth at a specified location. + Attributes + ---------- + ctrLat : float + Latitude of the center of the local tangent plane in decimal degrees North of the equator. + ctrLon : float + Longitude of the center of the local tangent plane in decimal degrees East of the Prime Meridian. + ctrAlt : float + Altitude of the center of the local tangent plane in meters above the ellipsoid. + centerECEF : numpy.ndarray + ECEF X, Y, Z coordinates of the center of the local tangent plane. + TransformToLocal : numpy.ndarray + Rotation matrix to convert from ECEF to local tangent plane coordinates. """ def __init__(self, ctrLat=0.0, ctrLon=0.0, ctrAlt=0.0): + """Initialize a TangentPlaneCartesianSystem object. + + Parameters + ---------- + ctrLat : float, default=0.0 + Latitude of the center of the local tangent plane in decimal degrees North of the equator. + ctrLon : float, default=0.0 + Longitude of the center of the local tangent plane in decimal degrees East of the Prime Meridian. + ctrAlt : float, default=0.0 + Altitude of the center of the local tangent plane in meters above the ellipsoid. + """ self.ctrLat = float(ctrLat) self.ctrLon = float(ctrLon) self.ctrAlt = float(ctrAlt) @@ -456,7 +892,32 @@ def __init__(self, ctrLat=0.0, ctrLon=0.0, ctrAlt=0.0): [z1, z2, z3]]).squeeze() def fromECEF(self, x, y, z): - """ Transforms 1D arrays of ECEF x, y, z to the local tangent plane system""" + """Converts x, y, z meters in the local tangent plane to Earth-Centered, Earth-Fixed (ECEF) X, Y, Z coordinates. + + Parameters + ---------- + x : float or array_like + x position in meters East of the tangent plane center. + y : float or array_like + y position in meters North of the tangent plane center. + z: float or array_like + z position in meters above the tangent plane center. + + Returns + ------- + X : float or array_like + ECEF X in meters from the center of the Earth. + Y : float or array_like + ECEF Y in meters from the center of the Earth. + Z : float or array_like + ECEF Z in meters from the center of the Earth. + + Warning + ------- + + The z coordinate input is **NOT** the altitude above the ellipsoid. It is the z position in the local tangent plane. + Due to the curvature of the Earth, the TPCS z position and altitude difference increases with distance from the center of the TPCS. + """ data = vstack((x, y, z)) tpXYZ = self.toLocal(data) if len(tpXYZ.shape) == 1: @@ -466,7 +927,35 @@ def fromECEF(self, x, y, z): return tpX, tpY, tpZ def toECEF(self, x, y, z): - """ Transforms 1D arrays of x, y, z in the local tangent plane system to ECEF""" + """Converts Earth-Centered, Earth-Fixed (ECEF) X, Y, Z to x, y, z meters in the local tangent plane. + + Parameters + ---------- + x : float or array_like + ECEF X in meters from the center of the Earth. + y : float or array_like + ECEF Y in meters from the center of the Earth. + z : float or array_like + ECEF Z in meters from the center of the Earth. + + Returns + ------- + x : float or array_like + x position in meters East of the tangent plane center. + y : float or array_like + y position in meters North of the tangent plane center. + z: float or array_like + z position in meters above the tangent plane center. + + Warnings + -------- + The x and z output coordinates are not the great circle distance (the distance along the surface of the Earth). Distances between points should be compared in their ECEF coordinates. + + Similarly, the z coordinate output is **NOT** the altitude above the ellipsoid. It is the z position in the local tangent plane. + Due to the curvature of the Earth, the TPCS z position and altitude difference increases with distance from the center of the TPCS. + + If you want to find the altitude of a point above the ellipsoid, use the `GeographicSystem` class. + """ data = vstack((x, y, z)) ecXYZ = self.fromLocal(data) if len(ecXYZ.shape) == 1: @@ -476,52 +965,80 @@ def toECEF(self, x, y, z): return ecX, ecY, ecZ def toLocal(self, data): - """Transforms 3xN array of data (position vectors) in the ECEF system to - the local tangent plane cartesian system. + """Transforms 3xN array of ECEF X, Y, Z coordinates to the local tangent plane cartesian system. - Returns another 3xN array. + Parameters + ---------- + data : array_like + (3, N) array of data (position vectors) in the ECEF system, representing X, Y, Z meters from the center of the Earth. + + Returns + ------- + local_data : array_like + (3, N) array of data (position vectors) in the local tangent plane cartesian system, representing x, y, z meters from the center of the tangent plane. """ - return array( [ dot(self.TransformToLocal, (v-self.centerECEF)[:,None]) + local_data = array( [ dot(self.TransformToLocal, (v-self.centerECEF)[:,None]) for v in data[0:3,:].transpose()] ).squeeze().transpose() + return local_data def fromLocal(self, data): - """Transforms 3xN array of data (position vectors) in the local tangent - plane cartesian system to the ECEF system. + """Transforms 3xN array of ECEF X, Y, Z coordinates to the local tangent plane cartesian system. - Returns another 3xN array. + Parameters + ---------- + data : array_like + (3, N) array of data (position vectors) in the ECEF system, representing X, Y, Z meters from the center of the Earth. + + Returns + ------- + ecef_data : array_like + (3, N) array of data (position vectors) in the local tangent plane cartesian system, representing x, y, z meters from the center of the tangent plane. """ #Transform from local to ECEF uses transpose of the TransformToLocal matrix - return array( [ (dot(self.TransformToLocal.transpose(), v) + self.centerECEF) + ecef_data = array( [ (dot(self.TransformToLocal.transpose(), v) + self.centerECEF) for v in data[0:3,:].transpose()] ).squeeze().transpose() + return ecef_data def semiaxes_to_invflattening(semimajor, semiminor): - """ Calculate the inverse flattening from the semi-major - and semi-minor axes of an ellipse""" + """ Calculate the inverse flattening from the semi-major and semi-minor axes of an ellipse" + + Parameters + ---------- + semimajor : float + Semi-major axis of the ellipse + semiminor : float + Semi-minor axis of the ellipse + + Returns + ------- + rf : float + Inverse flattening of the ellipse + """ rf = semimajor/(semimajor-semiminor) return rf def centers_to_edges(x): - """ - Create an array of length N+1 edge locations from an - array of lenght N grid center locations. - - In the interior, the edge positions set to the midpoints - of the values in x. For the outermost edges, half the - closest dx is assumed to apply. + """ Create an array of length N+1 edge locations from an array of lenght N grid center locations. Parameters ---------- - x : array, shape (N) - Locations of the centers + x : array_like + (N,) array of locations of the centers Returns ------- - xedge : array, shape (N+1,M+1) - + xedge : numpy.ndarray + (N+1,) array of locations of the edges + + Notes + ----- + In the interior, the edge positions set to the midpoints + of the values in x. For the outermost edges, half the + closest dx is assumed to apply. """ xedge=zeros(x.shape[0]+1) xedge[1:-1] = (x[:-1] + x[1:])/2.0 @@ -531,26 +1048,28 @@ def centers_to_edges(x): def centers_to_edges_2d(x): - """ - Create a (N+1, M+1) array of edge locations from a - (N, M) array of grid center locations. + """Create a (N+1)x(M+1) array of edge locations from a + NxM array of grid center locations. - In the interior, the edge positions set to the midpoints - of the values in x. For the outermost edges, half the - closest dx is assumed to apply. This matters for polar - meshes, where one edge of the grid becomes a point at the - polar coordinate origin; dx/2 is a half-hearted way of - trying to prevent negative ranges. Parameters ---------- - x : array, shape (N,M) - Locations of the centers + x : array_like + (N,M) array locations of the centers. Returns ------- - xedge : array, shape (N+1,M+1) + xedge : array_like + (N+1,M+1) array of locations of the edges. + Notes + ----- + In the interior, the edge positions set to the midpoints + of the values in x. For the outermost edges, half the + closest dx is assumed to apply. This matters for polar + meshes, where one edge of the grid becomes a point at the + polar coordinate origin; dx/2 is a half-hearted way of + trying to prevent negative ranges. """ xedge = zeros((x.shape[0]+1,x.shape[1]+1)) # interior is a simple average of four adjacent centers diff --git a/pyxlma/lmalib/flash/cluster.py b/pyxlma/lmalib/flash/cluster.py index 871183a..0fc96f7 100644 --- a/pyxlma/lmalib/flash/cluster.py +++ b/pyxlma/lmalib/flash/cluster.py @@ -6,14 +6,27 @@ import pyxlma.lmalib.io.cf_netcdf as cf_netcdf def cluster_dbscan(X, Y, Z, T, min_points=1): - """ Identify clusters in spatiotemporal data X, Y, Z, T. + """Identify clusters in spatiotemporal data X, Y, Z, T. - min_points is used as min_samples in the call to DBSCAN. + Parameters + ---------- + X : array_like + The x coordinate of the data points. + Y : array_like + The y coordinate of the data points. + Z : array_like + The z coordinate of the data points. + T : array_like + The time coordinate of the data points. + min_points : int, default=1 + Used as the min_samples parameter in the call to [DBSCAN](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.DBSCAN.html). - Returns an unsigned 64 bit integer array of cluster labels. - Noise points identified by DBSCAN are assigned the maximum value - that can be represented by uint64, i.e., np.iinfo(np.uint64).max or - 18446744073709551615. + Returns + ------- + labels : numpy.ndarray + an unsigned 64 bit integer array of cluster labels. + Noise points identified by DBSCAN are assigned the maximum value that can be represented by uint64, + i.e., np.iinfo(np.uint64).max or 18446744073709551615. """ from sklearn.cluster import DBSCAN coords = np.vstack((X, Y, Z, T)).T @@ -25,12 +38,25 @@ def cluster_dbscan(X, Y, Z, T, min_points=1): return labels def cluster_flashes(events, distance=3000.0, time=0.15): - """ + """Cluster LMA VHF sources into flashes. + + Parameters + ---------- + events : xarray.Dataset + LMA dataset with event position and time and network center position. + distance : float, default=3000.0 + Spatial separation in meters. Used for normalization of space. + time : float, default=0.15 + Temporal separation in seconds. Used for normalization of time. - events: xarray dataset conforming to the pyxlma CF NetCDF format + Returns + ------- + ds : xarray.Dataset + LMA dataset with added flash_id and event_parent_flash_id variables. - distance: spatial separation in meters. Used for normalization of - time: number + Notes + ----- + Additional data variables for flash properties are created, but are filled with NaN. To compute these properties, use the `pyxlma.lmalib.flash.properties.flash_stats` function. """ geoCS = GeographicSystem() diff --git a/pyxlma/lmalib/flash/properties.py b/pyxlma/lmalib/flash/properties.py index eec9384..cd25c82 100644 --- a/pyxlma/lmalib/flash/properties.py +++ b/pyxlma/lmalib/flash/properties.py @@ -5,6 +5,32 @@ from pyxlma.lmalib.traversal import OneToManyTraversal def local_cartesian(lon, lat, alt, lonctr, latctr, altctr): + """Converts lat, lon, altitude points to x, y, z distances from a center point. + + Parameters + ---------- + lon : array_like + Longitude of points in degrees. + lat : array_like + Latitude of points in degrees. + alt : array_like + Altitude of points in meters. + lonctr : float + Longitude of the center point in degrees. + latctr : float + Latitude of the center point in degrees. + altctr : float + Altitude of the center point in meters. + + Returns + ------- + x : array_like + x distance in meters from the center point. + y : array_like + y distance in meters from the center point. + z : array_like + z distance in meters from the center point. + """ Re = 6378.137e3 #Earth's radius in m latavg, lonavg, altavg = latctr, lonctr, altctr x = Re * (np.radians(lon) - np.radians(lonavg)) * np.cos(np.radians(latavg)) @@ -13,9 +39,30 @@ def local_cartesian(lon, lat, alt, lonctr, latctr, altctr): return x,y,z def hull_volume(xyz): - """ Calculate the volume of the convex hull of 3D (X,Y,Z) LMA data. - xyz is a (N_points, 3) array of point locations in space. """ - assert xyz.shape[1] == 3 + """Calculate the volume of the convex hull of 3D (X,Y,Z) LMA data. + + Parameters + ---------- + xyz : array_like + A (N_points, 3) array of point locations in space. + + Returns + ------- + volume : float + The volume of the convex hull. + vertices : array_like + The vertices of the convex hull. + simplex_volumes : array_like + The volumes of the simplices that make up the convex hull. + + Raises + ------ + QhullError + If the convex hull cannot be computed. This is usually because too few points with little spacing are provided. See `perturb_vertex`. + + """ + if xyz.shape[1] != 3: + raise ValueError("Input must be an array of shape (N_points, 3).") tri = Delaunay(xyz[:,0:3]) vertices = tri.points[tri.simplices] @@ -37,11 +84,38 @@ def hull_volume(xyz): volume=np.sum(np.abs(simplex_volumes)) return volume, vertices, simplex_volumes + def perturb_vertex(x,y,z,machine_eps=1.0): - """ With only a few points, QHull can error on a degenerate first simplex - + """Add a random, small perturbation to an x, y, z point. + + With only a few points, QHull can error on a degenerate first simplex - all points are coplanar to machine precision. Add a random perturbation no greater than machine_eps to the first point in x, y, and z. Rerun QHull with the returned x,y,z arrays. + + Parameters + ---------- + x : array_like + x coordinates of points. + y : array_like + y coordinates of points. + z : array_like + z coordinates of points. + machine_eps : float, default=1.0 + The maximum absolute value of the perturbation. (Perturbation can be positive or negative.) + + Returns + ------- + x : array_like + x coordinates of points with perturbation added to the first point. + y : array_like + y coordinates of points with perturbation added to the first point. + z : array_like + z coordinates of points with perturbation added to the first point. + + Notes + ----- + This function is provided to work around the QHullError that can be raised by `pyxlma.lmalib.flash.properties.hull_volume` when the input points are coplanar. """ perturb = 2*machine_eps*np.random.random(size=3)-machine_eps x[0] += perturb[0] @@ -50,6 +124,26 @@ def perturb_vertex(x,y,z,machine_eps=1.0): return (x,y,z) def event_hull_area(x,y,z): + """Compute the 2D area of the convex hull of a set of x, y points. + + Parameters + ---------- + x : array_like + (N_points, 1) array of x coordinates of points. + y : array_like + (N_points, 1) array of y coordinates of points. + z : array_like + (N_points, 1) array of z coordinates of points [unused]. + + Returns + ------- + area : float + The area of the convex hull of the points. + + Notes + ----- + For more info on convex hull area and flash size, see [Bruning and MacGorman 2013](https://doi.org/10.1175/JAS-D-12-0289.1). + """ pointCount = x.shape[0] area = 0.0 if pointCount > 3: @@ -79,6 +173,22 @@ def event_hull_area(x,y,z): return area def event_hull_volume(x,y,z): + """Compute the 3D volume of the convex hull of a set of x, y points. + + Parameters + ---------- + x : array_like + (N_points, 1) array of x coordinates of points. + y : array_like + (N_points, 1) array of y coordinates of points. + z : array_like + (N_points, 1) array of z coordinates of points. + + Returns + ------- + volume : float + The volume of the convex hull of the points. + """ pointCount = x.shape[0] volume = 0.0 if pointCount > 4: @@ -96,12 +206,20 @@ def event_hull_volume(x,y,z): def rrbe(zinit): - """ - Compute the runway breakeven threshold electric fields given - an initiation altitude for a lightning flash assuming a - surface breakdown electric field threshold of 281 kVm^-1 [Marshall et al. 2005]. + """Compute the runway breakeven threshold electric field. - Returns e_init in kVm^1. + Uses a given initiation altitude for a lightning flash assuming a surface breakdown electric field + threshold of 281 kV/m, following [Marshall et al. 2005](https://doi.org/10.1029/2004GL021802). + + Parameters + ---------- + zinit : float + The altitude of the flash initiation point in meters. + + Returns + ------- + e_init : float + The electric field at the initiation point in kV/m. """ #Compute scaled air density with height. rho_air = 1.208 * np.exp(-(zinit/8.4)) @@ -109,17 +227,29 @@ def rrbe(zinit): return(e_init) def event_discharge_energy(z,area): - """ - Estimate the electrical energy discharged by lightning flashes - using a simple capacitor model. - - Note: -) Model assumes plates area defined by convex hull area, and separation - between the 73 and 27th percentiles of each flash's vertical source - distributions. - -) Only the initiation electric field (e_init) between the plates at the height - of flash initiation to find the critical charge density (sigma_crit) on the plates, - sigma_crit = epsilon * e_init -> epsilon = permitivity of air - -) Model considers the ground as a perfect conductor (image charging), + """Estimate the electrical energy discharged by lightning flashes using a simple capacitor model. + + Parameters + ---------- + z : array_like + The altitude of the flash initiation points. + area : array_like + The area of the convex hull of the flash initiation points. + + Returns + ------- + energy : array_like + The energy discharged by the flash in Joules. + + Notes + ----- + - Model assumes plates area defined by convex hull area, and separation + between the 73 and 27th percentiles of each flash's vertical source + distributions. + - Only the initiation electric field (e_init) between the plates at the height + of flash initiation to find the critical charge density (sigma_crit) on the plates, + sigma_crit = epsilon * e_init -> epsilon = permitivity of air + - Model considers the ground as a perfect conductor (image charging), """ #Compute separation between charge plates (d) and critial charge density (sigma_crit): @@ -139,7 +269,39 @@ def event_discharge_energy(z,area): def flash_stats(ds, area_func=None, volume_func=None): - + """Compute flash statistics from LMA data. + + Calculates the following variables for each flash in the dataset: + + - flash_time_start + - flash_time_end + - flash_duration + - flash_init_latitude + - flash_init_longitude + - flash_init_altitude + - flash_area + - flash_volume + - flash_energy + - flash_center_latitude + - flash_center_longitude + - flash_center_altitude + - flash_power + - flash_event_count + + Parameters + ---------- + ds : xarray.Dataset + An LMA dataset that has flash clustering applied (i.e., has a `flash_id` and `event_parent_flash_id` variable). + area_func : callable, optional + A function that computes the area of the convex hull of a set of points. If None, `event_hull_area` is used. + volume_func : callable, optional + A function that computes the volume of the convex hull of a set of points. If None, `event_hull_volume` is used. + + Returns + ------- + ds : xarray.Dataset + LMA dataset with the computed flash statistics added as variables. + """ if area_func is None: area_func = lambda df: event_hull_area(df['event_x'].array, df['event_y'].array, @@ -274,14 +436,28 @@ def flash_stats(ds, area_func=None, volume_func=None): return ds -def filter_flashes(ds, **kwargs): - """ each kwarg is a flash variable name, with tuple of minimum and maximum - values for that kwarg. min and max are inclusive (<=, >=). If either - end of the range can be None to skip it. - - Also removes events not associated with any flashes if prune=True (default) +def filter_flashes(ds, prune=True, **kwargs): + """Filter flashes by their properties. + + Allows removing unwanted flashes from an LMA dataset based on their properties. + After the flashes are removed by the criteria, the dataset can be pruned to remove any events that + are not assoicated with any flashes (or events that were associated with now-removed flashes). + + Parameters + ---------- + ds : xarray.Dataset + An LMA dataset that has flash clustering applied (i.e., has a `flash_id` and `event_parent_flash_id` variable). + prune : bool, default=True + If True, remove events not associated with any flashes. + **kwargs + Variable names and ranges to filter by. The name of the keyword argument is used as the variable name, + and ranges are given as a tuple of (min, max) values. Either end of the range can be None to skip it. + + Returns + ------- + ds : xarray.Dataset + LMA dataset with the flashes filtered by the given criteria. """ - prune = kwargs.pop('prune', True) # keep all points good = np.ones(ds.flash_id.shape, dtype=bool) # print("Starting flash count: ", good.sum()) diff --git a/pyxlma/lmalib/grid.py b/pyxlma/lmalib/grid.py index fef41ad..37f2ea1 100644 --- a/pyxlma/lmalib/grid.py +++ b/pyxlma/lmalib/grid.py @@ -3,29 +3,32 @@ import xarray as xr def discretize(x, x0, dx, int_type='uint64', bounds_check=True): - """ Calculate a unique location ID given some - discretization interval and allowed range. - - Values less than x0 raise an exception of bounds_check=True, - othewise the assigned index may wrap according when integer - casting is performed. - - Arguments: - x: coordinates, float array - x0: minimum x value - dx: discretization interval - - Keyword arguments: - int_type: numpy dtype of x_id. 64 bit unsigned int by default, since 32 - bit is limited to a 65536 pixel square grid. - - Returns: - x_id = unique pixel ID - - assert (np.array([0, 0, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3], dtype='uint64') - == discretize(np.asarray([-2.0, -1.5, -1.0, -0.1, 0.0, 0.1, 0.4, - 0.5, 0.6, 0.9, 1.0, 1.1]), -2.00, 1.0)).all() + """Calculate a unique location ID given some discretization interval and allowed range. + + Values less than x0 raise an exception of bounds_check=True, + othewise the assigned index may wrap according when integer + casting is performed. + + Parameters + ---------- + x : array_like + coordinates of the points + x0 : float + minimum x value + dx : float + discretization interval + + int_type : str, default='uint64' + numpy dtype of x_id. 64 bit unsigned int by default, since 32 bit is limited to a 65536 pixel square grid. + + Returns + ------- + x_id : array_like + unique pixel ID """ + # assert (np.array([0, 0, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3], dtype='uint64') + # == discretize(np.asarray([-2.0, -1.5, -1.0, -0.1, 0.0, 0.1, 0.4, + # 0.5, 0.6, 0.9, 1.0, 1.1]), -2.00, 1.0)).all() if bounds_check: if (x= kwargs['stn']) & \ (self.lma.chi2 <= kwargs['chi2']) & (self.lma.alt < kwargs['alt']) good2 = np.logical_and( diff --git a/pyxlma/lmalib/io/read.py b/pyxlma/lmalib/io/read.py index 4cefe2f..6ad2c50 100644 --- a/pyxlma/lmalib/io/read.py +++ b/pyxlma/lmalib/io/read.py @@ -5,6 +5,11 @@ import datetime as dt class open_gzip_or_dat: + """Helper class open a file with gzip if necessary or as binary if already decompressed. + + Use as a context manager in the same way as `with open(filename) as f:`. If the filename ends with + '.gz', the file will be decompressed with gzip. Otherwise, the file will be opened as binary. + """ def __init__(self, filename): self.filename = filename @@ -19,7 +24,22 @@ def __exit__(self, exc_type, exc_value, traceback): self.file.close() def mask_to_int(mask): - """ Convert object array of mask strings to integers""" + """Convert object array of mask strings to integers. + + Parameters + ---------- + mask : array_like + An array of strings representing the station mask. Each string should be a hexidecimal representing a binary station mask. + + Returns + ------- + mask_int : ndarray + An array of integers representing the station mask. Each integer is the binary representation of the station mask. + + Notes + ----- + Information on the station mask can be found in the [LMA Hardware section of lmaworkshop](https://github.com/deeplycloudy/lmaworkshop/blob/master/LEE-2023/LMAsourcefiles.ipynb) + """ if len(mask.shape) == 0: mask_int = np.asarray([], dtype=int) else: @@ -32,9 +52,17 @@ def mask_to_int(mask): return mask_int def combine_datasets(lma_data): - """ lma_data is a list of xarray datasets of the type returned by - pyxlma.lmalib.io.cf_netcdf.new_dataset or - pyxlma.lmalib.io.read.to_dataset + """Combine a list of LMA datasets. + + Parameters + ---------- + lma_data : iterable[xarray.Dataset] + Collection of LMA datasets of the type returned by pyxlma.lmalib.io.cf_netcdf.new_dataset or pyxlma.lmalib.io.read.to_dataset + + Returns + ------- + ds : xarray.Dataset + A combined dataset of the LMA data """ # Get a list of all the global attributes from each dataset attrs = [d.attrs for d in lma_data] @@ -94,8 +122,25 @@ def combine_datasets(lma_data): return ds def dataset(filenames, sort_time=True): - """ Create an xarray dataset of the type returned by - pyxlma.lmalib.io.cf_netcdf.new_dataset for each filename in filenames + """Read LMA .dat or .dat.gz file(s) into an xarray dataset. + + The dataset returned is the same type as is returned by + pyxlma.lmalib.io.cf_netcdf.new_dataset for each filename in filenames. + + Parameters + ---------- + filenames : str or iterable[str] + The file or files to read in + sort_time : bool, default=True + Whether to sort the VHF events by time after reading in the data. + + Returns + ------- + ds : xarray.Dataset + An xarray dataset of the LMA data. + starttime : datetime + The start of the period of record of the LMA data. + """ if type(filenames) == str: filenames = [filenames] @@ -140,7 +185,7 @@ def dataset(filenames, sort_time=True): return ds, starttime def to_dataset(lma_file, event_id_start=0): - """ lma_file: an instance of an lmafile object + """lma_file: an instance of an lmafile object returns an xarray dataset of the type returned by pyxlma.lmalib.io.cf_netcdf.new_dataset @@ -222,32 +267,33 @@ def to_dataset(lma_file, event_id_start=0): def nldn(filenames): """ - Read Viasala NLDN data + Read [Viasala NLDN](https://www.vaisala.com/en/products/national-lightning-detection-network-nldn) data Reads in one or multiple NLDN files and and returns a pandas dataframe with appropriate column names Parameters ---------- - filenames : str or list of str + filenames : str or iterable[str] The file or files to read in Returns ------- - full_df : `pandas.DataFrame` + full_df : pandas.DataFrame A pandas dataframe of entln data, with columns: - 'latitude' - the latitude of the event - 'longitude' - the longitude of the event - 'peak_current_kA' - the peak current in kA - 'multiplicity' - the number of strokes in the event - 'semimajor' - the semimajor axis length in km of the 50% confidence ellipse - 'semiminor' - the semiminor axis length in km of the 50% confidence ellipse - 'ellipseangle' - the angle of the 50% confidence ellipse - 'chi2' - the reduced chi-squared value of the event - 'num_stations' - the number of stations contributing to the event - 'type' - 'IC' or 'CG' for intracloud or cloud-to-ground - 'datetime' - the time of the event + + - 'latitude' - the latitude of the event + - 'longitude' - the longitude of the event + - 'peak_current_kA' - the peak current in kA + - 'multiplicity' - the number of strokes in the event + - 'semimajor' - the semimajor axis length in km of the 50% confidence ellipse + - 'semiminor' - the semiminor axis length in km of the 50% confidence ellipse + - 'ellipseangle' - the angle of the 50% confidence ellipse + - 'chi2' - the reduced chi-squared value of the event + - 'num_stations' - the number of stations contributing to the event + - 'type' - 'IC' or 'CG' for intracloud or cloud-to-ground + - 'datetime' - the time of the event Notes @@ -280,31 +326,32 @@ def nldn(filenames): def entln(filenames): """ - Read Earth Networks Total Lightning Network data + Read [aem/Earth Networks Total Lightning Network](https://aem.eco/product/earth-networks-total-lightning-network/) data Reads in one or multiple ENTLN files and and returns a pandas dataframe with appropriate column names Parameters ---------- - filenames : str or list of str + filenames : str or iterable[str] The file or files to read in Returns ------- - full_df : `pandas.DataFrame` + full_df : pandas.DataFrame A pandas dataframe of entln data, with columns: - 'type' - 'IC' or 'CG' for intracloud or cloud-to-ground - 'datetime' - the time of the event - 'latitude' - the latitude of the event - 'longitude' - the longitude of the event - 'peak_current_kA' - the peak current in kA - 'icheight' - the height of the IC event in meters - 'num_stations' - the number of stations contributing to the event - 'ellipseangle' - the angle of the 50% confidence ellipse - 'semimajor' - the semimajor axis length in km of the 50% confidence ellipse - 'semiminor' - the semiminor axis length in km of the 50% confidence ellipse + + - 'type' - 'IC' or 'CG' for intracloud or cloud-to-ground + - 'datetime' - the time of the event + - 'latitude' - the latitude of the event + - 'longitude' - the longitude of the event + - 'peak_current_kA' - the peak current in kA + - 'icheight' - the height of the IC event in meters + - 'num_stations' - the number of stations contributing to the event + - 'ellipseangle' - the angle of the 50% confidence ellipse + - 'semimajor' - the semimajor axis length in km of the 50% confidence ellipse + - 'semiminor' - the semiminor axis length in km of the 50% confidence ellipse Notes ----- @@ -331,20 +378,28 @@ def entln(filenames): class lmafile(object): + """Class representing an lmafile object for data being read in. To read in the data, use the `to_dataset` method.""" def __init__(self,filename): - """ - Pull the basic metadata from a '.dat.gz' LMA file - - startday : the date (datetime format) - station_info_start : the line number (int) where the station information starts - station_data_start : the line number (int) where the summarized station data starts - station_data_end : the line number (int) end of the summarized station data - maskorder : the order of stations in the station mask (str) - names : column header names - data_starts : the line number (int) where the VHF source data starts - - overview : summarized station data from file header (DataFrame, assumes fixed-width format) - stations : station information from file header (DataFrame, assumes fixed-width format) + """Pull the basic metadata from a '.dat.gz' LMA file + + Attributes + ---------- + startday : datetime + the date of the start of the period of record. + station_info_start : int + the line number in the decompressed file where the station information starts + station_data_start : int + the line number in the decompressed file where the summarized station data starts + station_data_end : int + the line number in the decompressed file of the end of the summarized station data + maskorder : str + the order of stations in the station mask + names : list + column header names + data_starts : int + the line number in the decompressed file where the VHF source data starts + stations : pd.Dataframe + station information from file header (DataFrame, assumes fixed-width format) """ self.file = filename diff --git a/pyxlma/lmalib/lma_intercept_rhi.py b/pyxlma/lmalib/lma_intercept_rhi.py index e7a810a..96b957e 100644 --- a/pyxlma/lmalib/lma_intercept_rhi.py +++ b/pyxlma/lmalib/lma_intercept_rhi.py @@ -5,30 +5,29 @@ def rcs_to_tps(radar_latitude, radar_longitude, radar_altitude, radar_azimuth): - """ - Find the unit vector coordinates (east, north, up) of the plane of a radar RHI scan. + """Find the unit vector coordinates (east, north, up) of the plane of a radar RHI scan. Creates a azimuth, elevation, range and tangent plane cartesian system at the radar's latitude and longitude, and converts the RHI azimuth direction to the tangent plane coordinate system. Parameters ---------- - radar_latitude : `float` + radar_latitude : float Latitude of the radar in degrees. - radar_longitude : `float` + radar_longitude : float Longitude of the radar in degrees. - radar_altitude : `float` + radar_altitude : float Altitude of the radar in meters. - radar_azimuth : `float` + radar_azimuth : float Azimuth of the RHI scan in degrees. Returns ---------- - X : `numpy.ndarray` + X : numpy.ndarray A 1x2 array representing the start and end points eastward component of the RHI scan. - Y : `numpy.ndarray` + Y : numpy.ndarray A 1x2 array representing the start and end points northward component of the RHI scan. - Z : `numpy.ndarray` + Z : numpy.ndarray A 1x2 array representing the start and end points upward component of the RHI scan. """ @@ -68,22 +67,22 @@ def geo_to_tps(event_longitude, event_latitude, event_altitude, tps_latitude, tp Parameters ---------- - event_longitude : `xarray.Dataset` - A pyxlma dataset containing latitude, longitude, and altitude of LMA VHF sources. - tps_latitude : `float` + event_longitude : xarray.Dataset + An LMA dataset containing latitude, longitude, and altitude of LMA VHF sources. + tps_latitude : float Latitude of the tangent plane in degrees. - tps_longitude : `float` + tps_longitude : float Longitude of the tangent plane in degrees. - tps_altitude : `float` + tps_altitude : float Altitude of the tangent plane in meters. Returns ---------- - Xlma : `numpy.ndarray` + Xlma : numpy.ndarray A 1xN array representing the eastward distance (in meters) of the tangent plane center to the LMA VHF sources. - Ylma : `numpy.ndarray` + Ylma : numpy.ndarray A 1xN array representing the northward distance (in meters) of the tangent plane center to the LMA VHF sources. - Zlma : `numpy.ndarray` + Zlma : numpy.ndarray A 1xN array representing the upward distance (in meters) of the tangent plane center to the LMA VHF sources. """ # GeographicSystem GEO - Lat, lon, alt @@ -115,20 +114,20 @@ def ortho_proj_lma(event_longitude, event_latitude, event_altitude, radar_latitu Parameters ---------- - lma_file : `xarray.Dataset` - A pyxlma dataset containing latitude, longitude, and altitude of N number of LMA VHF sources. - radar_latitude : `float` + lma_file : xarray.Dataset + An LMA dataset containing latitude, longitude, and altitude of N number of LMA VHF sources. + radar_latitude : float Latitude of the radar in degrees. - radar_longitude : `float` + radar_longitude : float Longitude of the radar in degrees. - radar_altitude : `float` + radar_altitude : float Altitude of the radar in meters. - radar_azimuth : `float` + radar_azimuth : float Azimuth of the RHI scan in degrees. Returns ---------- - lma_file_loc : `numpy.ndarray` + lma_file_loc : numpy.ndarray A Nx3 array representing the distance along, distance from, and height above the ground (in m) of the LMA VHF sources. """ @@ -203,39 +202,39 @@ def find_points_near_rhi(event_longitude, event_latitude, event_altitude, event_ Parameters ---------- - event_longitude : array-like + event_longitude : array_like An array of the latitudes of events to be transformed. - event_latitude : array-like + event_latitude : array_like An array of the latitudes of events to be transformed. - event_altitude : array-like + event_altitude : array_like An array of the altitudes of events to be transformed. - event_time : array-like + event_time : array_like An array of the times of events to be transformed. - radar_latitude : `float` + radar_latitude : float Latitude of the radar in degrees. - radar_longitude : `float` + radar_longitude : float Longitude of the radar in degrees. - radar_altitude : `float` + radar_altitude : float Altitude of the radar in meters. - radar_azimuth : `float` + radar_azimuth : float Azimuth of the RHI scan in degrees. - radar_scan_time : `datetime.datetime` or `numpy.datetime64` or `pandas.Timestamp` + radar_scan_time : datetime.datetime or numpy.datetime64 or pandas.Timestamp Time of the RHI scan. - distance_threshold : `float` + distance_threshold : float Maximum distance from the radar to the LMA VHF sources in meters. Default is 1000. - time_threshold : `float` + time_threshold : float Number of seconds before and after the RHI scan time to include LMA VHF sources. Default is 30. (total length: 1 minute) Returns ---------- - lma_range : `numpy.ndarray` + lma_range : numpy.ndarray A 1D array representing the distance along the tangent plane in the direction of the RHI scan. - lma_dist : `numpy.ndarray` + lma_dist : numpy.ndarray A 1D array representing the distance from the radar RHI scan plane to each filtered LMA point. - lma_alt : `numpy.ndarray` + lma_alt : numpy.ndarray A 1D array representing the height above the tangent plane centered at radar level of each filtered LMA point. - point_mask : `numpy.ndarray` + point_mask : numpy.ndarray A 1D array of booleans representing the VHF points that were included in the return. """ @@ -274,33 +273,33 @@ def find_lma_points_near_rhi(lma_file, radar_latitude, radar_longitude, radar_al Parameters ---------- - lma_file : `xarray.Dataset` - A pyxlma dataset containing latitude, longitude, and altitude, and event_id of N number of LMA VHF sources. - radar_latitude : `float` + lma_file : xarray.Dataset + An LMA dataset containing latitude, longitude, and altitude, and event_id of N number of LMA VHF sources. + radar_latitude : float Latitude of the radar in degrees. - radar_longitude : `float` + radar_longitude : float Longitude of the radar in degrees. - radar_altitude : `float` + radar_altitude : float Altitude of the radar in meters. - radar_azimuth : `float` + radar_azimuth : float Azimuth of the RHI scan in degrees. - radar_scan_time : `datetime.datetime` or `numpy.datetime64` or `pandas.Timestamp` + radar_scan_time : datetime.datetime or numpy.datetime64 or pandas.Timestamp Time of the RHI scan. - distance_threshold : `float` + distance_threshold : float Maximum distance from the radar to the LMA VHF sources in meters. Default is 1000. - time_threshold : `float` + time_threshold : float Number of seconds before and after the RHI scan time to include LMA VHF sources. Default is 30. (total length: 1 minute) Returns ---------- - lma_range : `numpy.ndarray` + lma_range : numpy.ndarray A 1D array representing the distance along the tangent plane in the direction of the RHI scan. - lma_dist : `numpy.ndarray` + lma_dist : numpy.ndarray A 1D array representing the distance from the radar RHI scan plane to each filtered LMA point. - lma_alt : `numpy.ndarray` + lma_alt : numpy.ndarray A 1D array representing the height above the tangent plane centered at radar level of each filtered LMA point. - point_mask : `numpy.ndarray` + point_mask : numpy.ndarray A 1D array of booleans representing the VHF points that were included in the return. """ return find_points_near_rhi(lma_file.event_longitude.data, lma_file.event_latitude.data, lma_file.event_altitude.data, diff --git a/pyxlma/lmalib/traversal.py b/pyxlma/lmalib/traversal.py index 54cdc01..9b5a551 100644 --- a/pyxlma/lmalib/traversal.py +++ b/pyxlma/lmalib/traversal.py @@ -1,32 +1,13 @@ import collections, itertools import numpy as np -""" Code from glmtools, where there are also unit tests for this class. - TODO: port to work as an xarray accessor? And move unit tests here. Adapt - to automatically use cf-tree metadata. -""" +# """ Code from glmtools, where there are also unit tests for this class. +# TODO: port to work as an xarray accessor? And move unit tests here. Adapt +# to automatically use cf-tree metadata. +# """ class OneToManyTraversal(object): - def __init__(self, dataset, entity_id_vars, parent_id_vars): - """ - dataset is an xarray.Dataset - - entity_id_vars is a list of variable names giving indices that are unique - along the dimension of that variable. The names should be given in - order from the parent to children, e.g. - ('child_id', 'grandchild_id', 'greatgrandchild_id') - which corresponds to - ('child_dim', 'grandchild_dim', 'greatgrandchild_dim') - - parent_id_vars is a list of variable names giving the child to parent - relationships, and is along the dimension of the child. The names - should be given in order from parent to children, and be one less in - length than entitiy_id_vars, e.g. - ('parent_id_of_grandchild', 'parent_id_of_greatgrandchild') - which corresponds to - ('grandchild_dim', 'greatgrandchild_dim') - - + """A class to allow traversal of a dataset where data variables have a one-to-many relationships. This object creates groupbys stored by dictionaries that make it convenient to look up the groupby given either the p @@ -53,6 +34,19 @@ def __init__(self, dataset, entity_id_vars, parent_id_vars): of virtual dimension that links set membership instead of coordinates. A future NetCDF spec could formalize this relationship structure to encourage the same library-level functionality as this class. + """ + def __init__(self, dataset, entity_id_vars, parent_id_vars): + """Initialize a OneToManyTraversal object. + + Parameters + ---------- + dataset : xarray.Dataset + The dataset to be traversed. + entity_id_vars : iterable[str] + The names of the N variables to be traversed, in order from grandparent -> parent -> child -> grandchild -> ... + Variables must be unique along the dimension of the variable. + parent_id_vars : iterable[str] + The names of the N-1 variables that link the entities in entity_id_vars, in order from (grandparent_id_of_parent) -> (parent_id_of_child) -> (child_id_of_grandchild) -> ... """ n_entities = len(entity_id_vars) @@ -107,7 +101,7 @@ def _ascend(self): yield from zip(self.entity_id_vars[::-1], self.parent_id_vars[::-1]) def count_children(self, entity_id_var, child_entity_id_var=None): - """ Count the children of entity_id_var. + """Count the children of entity_id_var. Optionally, accumulate counts of children down to and including the level of child_entity_id_var. These are the counts from parent @@ -116,14 +110,27 @@ def count_children(self, entity_id_var, child_entity_id_var=None): If replicate_parent_ids has been used to create 'bottom_parent_top_id', where the top and bottom are separated by a few levels, it is possible to get an aggregate count (matching the top parent dimension) of the - children many generations below by doing: - > grouper = dataset.groupby('bottom_parent_top_id').groups - > count = [len(grouper[eid]) if (eid in grouper) else 0 + children many generations below: + ```py + grouper = dataset.groupby('bottom_parent_top_id').groups + count = [len(grouper[eid]) if (eid in grouper) else 0 for eid in d['top_id'].data] - > assert_equal(storm_child_trig_count, count) - - Returns a list of counts of the children of entity_id_var and - (optionally) its children. + assert_equal(storm_child_trig_count, count) + ``` + + Parameters + ---------- + entity_id_var : str + The name of the variable to count children of. + child_entity_id_var : str, optional + The name of the lowest variable in the hierarchy to count children of. + If None, only the immediate children are counted. + + Returns + ------- + all_counts : tuple of numpy.ndarray + A tuple of arrays, where each array is the count of children at the + corresponding level in the hierarchy. """ count_next = False all_counts = [] @@ -142,30 +149,38 @@ def count_children(self, entity_id_var, child_entity_id_var=None): return all_counts def replicate_parent_ids(self, entity_id_var, parent_id_var): - """ Ascend from the level of parent_id_var to entity_id_var, - and replicate the IDs in entity_id_var, one each for the - number of rows in the parent_id_var dimension. - - entity_id_var is one of the variables originally - given by entity_id_vars upon class initialization. - - parent_id_var is one of the variables originally - given by parent_id_vars upon class initialization. - - This function is not strictly needed for queries up one level - (where parent_id_var can be used directly), but is needed to ascend - the hierarchy by two or more levels. - - Once the parent entity_ids have been replicated, they can be used - with xarray's indexing functions to replicate any other variable - along that parent entity's dimension. - - returns replicated_parent_ids - - TODO: args should really be: - replicate_to_dim - replicate_var -> replicate_from_dim + """Replicate the IDs of the ancestors at the level of a child entity. + + If given a mapping of child->parent, this function can find the grandparents, great-grandparents, etc. + + Parameters + ---------- + entity_id_var : str + The name of the ancestor entity to find. Must be a variable originally specificied to `entity_id_vars` the class initialization. + parent_id_var : str + The name of initial the child->parent relationship to start ascending through the heirarchy with. + Must be a variable originally specified to the `parent_id_vars` and must be lower than the value specified to `entity_id_var` + + Returns + ------- + last_replicated_p_ids : numpy.ndarray + The replicated IDs of the ancestors (at the level of `entity_id_var`) replicated to the level of the child in the child->parent + relationship specified to `parent_id_var` + + Notes + ----- + This function is not strictly needed for queries up one level + (where parent_id_var can be used directly), but is needed to ascend + the hierarchy by two or more levels. + + Once the parent entity_ids have been replicated, they can be used + with xarray's indexing functions to replicate any other variable + along that parent entity's dimension. """ + # TODO: args should really be: + # replicate_to_dim + # replicate_var -> replicate_from_dim + # Work from bottom up. # First, wait until the parent_id_var is reached # then use the groupby corresponding to the next level up @@ -201,8 +216,23 @@ def replicate_parent_ids(self, entity_id_var, parent_id_var): # last_e_var, last_p_var = e_var, p_var def reduce_to_entities(self, entity_id_var, entity_ids): - """ Reduce the dataset to the children and parents of entity_id_var - given entity_ids that are within entity_id_var + """Reduce the dataset to the ancestors and descendents of the given entity_ids that are within the given entity_id_var. + + Finds all ancestors and descendents of the specified IDs on the specified variable, and returns a filtered dataset containing only the requested + entity IDs and their ancestors and children. + + Parameters + ---------- + entity_id_var : str + The variable to find the ancestors and descendants of. + + entity_ids : array_like + The IDs to filter the dataset by. + + Returns + ------- + dataset : xarray.Dataset + original dataset filtered to the requested entity_ids along the entity_id_var and all ancestors and descendants. """ entity_ids = np.asarray(entity_ids) diff --git a/pyxlma/plot/interactive.py b/pyxlma/plot/interactive.py index c8038d0..7cc64b6 100644 --- a/pyxlma/plot/interactive.py +++ b/pyxlma/plot/interactive.py @@ -25,25 +25,32 @@ def __init__(self, a, ax): class Accumulator(object): - """ Provides for event callbacks for matplotlib drag/release events and - axis limit changes by accumulating a series of event occurrences. - Produces a single call to func after a user interacts with the plot. - Also stores the axes that got the event, and passes them to func. - Sample usage: - from pylab import figure, show - def simple(axes): - print "update ", axes - a = Accumulator(simple) - f=figure() - ax=f.add_subplot(111) - plt=ax.plot(range(10)) - f.canvas.mpl_connect('draw_event', a.draw_event) - f.canvas.mpl_connect('button_release_event', a.mouse_up_event) - f.canvas.mpl_connect('button_press_event', a.mouse_down_event) - ax.callbacks.connect('xlim_changed', a.axis_limit_changed) - ax.callbacks.connect('ylim_changed', a.axis_limit_changed) - show() - """ + """Provides for event callbacks for matplotlib drag/release events and axis limit changes by accumulating a series of event occurrences. + + Produces a single call to func after a user interacts with the plot. + Also stores the axes that got the event, and passes them to func. + + Example + ------- + ```py + from pyxlma.plot.interactive import Accumulator + from matplotlib import pyplot as plt + + def simple(axes): + print("update ", axes) + + a = Accumulator(simple) + f=plt.figure() + ax=f.add_subplot(111) + plt=ax.plot(range(10)) + f.canvas.mpl_connect('draw_event', a.draw_event) + f.canvas.mpl_connect('button_release_event', a.mouse_up_event) + f.canvas.mpl_connect('button_press_event', a.mouse_down_event) + ax.callbacks.connect('xlim_changed', a.axis_limit_changed) + ax.callbacks.connect('ylim_changed', a.axis_limit_changed) + plt.show() + ``` + """ def __init__(self, func): self.func=func @@ -52,21 +59,27 @@ def __init__(self, func): # print('did init') def reset(self): - """ Reset flags after the update function is called. - Mouse is tracked separately. - """ + """Reset flags after the update function is called. Mouse is tracked separately.""" # print('reset') self.limits_changed = 0 self.got_draw = False self.axes = None def axis_limit_changed(self, ax): + """Bindable function for matplotlib axis limit change events. + + Parameters + ---------- + ax : matplotlib.axes.Axes + The axes that had its limits changed. + """ # print('ax limits') self.limits_changed += 1 self.axes = ax self.check_status() def draw_event(self, event): + """Bindable function for matplotlib draw events.""" # print('draw event') if self.limits_changed > 0: # only care about draw events if one of the axis limits has changed @@ -74,6 +87,7 @@ def draw_event(self, event): self.check_status() def mouse_up_event(self, event): + """Bindable function for matplotlib mouse up events.""" # print('mouse up') self.mouse_up = True self.check_status() @@ -83,7 +97,7 @@ def mouse_down_event(self, event): self.mouse_up = False def both_limits_changed(self): - """ Both x and y limits changed and the mouse is up (not dragging) + """Both x and y limits changed and the mouse is up (not dragging) This condition takes care of the limits being reset outside of a dragging context, such as the view-reset (home) button on the Matplotlib standard toolbar. @@ -114,6 +128,24 @@ def check_status(self): def event_space_time_limits(ds): + """Get the limits of the event locations and times on an LMA dataset. + + Parameters + ---------- + ds : xarray.Dataset + An LMA dataset. + + Returns + ------- + xlim : tuple + The minimum and maximum longitude of the events. + ylim : tuple + The minimum and maximum latitude of the events. + zlim : tuple + The minimum and maximum altitude of the events. + tlim : tuple + The minimum and maximum time of the events. + """ xlim = ds.event_longitude.min().values, ds.event_longitude.max().values ylim = ds.event_latitude.min().values, ds.event_latitude.max().values zlim = ds.event_altitude.min().values/1000, ds.event_altitude.max().values/1000 @@ -124,14 +156,77 @@ def event_space_time_limits(ds): class InteractiveLMAPlot(object): - + """Class representing an ipywidgets interactive LMA plot. + + Visualization handled by matplotlib, interactivity handled by ipywidgets. Works with jupyter notebooks. + + Attributes + ---------- + ds : xarray.Dataset + The LMA dataset being plotted. + stationmin : int + The minimum number of stations a source must be received by to be plotted. + plot_camp : str + The colormap to use for the plot. + point_size : int + The size of the points in the plot. + clon : float + The longitude of the center of the 40km range ring. + clat : float + The latitude of the center of the 40km range ring. + widget_output : ipywidgets.widgets.widget_output.Output + Handle to the output widget for the plot. + data_artists : list + List of all artists in the plot that change when the view subset changes. + bounds : dict + Dictionary of the current plot limits for x, y, z, and t. + lma_plot : BlankPlot + The BlankPlot object representing the current plot. + inset_ax : matplotlib.axes.Axes or cartopy.mpl.geoaxes.GeoAxes + The inset axes showing the current view subset. + this_lma_lon : numpy.ndarray + The longitude of the points of the current subset. + this_lma_lat : numpy.ndarray + The latitude of the points of the current subset. + this_lma_alt : numpy.ndarray + The altitude of the points of the current subset. + this_lma_time : numpy.ndarray + The time of the points of the current subset. + this_lma_sel : numpy.ndarray + A boolean mask representing the current subset. + + + + """ def __init__(self, ds, chimax=1.0, stationmin=6, plot_cmap='plasma', point_size=5, clon=-101.5, clat=33.5, xlim=None, ylim=None, zlim=None, tlim=None): - """ lma_plot is as returned by pyxlma's BlankPlot - - xlim, ylim, zlim, and tlim are initial limits to use with the dataset; - otherwise, the limits are taken from the limits of the data in ds. + """Create an interactive LMA plot of some data. + + Parameters + ---------- + ds : xarray.Dataset + An LMA dataset. + chimax : float, default=1.0 + Sources with a chi2 value greater/worse than this value are excluded from the plot. + stationmin : int, default=6 + Sources received by fewer than this number of stations are excluded from the plot. + plot_cmap : str, default='plasma' + The colormap to use for the plot. + point_size : int, default=5 + The size of the points in the plot. + clon : float, default=-101.5 + The longitude of the center of the 40km range ring. + clat : float, default=33.5 + The latitude of the center of the 40km range ring. + xlim : tuple, optional + The initial longitude limits of the plot. + ylim : tuple, optional + The initial latitude limits of the plot. + zlim : tuple, optional + The initial altitude limits of the plot. + tlim : tuple, optional + The initial time limits of the plot. """ xlim_ds, ylim_ds, zlim_ds, tlim_ds = event_space_time_limits(ds) @@ -189,8 +284,12 @@ def _axes_spacetime_limits(self): @output.capture() def limits_changed(self, axes): - """ updates self.bounds from the limits of all axes in the plot making sure to - include changes from axes + """updates self.bounds from the limits of all axes in the plot to a changed axis + + Parameters + ---------- + axes : matplotlib.axes.Axes + The axes that had its limits changed. """ # When we start out, we get the old axis limits, naively from the xy and tz @@ -239,6 +338,7 @@ def limits_changed(self, axes): @output.capture() def make_plot(self): + """Draw the LMA plot.""" # pull out relevant plot params from object attributes ds = self.ds plot_s = self.point_size diff --git a/pyxlma/plot/xlma_base_plot.py b/pyxlma/plot/xlma_base_plot.py index c2f14b7..69e074e 100644 --- a/pyxlma/plot/xlma_base_plot.py +++ b/pyxlma/plot/xlma_base_plot.py @@ -35,7 +35,15 @@ COORD_HIST = [0.8, 0.65, 0.13, 0.1] class FractionalSecondFormatter(Formatter): + """Helper class to format fractional seconds in time labels on BlankPlot.""" def __init__(self, axis): + """Create a FractionalSecondFormatter object. + + Parameters + ---------- + axis : matplotlib.axis.Axis + The axis to format. + """ self._axis = axis def __call__(self, x, pos=None): @@ -86,18 +94,57 @@ def __call__(self, x, pos=None): class BlankPlot(object): - """ - Generate a matching plot setup with no data - - Requires: - stime = starting time - xlim, ylim = bounds of the domain - title = title string - tlim = list of start and end times - - Will include map information if bkgmap==True + """Generate LMA plot setup with no data. + + Generates a matplotlib figure in the style of the base plot from XLMA (plan view, cross sections, alt histogram, alt-timeseries axes). + + Attributes + ---------- + zlim : iterable[float] + Altitude limits [min, max] + tlim : iterable[datetime] + Time limits [start, end] + ylim : iterable[float] + Latitude limits [min, max] + xlim : iterable[float] + Longitude limits [min, max] + bkgmap : bool + Whether ax_plan is a GeoAxes or standard matplotlib axis. + ax_plan : matplotlib.axes.Axes or cartopy.mpl.geoaxes.GeoAxes + Plan view axis + ax_th : matplotlib.axes.Axes + Time-altitude axis + ax_lon : matplotlib.axes.Axes + Longitude-altitude axis + ax_lat : matplotlib.axes.Axes + Latitude-altitude axis + ax_hist : matplotlib.axes.Axes + Altitude histogram axis + fig : matplotlib.figure.Figure + The figure object """ def __init__(self, stime, bkgmap=True, **kwargs): + """Create a BlankPlot object. + + Parameters + ---------- + stime : datetime + Starting time for the plot + bkgmap : bool, default=True + Whether to include background map projection and add country/state borders. + If False, the axis will be a standard matplotlib axis. If True, the axis will be a cartopy GeoAxes with the PlateCarree projection. + If MetPy is installed and importable, the plot will include US county borders. + zlim : iterable[float] + Altitude limits [min, max] + tlim : iterable[datetime] + Time limits [start, end] + ylim : iterable[float] + Latitude limits [min, max] + xlim : iterable[float] + Longitude limits [min, max] + title : str + The title to use for the plot. + """ self.zlim = kwargs['zlim'] self.tlim = kwargs['tlim'] self.ylim = kwargs['ylim'] @@ -110,6 +157,10 @@ def __init__(self, stime, bkgmap=True, **kwargs): self.plot(**kwargs) def set_ax_plan_labels(self): + """Sets the labels on the planview axis to match the cross section axes. + + If the background map is enabled, ensure the plan view axis has the same extent as the cross section axes. + """ self.ax_plan.set_xticks(self.ax_lon.get_xticks()) self.ax_plan.set_yticks(self.ax_lat.get_yticks()) if self.bkgmap==True: @@ -117,6 +168,7 @@ def set_ax_plan_labels(self): self.ylim[0], self.ylim[1]]) def plot(self, **kwargs): + """Draw the Blank Plot""" self.fig = plt.figure(figsize=(8.5, 11)) self.ax_th = self.fig.add_axes(COORD_TH) if self.bkgmap == True: @@ -206,10 +258,21 @@ def plot(self, **kwargs): self.ax_plan.set_aspect('auto') def subplot_labels(plot): - """ - Place letters on each subplot panel. - - Returns a list of matplotlib text artists + """Place letters on each panel of a BlankPlot. + + Parameters + ---------- + plot : BlankPlot + The BlankPlot object to label. + + Returns + ------- + list of matplotlib.text.Text + Handles of the text objects for each label. + + Notes + ----- + This function is useful for creating publication-quality figures with multiple panels which require labeling. """ a = plt.text(0.05, 0.8, '(a)', fontsize='x-large', weight='bold', horizontalalignment='center', verticalalignment='center', @@ -231,12 +294,37 @@ def subplot_labels(plot): def inset_view(plot, lon_data, lat_data, xlim, ylim, xdiv, ydiv, buffer=0.5, inset_size=0.15, plot_cmap = 'magma', bkgmap = True): - """ - Overlay an inset panel of size 'inset_size' showing a plan-view histogram - of sources at xdiv, ydiv intervals and outlining a box over xlim and ylim - with buffer of 'buffer' lat/lon degrees in the image. - - Add background map features if 'bkgmap' == True + """Overlay an inset panel showing a plan-view 2D histogram of sources. + + Parameters + ---------- + plot : BlankPlot + The BlankPlot object to add the inset to. + lon_data : array_like + longitudes of sources to be added to the histogram + lat_data : array_like + latitudes of sources to be added to the histogram + xlim : iterable[float] + x (or longitude) limits [min, max] + ylim : iterable[float] + y (or latitude) limits [min, max] + xdiv : float + x (or longitude) bin width for histogram + ydiv : float + y (or latitude) bin width for histogram + buffer : float, default=0.5 + x/y buffer to be added/subtracted from xlim/ylim in the histogram bin edges + inset_size : float, default=0.15 + Size of the inset panel as a fraction of the figure size + plot_cmap : str, default='magma' + Colormap to use for the histogram + bkgmap : bool, default=True + Whether to include background map projection and add country/state borders. + + Returns + ------- + inset : cartopy.mpl.geoaxes.GeoAxes or matplotlib.axes.Axes + The axis object containing the inset histogram. """ inset = plot.fig.add_axes([0.02, 0.01, 0.02+inset_size, 0.01+inset_size],projection=ccrs.PlateCarree()) diff --git a/pyxlma/plot/xlma_plot_feature.py b/pyxlma/plot/xlma_plot_feature.py index 1f01800..21c4002 100644 --- a/pyxlma/plot/xlma_plot_feature.py +++ b/pyxlma/plot/xlma_plot_feature.py @@ -4,15 +4,52 @@ import matplotlib.dates as md from matplotlib.patches import Polygon from matplotlib.collections import PatchCollection +from pyxlma.plot.xlma_base_plot import BlankPlot def subset(lon_data, lat_data, alt_data, time_data, chi_data,station_data, xlim, ylim, zlim, tlim, xchi, stationmin): - """ - Generate a subset of x,y,z,t of sources based on maximum - reduced chi squared and given x,y,z,t bounds + """Generate a subset of x,y,z,t of sources based on maximum reduced chi squared and given x,y,z,t bounds - Returns: longitude, latitude, altitude, time and boolean arrays + Parameters + ---------- + lon_data : array_like + longitudes of the sources + lat_data : array_like + latitudes of the sources + alt_data : array_like + altitudes of the sources + time_data : array_like + times of the sources + chi_data : array_like + reduced chi squared values of the sources + station_data : array_like + the number of stations receiving each source + xlim : iterable + [min, max] longitude bounds + ylim : iterable + [min, max] latitude bounds + zlim : iterable + [min, max] altitude bounds + tlim : iterable + [min, max] time bounds + xchi : float + maximum reduced chi squared value to be allowed + stationmin : int + minimum number of stations to be allowed + + Returns + ------- + lon_data : numpy.ndarray + longitudes of the sources + lat_data : numpy.ndarray + latitudes of the sources + alt_data : numpy.ndarray + altitudes of the sources + time_data : numpy.ndarray + times of the sources + selection : numpy.ndarray + boolean array of the sources that meet the criteria """ selection = ((alt_data>zlim[0])&(alt_dataxlim[0])&(lon_data