diff --git a/.gitmodules b/.gitmodules index 07c36556..318e0694 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,3 +1,4 @@ [submodule "hdw"] path = pydarn/utils/hdw url = https://github.com/superdarn/hdw + branch = main \ No newline at end of file diff --git a/.zenodo.json b/.zenodo.json index 483724ea..daf2c515 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -1,21 +1,26 @@ { "creators":[ { - "name": "SuperDARN Data Analysis Working Group" + "name": "SuperDARN Data Visualization Working Group" }, { "affiliation": "University of Saskatchewan", "name": "Schmidt, M.T.", "orcid": "0000-0002-3265-977X" }, - { "affiliation": "University of Scranton", - "name": "Tholley, F." - }, { "affiliation": "University of Saskatchewan", "name": "Martin, C.J.", "orcid": "0000-0002-8278-9783" }, + { + "affiliation": "Virginia Tech", + "name": "Shi, X.", + "orcid": "0000-0001-8425-8241" + }, + { "affiliation": "University of Scranton", + "name": "Tholley, F." + }, { "affiliation": "University of Saskatchewan", "name": "Billett, D.D.", @@ -31,32 +36,24 @@ "name": "Coyle, S.", "orcid": "0000-0003-1730-2753" }, - { - "affiliation": "University of Saskatchewan", - "name": "Detwiller, M.H." - }, { "affiliation": "University of Scranton", "name": "Frissell, N.", "orcid": "0000-0002-8398-4222" }, - { - "affiliation": "The University Centre in Svalbard", - "name": "Herlingshaw, K.", - "orcid": "0000-0001-6861-1914" - }, { "affiliation": "University of Saskatchewan", "name": "Huyghebaert, D.", "orcid": "0000-0002-4257-4235" }, + { + "affiliation": "University of Alabama", + "name": "Khanal, K.", + "orcid": "" + }, { "affiliation": "University of Bergen", "name": "Reistad, J.P.", "orcid": "0000-0003-3509-5479" - }, - { - "affiliation": "University of Saskatchewan", - "name": "Roberston, C.R." }] } diff --git a/README.md b/README.md index ad6e9713..e62bde68 100644 --- a/README.md +++ b/README.md @@ -9,18 +9,13 @@ Python data visualization library for the Super Dual Auroral Radar Network (Supe ## Changelog -## Version 2.2.1 - Release! +## Version 3.0 - Release! -**New Requirement**: pyDARN 2.2.1 requires minimum matplotlib version of 3.3.4 +**New Requirement**: pyDARN 3.0 requires minimum matplotlib version of 3.3.4 -pyDARN release v2.2.1 includes the following features: -* Bug fix: fixed but in Summary plot not rendering -* Bug fix: fixed axes object not being properly handled in plot_fan -* Enhancement: Added `__version__` ability in pydarn -* Updated: citing warning printing on `import` instead of plotting -* Updated: pyDARNio 1.1.0 is requirement - -**Deprecation**: `slant` option in `plot_range-time` and `plot_summary` is deprecated now uses `coords` +pyDARN release v3.0 includes the following features: +- **New** optional cartopy dependency +- **New** convection map plotting ## Documentation diff --git a/docs/imgs/fan_4.png b/docs/imgs/fan_4.png new file mode 100644 index 00000000..078f2109 Binary files /dev/null and b/docs/imgs/fan_4.png differ diff --git a/docs/imgs/fov_5.png b/docs/imgs/fov_5.png new file mode 100644 index 00000000..42781809 Binary files /dev/null and b/docs/imgs/fov_5.png differ diff --git a/docs/imgs/fov_6.png b/docs/imgs/fov_6.png new file mode 100644 index 00000000..4db0de4b Binary files /dev/null and b/docs/imgs/fov_6.png differ diff --git a/docs/imgs/fov_7.png b/docs/imgs/fov_7.png new file mode 100644 index 00000000..e08edc2f Binary files /dev/null and b/docs/imgs/fov_7.png differ diff --git a/docs/imgs/map_1.png b/docs/imgs/map_1.png new file mode 100644 index 00000000..8cff69f4 Binary files /dev/null and b/docs/imgs/map_1.png differ diff --git a/docs/imgs/map_2.png b/docs/imgs/map_2.png new file mode 100644 index 00000000..82251220 Binary files /dev/null and b/docs/imgs/map_2.png differ diff --git a/docs/index.md b/docs/index.md index b454fe39..296552ea 100644 --- a/docs/index.md +++ b/docs/index.md @@ -25,6 +25,7 @@ If you have any questions or concerns please submit an **Issue** on the SuperDAR - [FOV plots](user/fov.md) - [Fan plots](user/fan.md) - [Grid plots](user/grid.md) + - [Convection Map Potential plots](user/map.md) - [Power plots](user/power.md) - [ACF plotting](user/acf.md) - [Logging](user/logging.md) diff --git a/docs/user/coordinates.md b/docs/user/coordinates.md index e543f38e..737075e1 100644 --- a/docs/user/coordinates.md +++ b/docs/user/coordinates.md @@ -15,13 +15,15 @@ the additional permissions listed below. # Coordinate Systems --- -pyDARN offers several different coordinate systems and units for plotting: +pyDARN offers several different measurement systems for various types of plotting: +- Range Estimation: the estimate of how far the target (echo) is from the radar +- Coordinate systems: Generic geographic coordinate system with conversions to AACGM and AACGM MLT ## Range-Time Plots -**Range Gates**: `Coords.RANGE_GATE` a rectangle section determined by beam width and set distance for each range (nominally 45 km). RAWACF and FITACF data give their parameter values with respect to range gates. Range gates are a unit-less measure of estimating distance. +**Range Gates**: `RangeEstimation.RANGE_GATE` a rectangle section determined by beam width and set distance for each range (nominally 45 km). RAWACF and FITACF data give their parameter values with respect to range gates. Range gates are a unit-less measure of estimating distance. -**Slant Range**: `Coords.SLANT_RANGE` is a conversion from range gates to km units. Slant range estimates the distance of ionospheric echos from the radar, using the time it takes for the radio wave to travel to the ionosphere and return, assuming the radio wave is travelling at the speed of light. +**Slant Range**: `RangeEstimation.SLANT_RANGE` is a conversion from range gates to km units. Slant range estimates the distance of ionospheric echos from the radar, using the time it takes for the radio wave to travel to the ionosphere and return, assuming the radio wave is travelling at the speed of light. !!! Note Slant range is calculated from the value of `frang`, the distance to the first range gate. In pyDARN, we assume @@ -32,7 +34,7 @@ pyDARN offers several different coordinate systems and units for plotting: the value given in hardware files. Currently, pyDARN has decided to use the values for `rxrise` given in the hardware files. We will amend or reconsider this approach as and when a solution to the differing values is found. -**Ground Scatter Mapped Range**: `Coords.GROUND_SCATTER_MAPPED_RANGE` uses echos from ground scatter to adjust slant-range coordinates to be more accurate based on [Dr. Bill Bristow's paper](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/93JA01470). Implemented by Dr. Nathaniel Frissell and Francis Tholley from University of Scranton. Measured in km. +**Ground Scatter Mapped Range**: `RangeEstimation.GSMR` uses echos from ground scatter to adjust slant-range coordinates to be more accurate based on [Dr. Bill Bristow's paper](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/93JA01470). Implemented by Dr. Nathaniel Frissell and Francis Tholley from University of Scranton. Measured in km. **Geographic Latitude** (Coming Soon) @@ -46,6 +48,10 @@ Geographic plots include: fan, grid and convection map (coming soon) plots **AACGM**: `Coords.AACGM` is [Altitude Adjusted Corrected Geogmagnetic Coordinates developed by Dr. Simon Shepherd](http://superdarn.thayer.dartmouth.edu/aacgm.html) are an extension of corrected geomagnetic (CGM) coordinates that more accurately represent the actual magnetic field. In AACGM coordinates points along a given magnetic field line are given the same coordinates and thus are a better reflection of magnetic conjugacy. pyDARN uses AACGM-V2 from the [aacgmv2 python library](https://pypi.org/project/aacgmv2/). Implemented by Dr. Daniel Billett from University of Saskatchewan. -**Ground Scatter Mapped Geographic** (Coming Soon) +**AACGM_MLT**: `Coords.AACGM_MLT` is `Coords.AACGM` with the geomagnetic longitude converted to magnetic local time + +`RangeEstimation` methods can be used with a `Coords` calculation. For example, using `Coords.GEOGRAPHIC` and `RangeEstimation.GSMR` together, will give a plot of ionospheric echoes at a distance from the radar calculated in ground scatter mapped range, in geographic coordinates. + +!!! Warning + You cannot use `RangeEstimation.RANGE_GATES` with any `Coords`, the default is `RangeEstimation.SLANT_RANGE` -**Ground Scatter Mapped Geomagnetic** (Coming Soon) diff --git a/docs/user/fan.md b/docs/user/fan.md index 4b8eddf3..03fe671f 100644 --- a/docs/user/fan.md +++ b/docs/user/fan.md @@ -18,10 +18,13 @@ the additional permissions listed below. Fan plots are a way to visualise data from the entire scan of a SuperDARN radar. -All beams and ranges for a given parameter (such as line-of-sight velocity, backscatter power, etc) and a particular scan are projected onto a polar format plot in [AACGMv2](http://superdarn.thayer.dartmouth.edu/aacgm.html) coordinates. +All beams and ranges for a given parameter (such as line-of-sight velocity, backscatter power, etc) and a particular scan can be projected onto a polar format plot in [AACGMv2](http://superdarn.thayer.dartmouth.edu/aacgm.html) coordinates, or projected onto a geographic plot in geographic coordinates. -Currently, fan plots in pyDARN get the geographic positions of a radar's range gates by reading in pre-generated files (found in the `/pydarn/radar_fov_files` folder), and then [converts](https://pypi.org/project/aacgmv2/) them to AACGMv2 coordinates. In the future, pyDARN will generate the geographical position, which will bring support for not standard range/beam layouts. -The mapping of the range gate corners was based on [rbpos in RST](https://github.com/SuperDARN/rst/blob/0aa1fffed4cc48c1eb6372dfc9effa688af95624/codebase/superdarn/src.idl/lib/legacy.1.6/rbposlib.pro). +!!! Warning + At present, AACGM coordinates cannot be plotted on a geographic projection. + + +The mapping of the range gate corners was based on [rbpos in RST](https://github.com/SuperDARN/rst/blob/0aa1fffed4cc48c1eb6372dfc9effa688af95624/codebase/superdarn/src.idl/lib/legacy.1.6/rbposlib.pro) that is coded in `/pydarn/geo.py`. To get the coordinates read [hardware docs](hardware.md). ### Basic usage pyDARN and pyplot need to be imported, as well as any FITACF file needs to be [read in](https://pydarn.readthedocs.io/en/latest/user/SDarnRead/): @@ -56,9 +59,12 @@ pydarn.Fan.plot_fan(cly_data, scan_index=datetime(2015, 3, 8, 15, 26), !!! Warning Do not include seconds, typically scans are 1 minute long so seconds may end in a error with no data. -Default plots also do not show groundscatter as grey. Set it to true to colour groundscatter this way: +Default plots also do not show groundscatter as grey. Set it to true to colour groundscatter: + ```python -fanplot = pydarn.Fan.plot_fan(fitacf_data, scan_index=27, groundscatter=1) +fanplot = pydarn.Fan.plot_fan(fitacf_data, + scan_index=27, + groundscatter=True) plt.show() ``` @@ -66,7 +72,8 @@ plt.show() You might have noticed that the variable `fanplot` in the examples above actually holds some information. This contains the AACGM latitude and longitude of the fan just plotted, as well as the data, ground scatter information, and datetime object. If you instead change `fanplot` to 5 separate variables, it will return the latitude, longitude, data, groundscatter and datetime info into seperate variables: ```python -lats,lons,data,grndsct,datetime=pydarn.Fan.plot_fan(fitacf_data, scan_index=27) +ax, lats, lons, data, grndsct, datetime = pydarn.Fan.plot_fan(fitacf_data, + scan_index=27) lats.shape @@ -90,21 +97,21 @@ Here is a list of all the current options than can be used with `plot_fan` | Option | Action | | ----------------------------- | ------------------------------------------------------------------------------------------------------- | +| ax=(Axes Object) | Matplotlib axes object than can be used for cartopy additions | | scan_index=(int or datetime) | Scan number or datetime, from start of records in file corresponding to channel if given | | channel=(int or 'all') | Specify channel number or choose 'all' (default = 'all') | +| parameter=(string) | See above table for options | | groundscatter=(bool) | True or false to showing ground scatter as grey | | ranges=(list) | Two element list giving the lower and upper ranges to plot, grabs ranges from hardware file (default [] | -| cmap=matplotlib.cm | A matplotlib color map object. Will override the pyDARN defaults for chosen parameter | +| cmap=string | A matplotlib color map string | | zmin=(int) | Minimum data value for colouring | | zmax=(int) | Maximum data value for colouring | | colorbar=(bool) | Set true to plot a colorbar (default: True) | | colorbar_label=(string) | Label for the colour bar (requires colorbar to be true) | +| title=(string) | Title for the fan plot, default auto generated one based on input information | | boundary=(bool) | Set false to not show the outline of the radar FOV (default: True) | -| fov_color=(string) | Sets the fill in color for the fov plot (default: transparency) | -| line_color=(string) | Sets the boundary line and radar location dot color (default: black) | -| alpha=(int) | Sets the transparency of the fill color (default: 0.5) | -| radar_location=(bool) | Places a dot in the plot representing the radar location (default: True) | -| radar_label=(bool) | Places the radar 3-letter abbreviation next to the radar location | +| coords=(Coords) | [Coordinates](coordinates.md) for the data to be plotted in | +| projs=(Projs) | Projections to plot the data on top of | | kwargs ** | Axis Polar settings. See [polar axis](axis.md) | @@ -139,3 +146,16 @@ plt.show() ``` ![](../imgs/fan_3.png) + +Using *cartopy* to plot underlaid coastline map: + +!!! Warning + The *cartopy* coastlines mapping feature is currently only available for geographic projections in geographic coordinates. Expansion of this feature is coming soon! + +```python +ax, _, _, _, _ = pydarn.Fan.plot_fan(data, scan_index=5, radar_label=True, groundscatter=True, coords=pydarn.Coords.GEOGRAPHIC, projs=pydarn. Projs.GEO, colorbar_label="Velocity m/s") +ax.coastlines() +plt.show() +``` + +![](../imgs/fan_4.png) diff --git a/docs/user/fov.md b/docs/user/fov.md index 6b599a31..d27bab13 100644 --- a/docs/user/fov.md +++ b/docs/user/fov.md @@ -1,3 +1,18 @@ + + # Field Of View Plots --- @@ -14,10 +29,7 @@ plt.show() ![](../imgs/fov_1.png) -A `datetime` object of the date is required to convert to AACGM MLT coordinates (see [coordinates](coordinates.md)). - -!!! Note - This will not be required for other coordinates once implemented. +A `datetime` object of the date is required to convert to `Coords.AACGM_MLT` (default) or `Coords.AACGM`. ### Additional options @@ -26,9 +38,15 @@ Here is a list of all the current options than can be used with `plot_fov` | Option | Action | | ----------------------- | ------------------------------------------------------------------------------------------------------- | | stid=(int) | Station id of the radar. Can be found using [SuperDARNRadars](hardware.md) | -| date=(datetime) | `datetime` object to determine the position the radar fov is in MLT +| date=(datetime) | `datetime` object to determine the position the radar fov AACGM or AACGM MLT coordinates | | ranges=(list) | Two element list giving the lower and upper ranges to plot, grabs ranges from hardware file (default [] | -| colorbar=(bool) | Set true to plot a colorbar (default: True) | +| ccrs=(object) | Cartopy axes object for plotting using Cartopy | +| rsep=(int) | Range Seperation (km) (default: 45 km) | +| frang=(int) | Frequency Range (km) (default: 180 km) | +| projs=(Projs) | Projection for the plot to be plotted on Polar and Geographic (GEO) (default: Projs.POLAR) | +| coords=(Coords) | Coorindates Geographic, AACGM, or AACGM MLT (default: Coords.AACGM_MLT) | +| grid=(bool) | Boolean to apply the grid lay of the FOV (default: False ) | +| colorbar=(bool) | Set true to plot a colorbar (default: True) | | colorbar_label=(string) | Label for the colour bar (requires colorbar to be true) | | boundary=(bool) | Set false to not show the outline of the radar FOV (default: True) | | grid=(bool) | Set true to show the outline of the range gates inside the FOV (default: False) | @@ -74,6 +92,25 @@ plt.show() ![](../imgs/fov_3.png) + +This example shows the use of *cartopy* and plotting in geographic coordinates. + +```python +_ , _, ax, ccrs = pydarn.Fan.plot_fov(stid=65, + date=dt.datetime(2022, 1, 8, 14, 5), + fov_color='red', + coords=pydarn.Coords.GEOGRAPHIC, + projs=pydarn.Projs.GEO) +ax.coastlines() +plt.show() +``` + +![](../imgs/fov_7.png) + +!!! Warning + Currently you cannot plot AACGM coordinates on a geographic plot as its not correctly transformed. Currently in development. + + `plot_fov` use two other plotting methods `plot_radar_position` and `plot_radar_label`, these methods have the following parameters: | Option | Action | @@ -95,4 +132,6 @@ pydarn.Fan.plot_fov(66, datetime(2021, 2, 5, 12, 5), boundary=False, ![](../imgs/fov_4.png) !!! Note - The radar labels will appear at the same longitude, but at -5 degrees of latitude to the position of the radar station. This may cause some to overlap. Users can plot their own labels using `plt.text(*lon psn in radians*, *lat psn in degrees*, *text string*)` + The radar labels will appear at the same longitude, but at -5 degrees of latitude to the position of the radar station. This may cause some to overlap. Users can plot their own labels using `plt.text(*lon psn in radians*, *lat psn in degrees*, *text string*)` + + diff --git a/docs/user/hardware.md b/docs/user/hardware.md index d51da1a6..b379d454 100644 --- a/docs/user/hardware.md +++ b/docs/user/hardware.md @@ -52,12 +52,13 @@ Other information a user can access from the `_HdwInfo` object is: | :---: | :--- | | `stid` | Station Id of the radar | | `abbrev` | 3 letter radar abbreviation | +| `date` | Date the hardware specifications were changed | | `geographic` | Geographic coordinates of the radar and altitude in meters (lat, long, alt) | -| `boresight` | Boresight of the centre beam | +| `boresight` | Boresight of the centre beam (physical, electronic) | | `beam_separation` | Angular separation between radar beams in degrees | | `velocity_sign` | To help identify backscatter velocities which the signs can be reversed based on receiver design | | `rx_attenuator` | Analog Rx attenuator step in dB | -| `tdiff` | propagation time from interferometer array antenna to phasing matrix input minus propagation time from main array antenna through transmitter in phasing matrix in microseconds | +| `tdiff` | propagation time from interferometer array antenna to phasing matrix input minus propagation time from main array antenna through transmitter in phasing matrix in microseconds (channel_a, channel_b)| | `phase_sign` | Account for cable error in analyzing data | | `interferometer_offset` | Cartesian coordinates (x,y,z) from midpoint interferometer array to midpoint main array in meters | | `rx_rise_time` | Analog Rx rise time measured in microseconds | @@ -68,6 +69,10 @@ Other information a user can access from the `_HdwInfo` object is: !!! Note For more detailed information on all the fields in a hardware file, please read the [hardware repo README.md](https://github.com/SuperDARN/hdw). +!!! Warning + Prior to version 3.0, pyDARN was built to use the old format of hardware files. However, versions 2.2.1 or lower of pyDARN will try to pull hardware files from the `master` branch of the hardware repository and this may cause some errors in use. + Version 3.0 uses the new format of hardware files, and pulls hardware files from the `main` hardware branch. Updating to pyDARN version 3.0 or higher will fix any hardware errors. + # Accessing Radar and Hardware Information Another way to access the hardware information, the radar's full name, the institution's name and the hemisphere that the radar is located in is by using the `SuperDARNRadars` class with the station id number (`stid` field in most files). @@ -82,7 +87,7 @@ print(radar_info) Expected output: ```python -_Radar(name='Prince George', institution='University of Saskatchewan', hemisphere=, hardware_info=_HdwInfo(stid=6, abbrev='pgr', geographic=_Coord(lat=53.98, lon=-122.59, alt=670.0), boresight=-5.0, beam_separation=3.24, velocity_sign=1.0, rx_attenuator=10.0, tdiff=0.0, phase_sign=1.0, interferometer_offset=_InterferometerOffset(x=0.0, y=-100.0, z=0.0), rx_rise_time=0.0, attenuation_stages=0.0, gates=225, beams=16)) +_Radar(name='Prince George', institution='University of Saskatchewan', hemisphere=, hardware_info=_HdwInfo(stid=6, status=, abbrev='pgr', date=datetime.datetime(2000, 3, 3, 0, 0), geographic=_Coord(lat=53.98, lon=-122.59, alt=670.0), boresight=_Boresight(physical=-5.0, electronic=0.0), beam_separation=3.24, velocity_sign=1.0, rx_attenuator=10.0, tdiff=_Tdiff(channel_a=0.0, channel_b=0.0), phase_sign=1.0, interferometer_offset=_InterferometerOffset(x=0.0, y=-100.0, z=0.0), rx_rise_time=0.0, attenuation_stages=0, gates=225, beams=16)) ``` !!! Warning @@ -97,7 +102,7 @@ Example code: import pydarn # Geographic coordinates for Clyde River (STID: 66) FOV -geo_lats, geo_lons=pydarn.radar_fov(66, coords='geo') +geo_lats, geo_lons=pydarn.Coords.GEOGRAPHIC(66) ``` You also have the option to set the `coords` keyword to `aacgm`. In this case, [Altitude adjusted corrected geomagnetic](http://superdarn.thayer.dartmouth.edu/aacgm.html) latitude and longitude are returned instead of geographic. Because AACGM requires a date to convert coordinates accurately, a python datetime object is also required to be passed in to `radar_fov` under this circumstance: @@ -106,11 +111,9 @@ import pydarn import datetime as dt # AACGMv2 coordinates for Dome C (STID: 96), valid for November 26th, 2005 -aacgm_lats, aacgm_lons=pydarn.radar_fov(96, coords='aacgm', date=dt.datetime(2005,11,26)) +aacgm_lats, aacgm_lons=pydarn.Coords.AACGM(96, date=dt.datetime(2005, 11, 26)) ``` -If the `coords` keyword is not given, `radar_fov` defaults to geographic coordinates. - -The outputs for `radar_fov` are two numpy arrays of latitude and longitude coordinates with dimensions (number_of_beams+1 x number_of_gates+1). They correspond to the corners of each range gate. +The `Coords` keyword points to the function to convert the radar's Field-of-View to the designed coordinate system. The outputs are two numpy arrays of latitude and longitude coordinates with dimensions (number_of_beams+1 x number_of_gates+1). They correspond to the corners of each range gate. # Updating Radar and Hardware Information diff --git a/docs/user/install.md b/docs/user/install.md index f68447f9..530c48f7 100644 --- a/docs/user/install.md +++ b/docs/user/install.md @@ -56,6 +56,23 @@ pyDARN's setup will download the following dependencies: - [pyDARNio](https://pydarnio.readthedocs.io/en/latest/user/install/) - [AACGMv2](https://pypi.org/project/aacgmv2/) +!!! Note + If you wish to plot coastlines you will need to install cartopy>=0.19 separately + +### Cartopy +[Cartopy](https://scitools.org.uk/cartopy/docs/latest/) is a Python package designed for geospatial data processing in order to produce maps and other geospatial data analyses. This library is used when invoking a projection system needing overlapped coastline maps in Fan, Grid and Map plots. + +For installing cartopy please follow the packages [installation](https://scitools.org.uk/cartopy/docs/latest/installing.html) instructions. For ubuntu here is good installation [link](https://techoverflow.net/2021/07/11/how-to-install-cartopy-on-ubuntu/) + +!!! Warning + For cartopy to work with pyDARN please make sure it version `>=0.19`. Otherwise pydarn will throw an exception if you try to use it. + + +!!! Note + cartopy is a challenging package to install so please provide any information on troubleshoot or solutions to common issue on [pyDARN github](https://github.com/SuperDARN/pydarn) page. + + + ## Virtual Environments It is recommended to install pyDARN in one of the suggested virtual environments if you have multiple python/pip 3 version on your computer, or do not want to affect the main system's python libraries. diff --git a/docs/user/map.md b/docs/user/map.md new file mode 100644 index 00000000..c1914629 --- /dev/null +++ b/docs/user/map.md @@ -0,0 +1,121 @@ + + +# Map plots +--- + +Map plots are a way to visualize data from a MAP file of SuperDARN radar data. Please read RST documentation on how to process [MAP files](https://radar-software-toolkit-rst.readthedocs.io/en/latest/user_guide/map_grid/) from GRID files. + +Map field descriptions can be found [here](https://radar-software-toolkit-rst.readthedocs.io/en/latest/references/general/map/). pyDARN uses a `enum` object to select different common parameters to plot in a MAP file: + +| Name | parameter name | Map field name | +| ------------------ | ----------------------------- | ----------------------------- | +| Fitted Velocity | `MapParams.FITTED_VELOCITIES` | see Fitted Velocities section | +| Modeled Velocities | `MapParams.MODEL_VELOCITIES` | `model.vel.median` | +| Raw Velocities | `MapParams.RAW_VELOCITIES` | `vector.vel.median` | +| Power | `MapParams.POWER` | `vector.pwr.median` | +| Spectral Width | `MapParams.SPECTRAL_WIDTH` | `vector.wdt.median` | + +for a given `start_time` or `record` number projected onto a polar plot in [AACGMv2](http://superdarn.thayer.dartmouth.edu/aacgm.html) coordinates. + +Currently, map plots in pyDARN get geomagnetic positions of the mapped data in [`mlon` and `mlat`](https://pypi.org/project/aacgmv2/) from the MAP file, which uses AACGMv2 coordinates. In the future, pyDARN will also generate the geographic position of the data points, which will bring support for non-standard gridded vector layouts. + +### Fitted Velocities + +Fitted velocities are velocity vectors which represent the fitted convection pattern. They are entirely in the direction of the fitted convection flow and so can be fairly different to line-of-sight velocities, but still ultimately constrained by them. + +Fitted velocity vectors are by default only calculated at the same positions of the line-of-sight vectors, but fit vectors at an arbitrary position can be obtained by using the `calculated_fitted_velocities` function in `map.py`. + +## Basic usage + +pyDARN and pyplot need to be imported and the desired MAP file needs to be [read in](https://pydarn.readthedocs.io/en/master/user/SDarnRead/): + +```python +import matplotlib.pyplot as plt +import pydarn + +#Read in Map file using SuperDARNRead, then read_map +file = "path/to/grid/file" +SDarn_read = pydarn.SuperDARNRead(file) +map_data = SDarn_read.read_map() + +``` +With the map data loaded as a list of dictionaries (`map_data` variable in above example), you may now call the `plot_mapdata` method. Make sure you tell the method what time, in [`datetime` format], or record number (numbered from first recorded in file, counting from 0): +```python +mapplot = pydarn.Map.plot_mapdata(map_data, record=150) +plt.show() + +``` +In this example, the record at 150 was plotted with the defaulted parameter, `MapParams.FITTED_VELOCITIES` (fitted velocities), being mapped: +![](../imgs/map_1.png) + +You might have noticed that the variable `mapplot` in the examples above actually holds some information. This contains the AACGM latitude and longitude of the mapped vectors plotted. If you instead change `mapplot` to 3 separate variables, it will return the latitude, longitude, and data info into separate variables: +```python +lats,lons,data=pydarn.Map.plot_mapdata(map_data, start_time=stime) +``` + +### Additional options + +Here is a list of all the current options than can be used with `plot_mapdata` + +| Option | Action | +| ------------------------------ | ------------------------------------------------------------------------------------- | +| record=(int) | Record number to plot | +| start_time=(datetime.datetime) | The start time of the record to plot | +| time_delta=(int) | How close to the start time to be to the start time of the record | +| parameter=(Enum) | Specify the required parameter described above i.e. `pydarn.MapParams.FITTED_VELOCITY` | +| lowlat=(int) | Control the lower latitude boundary of the plot (default 30/-30 AACGM lat) | +| len_factor=(int) | Normalisation factor for the length of the vectors on the plot | +| boundary=(bool) | Boolean to show the Field-of-View of the radar(s) | +| cmap=matplotlib.cm | A matplotlib color map object. Will override the pyDARN defaults for chosen parameter | +| zmin=(int) | Minimum data value for colouring | +| zmax=(int) | Maximum data value for colouring | +| colorbar=(bool) | Set true to plot a colorbar (default: True) | +| colorbar_label=(string) | Label for the colour bar (requires colorbar to be true) | +| title=(str) | To add a title to the plot | +| hmb=(bool) | Set to True to include the Heppnar-Maynard Boundary on the plot | +| imf_dial=(bool) | Show the IMF data as a clock angle dial (default True) | +| map_info=(bool) | Show selected map file data (e.g. cross polar cap potential) (default True) | +| contour_levels=(list: floats) | Set custom levels for contours to appear at | +| contour_color=(str) | Set custom color for contour lines (default 'dim grey') | +| contour_linewidths=(float) | Set custom line width for contour lines (default 0.8) | +| contour_fill=(bool) | Choose to use filled contours, default is False to plot contour lines only | +| contour_colorbar=(bool) | If contour_fill=True, color bar will be displayed (default True) | +| contour_fill_cmap=(str) | If contour_fill=True, color map can be selected (default 'RdBu') | +| contour_colorbar_label=(str) | If contour_fill and contour_colorbar= True, set custom contour colorbar label | +| pot_minmax_color=(str) | Choose color of minimum and maximum potential markers | + +More `**kwargs` can be used to customise the display of the radars field-of-view if `boundary=True` + +The following is an example of the customization available: +```python +import matplotlib.pyplot as plt +import pydarn + +map_file = "20150310.n.map" +map_data = pydarn.SuperDARNRead().read_dmap(map_file) + +pydarn.Maps.plot_mapdata(map_data, record=150, + parameter=pydarn.MapParams.FITTED_VELOCITY, + lowlat=60, colorbar_label='Velocity m/s', + contour_fill = True, + contour_fill_cmap= 'RdBu', + contour_colorbar = True, + contour_colorbar_label='Potential (kV)', + pot_minmax_color = 'r', + map_info=True, imf_dial=True, hmb=True) + +plt.show() +``` +![](../imgs/map_2.png) diff --git a/docs/user/range_time.md b/docs/user/range_time.md index 617101f5..9ed484b2 100644 --- a/docs/user/range_time.md +++ b/docs/user/range_time.md @@ -56,7 +56,7 @@ To specify which beam to look at, add the option: As an example, taking a look at some `v` data from the first record of Clyde River radar FITACF file: ```python -pydarn.RTP.plot_range_time(fitacf_data, beam_num=fitacf_data[0]['bmnum'], coords=pydarn.Coords.RANGE_GATE) +pydarn.RTP.plot_range_time(fitacf_data, beam_num=fitacf_data[0]['bmnum'], range_estimation=pydarn.RangeEstimation.RANGE_GATE) plt.title("Radar {:d}, Beam {:d}".format(fitacf_data[0]['stid'], fitacf_data[0]['bmnum'])) plt.show() @@ -69,12 +69,12 @@ which produces: Notice that the velocity scale on the right is a bit larger than we need, and also ground scatter isn't coloured grey by default. Showing the dates on the x axis is also a bit redundant, because it's data from a single day. Below, there are some additional parameters you can set to address these and more. -In addition, we use `coords=pydarn.Coords.RANGE_GATE` to set the y-axis to plot in range gates. +In addition, we use `range_estimation=pydarn.RangeEstimation.RANGE_GATE` to set the y-axis to plot in range gates. Ground-Scatter Mapped Range is another type of axis you can use with range-time plots: ```python pydarn.RTP.plot_range_time(fitacf_data, beam_num=3, parameter='p_l', - coords=pydarn.Coords.GROUND_SCATTER_MAPPED_RANGE, + range_estimation=pydarn.RangeEstimation.GSMR, colorbar_label='SNR') plt.title("Clyde 20150308 14:00 UTC Fitacf 2.5") plt.ylabel('Ground Scatter Mapped (km)') @@ -97,14 +97,14 @@ To see all the customisation options, check out all the parameters listed in 'rt | date_fmt=(string) | How the x-tick labels look. Default is ('%y/%m/%d\n %H:%M') | | zmin=(int) | Minimum data value to be plotted | | zmax=(int) | Maximum data value to be plotted | -| coords=(Coords) | Coordinates to use for the y-axis (See [Coordinates](coordinates.md)) | +| range_estimation=(RangeEstimation) | Coordinates to use for the y-axis (See [Coordinates](coordinates.md)) | For instance, code for a velocity RTP showing the same beam of Clyde river radar as above, but with ground scatter plotted in grey, date format as `hh:mm`, custom min and max values and a colour bar label could look something like: ```python pydarn.RTP.plot_range_time(fitacf_data, beam_num=fitacf_data[0]['bmnum'], groundscatter=True, zmax=500, zmin=-500, date_fmt='%H:%M', colorbar_label='Line-of-Sight Velocity (m s$^{-1}$)', - coords=pydarn.Coords.RANGE_GATE) + range_estimation=pydarn.RangeEstimation.RANGE_GATE) ``` which outputs: @@ -117,7 +117,7 @@ Because the default parameter plotted is line-of-sight velocity, there is also a To change the colormap, use the 'cmap' parameter with the string name of a matplotlib color map ([found here](https://matplotlib.org/tutorials/colors/colormaps.html)). For example, plotting the power along the beam above using the colormap 'viridis': ```python -pydarn.RTP.plot_range_time(fitacf_data, beam_num=7, parameter='p_l', zmax=50, zmin=0, date_fmt='%H%M', colorbar_label='Power (dB)', coords=pydarn.Coords.RANGE_GATE, cmap='viridis') +pydarn.RTP.plot_range_time(fitacf_data, beam_num=7, parameter='p_l', zmax=50, zmin=0, date_fmt='%H%M', colorbar_label='Power (dB)', range_estimation=pydarn.RangeEstimation.RANGE_GATE, cmap='viridis') ``` produces: diff --git a/docs/user/summary.md b/docs/user/summary.md index 76eff0bd..d5230e55 100644 --- a/docs/user/summary.md +++ b/docs/user/summary.md @@ -44,7 +44,7 @@ darn_read = pydarn.SuperDARNRead(fitacf_file) fitacf_data = darn_read.read_fitacf() pydarn.RTP.plot_summary(fitacf_data, beam_num=2, - coords=pydarn.Coords.RANGE_GATE) + range_estimation=pydarn.RangeEstimation.RANGE_GATE) plt.show() ``` which gives: @@ -53,12 +53,12 @@ which gives: Note that ground scatter is displayed as default in the velocity panel. To disable it, you can use the option `groundscatter=false`. -In addition, we use `coords=pydarn.Coords.RANGE_GATE` to set the y-axis to plot in range gates. +In addition, we use `range_estimation=pydarn.RangeEstimation.RANGE_GATE` to set the y-axis to plot in range gates. Ground-Scatter Mapped Range is another type of axis you can use with range-time plots: ```python pydarn.RTP.plot_summary(fitacf_data, beam_num=3, - coords=pydarn.Coords.GROUND_SCATTER_MAPPED_RANGE) + range_estimation=pydarn.RangeEstimation.GSMR) ``` ![](../imgs/summary_1.png) @@ -86,5 +86,5 @@ Other common options include: | channel=(int) | Specify channel number (default=all) | | watermark=(bool) | True adds a 'not for publication' watermark | | cmaps=(dict/str) | Specifies the colour maps used in plotting | -| coords=(Coords) | Coordinates to use for the y-axis (See [Coordinates](coordinates.md)) | +| range_estimation=(RangeEstimation) | Coordinates to use for the y-axis (See [Coordinates](coordinates.md)) | For more options on how to modify plot_summary, take a look at the method in `rtp.py`. diff --git a/mkdocs.yml b/mkdocs.yml index 9fd0b689..dcdb4974 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -22,6 +22,7 @@ nav: - FOV plots: user/fov.md - Fan plots: user/fan.md - Grid plots: user/grid.md + - Convection Map plots: user/map.md - Power plots: user/power.md - ACF plots: user/acf.md - Logging: user/logging.md diff --git a/pydarn/__init__.py b/pydarn/__init__.py index d3cf6d88..e95d2d57 100644 --- a/pydarn/__init__.py +++ b/pydarn/__init__.py @@ -1,11 +1,18 @@ +#Copyright 2018 SuperDARN Canada, Univeristy Saskatchewan +#Author(s): Marina Schmidt +# +# Disclaimer: +# pyDARN is under the LGPL v3 license found in the root directory LICENSE.md +# Everyone is permitted to copy and distribute verbatim copies of this license +# document, but changing it is not allowed. +# +# This version of the GNU Lesser General Public License incorporates the terms +# and conditions of version 3 of the GNU General Public License, +# supplemented by the additional permissions listed below. +# +# Modifications: +# 2022-03-10 MTS - removed radar_fov from the __init__ file """ -Copyright 2018 SuperDARN Canada, Univeristy Saskatchewan -Author(s): Marina Schmidt - -Licensed under GNU v3.0 - -__init__.py -2018-11-05 Init file to setup the logging configuration and linking pyDARN's module, classes, and functions. """ @@ -27,33 +34,37 @@ from .exceptions.warning_formatting import only_message_warning_format from .exceptions.warning_formatting import citing_warning from .exceptions.warning_formatting import partial_record_warning +from .exceptions.warning_formatting import cartopy_warning +from .exceptions.warning_formatting import cartopy_print_warning # importing utils -from .utils.coordinates import Coords from .utils.constants import Re from .utils.constants import EARTH_EQUATORIAL_RADIUS from .utils.constants import C +from .utils.range_estimations import RangeEstimation +from .utils.virtual_heights import VHModels from .utils.conversions import dmap2dict -from .utils.conversions import gate2slant -from .utils.conversions import gate2GroundScatter +from .utils.plotting import MapParams from .utils.plotting import check_data_type from .utils.plotting import time2datetime +from .utils.plotting import find_record from .utils.superdarn_radars import SuperDARNRadars from .utils.superdarn_cpid import SuperDARNCpids from .utils.superdarn_radars import Hemisphere from .utils.superdarn_radars import read_hdw_file from .utils.superdarn_radars import get_hdw_files from .utils.scan import build_scan -from .utils.radar_pos import radar_fov -from .utils.radar_pos import geographic_cell_positions +from .utils.geo import geocentric_coordinates +from .utils.coordinates import Coords # import plotting from .plotting.color_maps import PyDARNColormaps -from .plotting.axis import Projections +from .plotting.projections import Projs from .plotting.rtp import RTP from .plotting.fan import Fan from .plotting.grid import Grid from .plotting.acf import ACF from .plotting.power import Power +from .plotting.maps import Maps citing_warning() diff --git a/pydarn/exceptions/plot_exceptions.py b/pydarn/exceptions/plot_exceptions.py index 812fedff..c102a280 100644 --- a/pydarn/exceptions/plot_exceptions.py +++ b/pydarn/exceptions/plot_exceptions.py @@ -3,6 +3,8 @@ # # Modifications: # 20210909: CJM - Added NoChannelError +# 20220308 MTS - Added PartialRecordsError() +# 2022-03-23 MTS - Added NotImplementedError() indicates when something is not implement in pyDARN # # Disclaimer: # pyDARN is under the LGPL v3 license found in the root directory LICENSE.md @@ -19,6 +21,52 @@ pydarn_log = logging.getLogger('pydarn') +class PartialRecordsError(Exception): + """ + Error given when a partial record prevents a plotting algorithm to work + """ + def __init__(self, missing_field): + self.message = "pyDARN could not plot the following due to the "\ + "following field {} was missing. This may be due to partial"\ + "records created in RST (data processing software for"\ + " SuperDARN)".format(missing_field) + super().__init__(self.message) + pydarn_log.error(self.message) + + +class CartopyMissingError(Exception): + """ + Error given when attmpting a cartopy style plot but cartopy isn't installed + """ + def __init__(self): + self.message = 'cartopy is independent library that the user must'\ + 'install to use. Please either change '\ + 'plot styles or install cartopy and dependencies: '\ + 'https://pydarn.readthedocs.io/en/main/user/install/.' + super().__init__(self.message) + pydarn_log.error(self.message) + + +class CartopyVersionError(Exception): + """ + Error given when attmpting a cartopy style plot but cartopy isn't installed + """ + def __init__(self, version): + self.message = 'cartopy is independent library that the user must'\ + 'install to use. Please insure the version number is '\ + '>= 0.19. Your current version is: {}'.format(version) + super().__init__(self.message) + pydarn_log.error(self.message) + +class NotImplemented(Exception): + """ + Error given when settings are not implemented right now. + """ + def __init__(self, msg): + self.message = msg + super().__init__(self.message) + pydarn_log.error(self.message) + class IncorrectPlotMethodError(Exception): """ This error is raised when the incorrect plotting @@ -192,6 +240,6 @@ def __init__(self, channel: int, opt_channel: int = None): self.message = "There is no data for channel {channel},"\ " try channel {opt_channel}."\ .format(channel=self.channel, - opt_channel=self.opt_channel) + opt_channel=self.opt_channel) super().__init__(self.message) - pydarn_log.error(self.message) \ No newline at end of file + pydarn_log.error(self.message) diff --git a/pydarn/exceptions/radar_exceptions.py b/pydarn/exceptions/radar_exceptions.py index 4e79890c..98b6b8c1 100644 --- a/pydarn/exceptions/radar_exceptions.py +++ b/pydarn/exceptions/radar_exceptions.py @@ -27,3 +27,15 @@ def __init__(self, abbrev): "".format(abv=self.abbreviation, path=self.path) super().__init__(self.message) pydarn_log.error(self.message) + + +class RangeEstimationError(Exception): + """ + This error is raised when there is an incorrect range estimation + used. + """ + def __init__(self, msg): + self.message = msg + super().__init__(self.message) + pydarn_log.error(self.message) + diff --git a/pydarn/exceptions/warning_formatting.py b/pydarn/exceptions/warning_formatting.py index 75da5a5a..ccd88a20 100644 --- a/pydarn/exceptions/warning_formatting.py +++ b/pydarn/exceptions/warning_formatting.py @@ -19,6 +19,30 @@ def citing_warning(): " for SuperDARN data is found at" " https://pydarn.readthedocs.io/en/master/user/citing/") +def cartopy_print_warning(): + """ + This warning prints on installation to inform users they + need to install cartopy and should read the docs before + using pyDARN if they wish to use geographical layouts. + """ + print() # add a new line to break up more of the text + print("IMPORTANT!: If you are going to use Fan, Grid, " + "and/or Convection Map " + "plots, then make sure cartopy is installed on your machine. " + "If you do not need to use cartopy for your plotting, ignore " + "this message.") + +def cartopy_warning(): + """ + This warning prints on installation to inform users they + need to install cartopy and should read the docs before + using pyDARN if they wish to use geographical layouts. + """ + warnings.warn("If you are going to use Fan, Grid, and/or Convection Map " + "plots, then make sure cartopy is installed on your machine. " + "If you do not need to use cartopy for your plotting, ignore " + "this message.") + def partial_record_warning(): """ prints a warning that the data chosen to be plotted is missing some diff --git a/pydarn/plotting/axis.py b/pydarn/plotting/axis.py deleted file mode 100644 index 530e9ac2..00000000 --- a/pydarn/plotting/axis.py +++ /dev/null @@ -1,75 +0,0 @@ -# Copyright (C) 2021 SuperDARN Canada, University of Saskatchewan -# Author: Daniel Billett -# -# Modifications: -# -# Disclaimer: -# pyDARN is under the LGPL v3 license found in the root directory LICENSE.md -# Everyone is permitted to copy and distribute verbatim copies of this license -# document, but changing it is not allowed. -# -# This version of the GNU Lesser General Public License incorporates the terms -# and conditions of version 3 of the GNU General Public License, -# supplemented by the additional permissions listed below. -""" -Code which generates axis objects for use in plotting functions -""" - -import matplotlib.pyplot as plt -import numpy as np - -from pydarn import Hemisphere - - -class Projections(): - """ - Methods for generating axis objects for use in plotting functions - Methods - ------- - axis_polar - """ - - def __str__(self): - return "This class is static class that provides"\ - " the following methods: \n"\ - " - axis_polar()\n" - - @classmethod - def axis_polar(cls, lowlat: int = 30, hemisphere: object = Hemisphere.North): - """ - Plots a radar's Field Of View (FOV) fan plot for the given data and - scan number - - Parameters - ----------- - lowlat: int - Lower AACGM latitude boundary for the polar plot - Default: 30 - hemiphere: enum - Hemisphere of polar projection. Can be Hemisphere.North or - Hemisphere.South for northern and southern hemispheres, - respectively - Default: Hemisphere.North - """ - - ax = plt.axes(polar=True) - - # Set upper and lower latitude limits (pole and lowlat) - if hemisphere == Hemisphere.North: - ax.set_ylim(90, lowlat) - ax.set_yticks(np.arange(lowlat, 90, 10)) - else: - # If hemisphere is South, lowlat must be negative - ax.set_ylim(-90, -abs(lowlat)) - ax.set_yticks(np.arange(-abs(lowlat), -90, -10)) - - # Locations of tick marks. Will be customisable in future - ax.set_xticks([0, np.radians(45), np.radians(90), np.radians(135), - np.radians(180), np.radians(225), np.radians(270), - np.radians(315)]) - - # Tick labels will depend on coordinate system - ax.set_xticklabels(['00', '', '06', '', '12', '', '18', '']) - ax.set_theta_zero_location("S") - - return ax diff --git a/pydarn/plotting/fan.py b/pydarn/plotting/fan.py index 5d5627d6..e387e4f7 100644 --- a/pydarn/plotting/fan.py +++ b/pydarn/plotting/fan.py @@ -9,7 +9,15 @@ # 2021-09-09: CJM - Included a channel option for plot_fan # 2021-09-08: CJM - Included individual gate and beam boundary plotting for FOV # 2021-11-22: MTS - pass in axes object to plot_fov -# +# 2021-11-18: MTS - Added Projectsion class for cartopy use +# 2021-02-02: CJM - Included rsep and frang options in plot_fov +# - Re-arranged logic for ranges option +# - Indexing fixed for give lower ranges that are non-zero +# 2022-03-10: MTS - switched coords involving range estimations to RangeEstimation +# - Removed GEO_COASTALINE, replaced with GEO on projections +# - reduced some parameters inputs to kwargs +# 2022-03-22: CJM - Set cmap bad values to transparent +# 2022-03-23: MTS - added the NotImplementedError for AACGM and GEO projection as this has yet to be figured out # Disclaimer: # pyDARN is under the LGPL v3 license found in the root directory LICENSE.md # Everyone is permitted to copy and distribute verbatim copies of this license @@ -34,23 +42,22 @@ # Third party libraries import aacgmv2 -from pydarn import (PyDARNColormaps, build_scan, radar_fov, - time2datetime, plot_exceptions, Coords, - SuperDARNRadars, Hemisphere, Projections, - partial_record_warning) +from pydarn import (PyDARNColormaps, build_scan, partial_record_warning, + time2datetime, plot_exceptions, SuperDARNRadars, + Projs, Coords, Hemisphere, RangeEstimation) class Fan(): """ - 'Fan', or 'Field-of-view' plots for SuperDARN FITACF data - This class inherits from matplotlib to generate the figures - Methods - ------- - plot_fan - plot_fov - plot_radar_position - plot_radar_label - """ + 'Fan', or 'Field-of-view' plots for SuperDARN FITACF data + This class inherits from matplotlib to generate the figures + Methods + ------- + plot_fan + plot_fov + plot_radar_position + plot_radar_label + """ def __str__(self): return "This class is static class that provides"\ @@ -59,13 +66,15 @@ def __str__(self): " - plot_fov()\n" @classmethod - def plot_fan(cls, dmap_data: List[dict], ax=None, + def plot_fan(cls, dmap_data: List[dict], ax=None, ranges: List = [], scan_index: Union[int, dt.datetime] = 1, parameter: str = 'v', cmap: str = None, groundscatter: bool = False, zmin: int = None, zmax: int = None, colorbar: bool = True, colorbar_label: str = '', title: bool = True, - channel = 'all', **kwargs): + boundary: bool = True, projs: Projs = Projs.POLAR, + coords: Coords = Coords.AACGM_MLT, + channel: int = 'all', **kwargs): """ Plots a radar's Field Of View (FOV) fan plot for the given data and scan number @@ -116,6 +125,15 @@ def plot_fan(cls, dmap_data: List[dict], ax=None, if true then will create a title, else user can define it with plt.title default: true + boundary: bool + if true then plots the FOV boundaries + default: true + projs: Enum + choice of projection for plot + default: Projs.POLAR (polar projection) + coords: Enum + choice of plotting coordinates + default: Coords.AACGM_MLT (Magnetic Lat and MLT) channel : int or str integer indicating which channel to plot or 'all' to plot all channels @@ -127,9 +145,11 @@ def plot_fan(cls, dmap_data: List[dict], ax=None, Returns ----------- beam_corners_aacgm_lats - n_beams x n_gates numpy array of AACGMv2 latitudes + n_beams x n_gates numpy array of latitudes + return values dependent on given coords enum beam_corners_aacgm_lons - n_beams x n_gates numpy array of AACGMv2 longitudes + n_beams x n_gates numpy array of longitudes or MLT + return values dependent on given coords enum scan n_beams x n_gates numpy array of the scan data (for the selected parameter) @@ -147,12 +167,8 @@ def plot_fan(cls, dmap_data: List[dict], ax=None, dmap_data = [rec for rec in dmap_data if rec['channel'] == channel] # If no records exist, advise user that the channel is not used if not dmap_data: - raise plot_exceptions.NoChannelError(channel,opt_channel) + raise plot_exceptions.NoChannelError(channel, opt_channel) - try: - ranges = kwargs['ranges'] - except KeyError: - ranges = [0, 75] # Get scan numbers for each record beam_scan = build_scan(dmap_data) scan_time = None @@ -163,7 +179,7 @@ def plot_fan(cls, dmap_data: List[dict], ax=None, scan_index = 0 found_match = False for rec in dmap_data: - rec_time = time2datetime(rec) + rec_time = time2datetime(rec) if abs(rec['scan']) == 1: scan_index += 1 # Need the abs since you cannot have negative seconds @@ -172,31 +188,70 @@ def plot_fan(cls, dmap_data: List[dict], ax=None, found_match = True break # handle datetimes out of bounds - if found_match == False: + if found_match is False: raise plot_exceptions.IncorrectDateError(rec_time, scan_time) # Locate scan in loaded data plot_beams = np.where(beam_scan == scan_index) - # Time for coordinate conversion if not scan_time: - date = time2datetime(dmap_data[plot_beams[0][0]]) + date = time2datetime(dmap_data[plot_beams[0][0]]) else: - date = scan_time + date = scan_time # Plot FOV outline - if ranges is None: - ranges = [0, dmap_data[0]['nrang']] + stid = dmap_data[0]['stid'] + if ranges == [] or ranges is None: + try: + # If not given, get ranges from data file + ranges = [0, dmap_data[0]['nrang']] + except KeyError: + # Otherwise, default to [0,75] + ranges = [0, SuperDARNRadars.radars[stid].range_gate_45] - beam_corners_aacgm_lats, beam_corners_aacgm_lons, thetas, rs, ax = \ - cls.plot_fov(dmap_data[0]['stid'], date, ax=ax, **kwargs) + # Get rsep and frang from data unless not there then take defaults + # of 180 km for frang and 45 km for rsep as these are most commonly + # used + try: + frang = dmap_data[0]['frang'] + except KeyError: + frang = 180 - fan_shape = beam_corners_aacgm_lons.shape + try: + rsep = dmap_data[0]['rsep'] + except KeyError: + rsep = 45 + + if coords != Coords.GEOGRAPHIC and projs == Projs.GEO: + raise plot_exceptions.NotImplemented("AACGM coordinates are"\ + " not implemented for "\ + " geographic projections right"\ + " now, if you would like to"\ + " see it sooner please help"\ + " out at "\ + "https://github.com"\ + "/SuperDARN/pyDARN") + + beam_corners_lats, beam_corners_lons =\ + coords(stid=dmap_data[0]['stid'], + rsep=rsep, frang=frang, + gates=ranges, date=date, + **kwargs) + + fan_shape = beam_corners_lons.shape + if ranges[0] < ranges[1] - fan_shape[0]: + ranges[0] = ranges[1] - fan_shape[0] + 1 + rs = beam_corners_lats + if projs != Projs.GEO: + thetas = np.radians(beam_corners_lons) + else: + thetas = beam_corners_lons # Get range-gate data and groundscatter array for given scan - scan = np.zeros((fan_shape[0] - 1, fan_shape[1]-1)) - grndsct = np.zeros((fan_shape[0] - 1, fan_shape[1]-1)) - + # fan_shape has no -1 as when given ranges we want to include the + # both ends for the ranges given + scan = np.zeros((fan_shape[0], fan_shape[1]-1)) + grndsct = np.zeros((fan_shape[0], fan_shape[1]-1)) # Colour table and max value selection depending on parameter plotted # Load defaults if none given if cmap is None: @@ -206,6 +261,10 @@ def plot_fan(cls, dmap_data: List[dict], ax=None, 'elv': PyDARNColormaps.PYDARN} cmap = plt.cm.get_cmap(cmap[parameter]) + # Set background to transparent - avoids carry over + # does not interfere with the fov color if chosen + cmap.set_bad(alpha=0.0) + # Setting zmin and zmax defaultzminmax = {'p_l': [0, 50], 'v': [-200, 200], 'w_l': [0, 250], 'elv': [0, 50]} @@ -227,37 +286,58 @@ def plot_fan(cls, dmap_data: List[dict], ax=None, # This is a temporary fix to manage inconsistencies between the # fitacf files and the hardware files. The issue will be # fully resolved when the `rpos` code is committed. - good_data=slist<(fan_shape[0] - 1) - slist=slist[good_data] - temp_data=dmap_data[i.astype(int)][parameter][good_data] - temp_ground=dmap_data[i.astype(int)]['gflg'][good_data] - - scan[slist, beam] = temp_data - grndsct[slist, beam] = temp_ground + good_data = np.where((slist >= ranges[0]) & + (slist < ranges[1])) + slist = slist[good_data] + temp_data = dmap_data[i.astype(int)][parameter][good_data] + temp_ground = dmap_data[i.astype(int)]['gflg'][good_data] + + scan[slist-ranges[0], beam] = temp_data + grndsct[slist-ranges[0], beam] = temp_ground # if there is no slist field this means partial record except KeyError: partial_record_warning() continue + # Begin plotting by iterating over ranges and beams - thetas = thetas[ranges[0]:ranges[1]] - rs = rs[ranges[0]:ranges[1]] - scan = scan[ranges[0]:ranges[1]-1] + + thetas = thetas[0:ranges[1]-ranges[0]+1] + rs = rs[0:ranges[1]-ranges[0]+1] + + scan = scan[0:ranges[1]-ranges[0]] + grndsct = grndsct[0:ranges[1]-ranges[0]] + # Set up axes in correct hemisphere + stid = dmap_data[0]['stid'] + kwargs['hemisphere'] = SuperDARNRadars.radars[stid].hemisphere + + ax, ccrs = projs(date=date, **kwargs) + + if ccrs is None: + transform = ax.transData + + else: + transform = ccrs.PlateCarree() ax.pcolormesh(thetas, rs, np.ma.masked_array(scan, ~scan.astype(bool)), - norm=norm, cmap=cmap) + norm=norm, cmap=cmap, transform=transform, + zorder=2) # plot the groundscatter as grey fill if groundscatter: - grndsct = grndsct[ranges[0]:ranges[1]-1] ax.pcolormesh(thetas, rs, np.ma.masked_array(grndsct, ~grndsct.astype(bool)), - norm=norm, cmap='Greys') + norm=norm, cmap='Greys', transform=transform) + if ccrs is None: + azm = np.linspace(0, 2 * np.pi, 100) + r, th = np.meshgrid(rs, azm) + plt.plot(azm, r, color='k', ls='none') + plt.grid() - azm = np.linspace(0, 2 * np.pi, 100) - r, th = np.meshgrid(rs, azm) - plt.plot(azm, r, color='k', ls='none') - plt.grid() + if boundary: + cls.plot_fov(stid=dmap_data[0]['stid'], date=date, ax=ax, + ccrs=ccrs, coords=coords, projs=projs, rsep=rsep, + frang=frang, ranges=ranges, **kwargs) # Create color bar if True if colorbar is True: @@ -276,16 +356,19 @@ def plot_fan(cls, dmap_data: List[dict], ax=None, end_time = time2datetime(dmap_data[plot_beams[-1][-1]]) title = cls.__add_title__(start_time, end_time) plt.title(title) - return beam_corners_aacgm_lats, beam_corners_aacgm_lons, scan, grndsct + return ax, beam_corners_lats, beam_corners_lons, scan, grndsct @classmethod def plot_fov(cls, stid: str, date: dt.datetime, - ax=None, ranges: List = [], boundary: bool = True, + ax=None, ccrs=None, ranges: List = [], boundary: bool = True, + rsep: int = 45, frang: int = 180, + projs: object = Projs.POLAR, + coords: Coords = Coords.AACGM_MLT, fov_color: str = None, alpha: int = 0.5, radar_location: bool = True, radar_label: bool = False, line_color: str = 'black', grid: bool = False, - line_alpha: int = 0.5 , **kwargs): + line_alpha: int = 0.5, **kwargs): """ plots only the field of view (FOV) for a given radar station ID (stid) @@ -305,6 +388,9 @@ def plot_fov(cls, stid: str, date: dt.datetime, Set to a two element list of the lower and upper ranges to plot If None, the max will be obtained by SuperDARNRadars Default: None + projs: Pojs object + Sets the projection type for the plot + Default: Projs.POLAR boundary: bool Set to false to not plot the outline of the FOV Default: True @@ -339,96 +425,137 @@ def plot_fov(cls, stid: str, date: dt.datetime, ------- beam_corners_aacgm_lats - list of beam corners AACGM latitudes beam_corners_aacgm_lons - list of beam corners AACGM longitudes - thetas - theta polar coordinates + beam_corners_lon - theta polar coordinates rs - radius polar coordinates """ if ranges == [] or ranges is None: ranges = [0, SuperDARNRadars.radars[stid].range_gate_45] - # Get radar beam/gate locations - beam_corners_aacgm_lats, beam_corners_aacgm_lons = \ - radar_fov(stid, ranges=ranges, date=date, **kwargs) - if not date: date = dt.datetime.now() - fan_shape = beam_corners_aacgm_lons.shape - if ranges == []: - ranges = [0, fan_shape[0]] - # Work out shift due in MLT - beam_corners_mlts = np.zeros((fan_shape[0], fan_shape[1])) - mltshift = beam_corners_aacgm_lons[0, 0] - \ - (aacgmv2.convert_mlt(beam_corners_aacgm_lons[0, 0], date) * 15) - beam_corners_mlts = beam_corners_aacgm_lons - mltshift + # Get radar beam/gate locations + beam_corners_lats, beam_corners_lons = \ + coords(stid=stid, gates=ranges, rsep=rsep, frang=frang, + date=date, **kwargs) - # Hold the beam positions - thetas = np.radians(beam_corners_mlts) - rs = beam_corners_aacgm_lats + if projs == Projs.POLAR: + beam_corners_lons = np.radians(beam_corners_lons) # Setup plot # This may screw up references + hemisphere = SuperDARNRadars.radars[stid].hemisphere if ax is None: # Get the hemisphere to pass to plotting projection - kwargs['hemisphere'] = SuperDARNRadars.radars[stid].hemisphere - # Get a polar projection using any kwarg input - ax = Projections.axis_polar(**kwargs) + kwargs['hemisphere'] = hemisphere + ax, ccrs = projs(date=date, **kwargs) + if ccrs is None: + transform = ax.transData + else: + transform = ccrs.Geodetic() if boundary: # left boundary line - plt.plot(thetas[0:ranges[1], 0], rs[0:ranges[1], 0], + plt.plot(beam_corners_lons[0:ranges[1]-ranges[0]+1, 0], + beam_corners_lats[0:ranges[1]-ranges[0]+1, 0], color=line_color, linewidth=0.5, - alpha=line_alpha) + alpha=line_alpha, transform=transform, zorder=3) # top radar arc - plt.plot(thetas[ranges[1] - 1, 0:thetas.shape[1]], - rs[ranges[1] - 1, 0:thetas.shape[1]], + plt.plot(beam_corners_lons[ranges[1]-ranges[0], + 0:beam_corners_lons.shape[1]], + beam_corners_lats[ranges[1]-ranges[0], + 0:beam_corners_lons.shape[1]], color=line_color, linewidth=0.5, - alpha=line_alpha) + alpha=line_alpha, transform=transform, zorder=3) # right boundary line - plt.plot(thetas[0:ranges[1], thetas.shape[1] - 1], - rs[0:ranges[1], thetas.shape[1] - 1], + plt.plot(beam_corners_lons[0:ranges[1]-ranges[0]+1, + beam_corners_lons.shape[1] - 1], + beam_corners_lats[0:ranges[1]-ranges[0]+1, + beam_corners_lons.shape[1] - 1], color=line_color, linewidth=0.5, - alpha=line_alpha) + alpha=line_alpha, transform=transform, zorder=3) # bottom arc - plt.plot(thetas[0, 0:thetas.shape[1] - 1], - rs[0, 0:thetas.shape[1] - 1], color=line_color, - linewidth=0.5, alpha=line_alpha) + plt.plot(beam_corners_lons[0, 0:beam_corners_lons.shape[1]], + beam_corners_lats[0, 0:beam_corners_lons.shape[1]], + color=line_color, linewidth=0.5, alpha=line_alpha, + transform=transform, zorder=3) + + fan_shape = beam_corners_lons.shape if grid: # This plots lines along the beams for bm in range(fan_shape[1]): - plt.plot(thetas[0:ranges[1], bm - 1], - rs[0:ranges[1], bm - 1], - color=line_color, linewidth=0.2, - alpha=line_alpha) + plt.plot(beam_corners_lons[0:ranges[1] + 1, bm - 1], + beam_corners_lats[0:ranges[1] + 1, bm - 1], + color=line_color, linewidth=0.2, + alpha=line_alpha, transform=transform, + zorder=3) # This plots arcs along the gates - for g in range(ranges[1]): - plt.plot(thetas[g-1, 0:thetas.shape[1]], - rs[g - 1, 0:thetas.shape[1]], - color=line_color, linewidth=0.2, - alpha=line_alpha) + for g in range(ranges[1] - ranges[0] + 1): + plt.plot(beam_corners_lons[g - 1, + 0:beam_corners_lons.shape[1]], + beam_corners_lats[g - 1, + 0:beam_corners_lons.shape[1]], + color=line_color, linewidth=0.2, + alpha=line_alpha, transform=transform, + zorder=3) if radar_location: - cls.plot_radar_position(stid, date, line_color, **kwargs) + cls.plot_radar_position(stid, date=date, line_color=line_color, + transform=transform, projs=projs, + coords=coords, ccrs=ccrs, **kwargs) if radar_label: - cls.plot_radar_label(stid, date, line_color, **kwargs) - + cls.plot_radar_label(stid, date=date, line_color=line_color, + transform=transform, projs=projs, + coords=coords, ccrs=ccrs, **kwargs) if fov_color is not None: - theta = thetas[0:ranges[1], 0] - theta = np.append(theta, thetas[ranges[1]-1, 0:thetas.shape[1]-1]) - theta = np.append(theta, np.flip(thetas[0:ranges[1], - thetas.shape[1]-1])) - theta = np.append(theta, np.flip(thetas[0, 0:thetas.shape[1]-1])) - - r = rs[0:ranges[1], 0] - r = np.append(r, rs[ranges[1]-1, 0:thetas.shape[1]-1]) - r = np.append(r, np.flip(rs[0:ranges[1], thetas.shape[1]-1])) - r = np.append(r, np.flip(rs[0, 0:thetas.shape[1]-1])) - ax.fill(theta, r, color=fov_color, alpha=alpha, zorder=0) - return beam_corners_aacgm_lats, beam_corners_aacgm_lons, thetas, rs, ax + theta = beam_corners_lons[0:ranges[1] + 1, 0] + theta =\ + np.append(theta, + beam_corners_lons[ranges[1]-ranges[0], + 0:beam_corners_lons.shape[1]-1]) + theta =\ + np.append(theta, + np.flip(beam_corners_lons[0:ranges[1] - + ranges[0] + 1, + beam_corners_lons. + shape[1] - 1])) + theta =\ + np.append(theta, + np.flip(beam_corners_lons[0, + 0:beam_corners_lons. + shape[1] - 1])) + + r = beam_corners_lats[0:ranges[1] + 1, 0] + r =\ + np.append(r, beam_corners_lats[ranges[1]-ranges[0], + 0:beam_corners_lons. + shape[1] - 1]) + r =\ + np.append(r, + np.flip(beam_corners_lats[0:ranges[1] - ranges[0]+1, + beam_corners_lons. + shape[1] - 1])) + r =\ + np.append(r, + np.flip(beam_corners_lats[0, + 0:beam_corners_lons. + shape[1] - 1])) + + theta = np.flip(theta) + r = np.flip(r) + + ax.fill(theta, r, color=fov_color, alpha=alpha, zorder=1, + transform=transform) + + return beam_corners_lats, beam_corners_lons, ax, ccrs @classmethod def plot_radar_position(cls, stid: int, date: dt.datetime, + transform: object = None, + coords: Coords = Coords.AACGM_MLT, + projs: Projs = Projs.POLAR, line_color: str = 'black', **kwargs): """ plots only a dot at the position of a given radar station ID (stid) @@ -452,19 +579,26 @@ def plot_radar_position(cls, stid: int, date: dt.datetime, lat = SuperDARNRadars.radars[stid].hardware_info.geographic.lat lon = SuperDARNRadars.radars[stid].hardware_info.geographic.lon # Convert to geomag coords - geomag_radar = aacgmv2.get_aacgm_coord(lat, lon, 250, date) - mltshift = geomag_radar[1] - (aacgmv2.convert_mlt(geomag_radar[1], - date) * 15) - # Convert to MLT then radians for theta - theta_lon = np.radians(geomag_radar[1] - mltshift) - r_lat = geomag_radar[0] + if coords == Coords.AACGM_MLT or coords == Coords.AACGM: + geomag_radar = aacgmv2.get_aacgm_coord(lat, lon, 250, date) + lat = geomag_radar[0] + lon = geomag_radar[1] + if coords == Coords.AACGM_MLT: + mltshift = geomag_radar[1] -\ + (aacgmv2.convert_mlt(geomag_radar[1], date) * 15) + lon = geomag_radar[1] - mltshift + if projs == Projs.POLAR: + lon = np.radians(lon) # Plot a dot at the radar site - plt.scatter(theta_lon, r_lat, c=line_color, s=5) + plt.scatter(lon, lat, c=line_color, s=5, transform=transform) return @classmethod def plot_radar_label(cls, stid: int, date: dt.datetime, - line_color: str = 'black', **kwargs): + coords: Coords = Coords.AACGM_MLT, + projs: Projs = Projs.POLAR, + line_color: str = 'black', transform: object = None, + **kwargs): """ plots only string at the position of a given radar station ID (stid) @@ -489,21 +623,27 @@ def plot_radar_label(cls, stid: int, date: dt.datetime, # Get location of radar lat = SuperDARNRadars.radars[stid].hardware_info.geographic.lat lon = SuperDARNRadars.radars[stid].hardware_info.geographic.lon + # Convert to geomag coords - geomag_radar = aacgmv2.get_aacgm_coord(lat, lon, 250, date) - mltshift = geomag_radar[1] - \ - (aacgmv2.convert_mlt(geomag_radar[1], date) * 15) - # Convert to MLT then radians for theta - theta_lon = np.radians(geomag_radar[1] - mltshift) - r_lat = geomag_radar[0] - - theta_text = theta_lon + if coords == Coords.AACGM_MLT or coords == Coords.AACGM: + geomag_radar = aacgmv2.get_aacgm_coord(lat, lon, 250, date) + lat = geomag_radar[0] + lon = geomag_radar[1] + if coords == Coords.AACGM_MLT: + mltshift = geomag_radar[1] -\ + (aacgmv2.convert_mlt(geomag_radar[1], date) * 15) + lon = geomag_radar[1] - mltshift + if projs == Projs.POLAR: + lon = np.radians(lon) + + theta_text = lon # Shift in latitude (dependent on hemisphere) if SuperDARNRadars.radars[stid].hemisphere == Hemisphere.North: - r_text = r_lat - 5 + r_text = lat - 5 else: - r_text = r_lat + 5 - plt.text(theta_text, r_text, label_str, ha='center', c=line_color) + r_text = lat + 5 + plt.text(theta_text, r_text, label_str, ha='center', + transform=transform, c=line_color) return @classmethod diff --git a/pydarn/plotting/grid.py b/pydarn/plotting/grid.py index 400ece82..52895c2a 100644 --- a/pydarn/plotting/grid.py +++ b/pydarn/plotting/grid.py @@ -2,6 +2,7 @@ # Author: Daniel Billett, Marina Schmidt # # Modifications: +# 20220308 MTS added partial record exception # # Disclaimer: # pyDARN is under the LGPL v3 license found in the root directory LICENSE.md @@ -31,8 +32,16 @@ from pydarn import (PyDARNColormaps, Fan, plot_exceptions, standard_warning_format) +try: + from cartopy.mpl import geoaxes + import cartopy.crs as ccrs + cartopyInstalled = True +except Exception: + cartopyInstalled = False + warnings.formatwarning = standard_warning_format + class Grid(): """ Grid plots for SuperDARN data @@ -159,11 +168,14 @@ def plot_grid(cls, dmap_data: List[dict], record: int = 0, with warnings.catch_warnings(): warnings.simplefilter("ignore") for stid in dmap_data[record]['stid']: - _, aacgm_lons, _, _, ax =\ + _, aacgm_lons, ax, _ =\ Fan.plot_fov(stid, date, ax=ax, **kwargs) - data_lons = dmap_data[record]['vector.mlon'] - data_lats = dmap_data[record]['vector.mlat'] + try: + data_lons = dmap_data[record]['vector.mlon'] + data_lats = dmap_data[record]['vector.mlat'] + except KeyError: + raise plot_exceptions.PartialRecordsError('vector.mlon') # Hold the beam positions shifted_mlts = aacgm_lons[0, 0] - \ diff --git a/pydarn/plotting/maps.py b/pydarn/plotting/maps.py new file mode 100644 index 00000000..683c3c28 --- /dev/null +++ b/pydarn/plotting/maps.py @@ -0,0 +1,927 @@ +# Copyright (C) 2021 SuperDARN Canada, University of Saskatchewan +# Author: Marina Schmidt +# Copyright (C) 2012 VT SuperDARN Lab +# Modifications: +# 2022-03-08: MTS - added partial records exception +# 2021-03-18: CJM - Included contour plotting and HMB +# 2021-03-31: CJM - Map info included +# 2021-03-31: CJM - IMF clock angle dial added +# 2021-04-01: CJM - Bug fix for lon shifting to MLT +# +# Disclaimer: +# pyDARN is under the LGPL v3 license found in the root directory LICENSE.md +# Everyone is permitted to copy and distribute verbatim copies of this license +# document, but changing it is not allowed. +# +# This version of the GNU Lesser General Public License incorporates the terms +# and conditions of version 3 of the GNU General Public License, +# supplemented by the additional permissions listed below. + + +""" +Grid plots, mapped to AACGM coordinates in a polar format +""" + +import datetime as dt +import matplotlib.pyplot as plt +import numpy as np +import warnings + +from enum import Enum +from matplotlib import ticker, cm, colors +from mpl_toolkits.axes_grid.inset_locator import InsetPosition +from scipy import special +from typing import List + +# Third party libraries +import aacgmv2 + +from pydarn import (PyDARNColormaps, plot_exceptions, + standard_warning_format, Re, Hemisphere, + time2datetime, find_record, Fan, MapParams) + +warnings.formatwarning = standard_warning_format + + +class Maps(): + """ + Maps plots for SuperDARN data + + Methods + ------- + plot_maps + calculated_fitted_velocities + """ + + def __str__(self): + return "This class is static class that provides"\ + " the following methods: \n"\ + " - plot_maps()\n" + + + @classmethod + def plot_mapdata(cls, dmap_data: List[dict], ax=None, + parameter: Enum = MapParams.FITTED_VELOCITY, + record: int = 0, start_time: dt.datetime = None, + time_delta: float = 1, alpha: float = 1.0, + len_factor: float = 150, cmap: str = None, + colorbar: bool = True, colorbar_label: str = '', + title: str = '', zmin: float = None, zmax: float = None, + hmb: bool = True, boundary: bool = False, + radar_location: bool = False, map_info: bool = True, + imf_dial: bool = True, **kwargs): + """ + Plots convection maps data points and vectors + + Parameters + ---------- + dmap_data : dict[List] + DMAP data returned from pyDARN.SuperDARNRead or pyDARNio + ax: object + matplotlib axis object + parameters: enum + enum object to determine what to plot + default: MapParams.FITTED_VELOCITY + record: int + record number to plot + default: 0 + start_time: datetime.datetime + datetime object as the start time of the record to plot + if none then record will be used + default: none + time_delta: int + How close the start_time has to be start_time of the record + in minutes + default: 1 + cmap: matplotlib.cm + matplotlib colour map + https://matplotlib.org/tutorials/colors/colormaps.html + Default: Official pyDARN colour map for given parameter + MapParams.FITTED_VELOCITY: 'plasma', + MapParams.MODEL_VELOCITY: 'plasma', + MapParams.RAW_VELOCITY: 'plasma', + MapParams.POWER: 'plasma_r', + MapParams.SPECTRAL_WIDTH: PyDARNColormaps.PYDARN_VIRIDIS + zmin: int + The minimum parameter value for coloring + Default: MapParams.FITTED_VELOCITY: [0], + MapParams.MODEL_VELOCITY: [0], + MapParams.RAW_VELOCITY: [0], + MapParams.POWER: [0], + MapParams.SPECTRAL_WIDTH: [0] + zmax: int + The maximum parameter value for coloring + Default: MapParams.FITTED_VELOCITY: [1000], + MapParams.MODEL_VELOCITY: [1000], + MapParams.RAW_VELOCITY: [1000], + MapParams.POWER: [250], + MapParams.SPECTRAL_WIDTH: [250] + colorbar: bool + Draw a colourbar if True + Default: True + colorbar_label: str + The label that appears next to the colour bar. + Requires colorbar to be true + Default: '' + title: str + Adds a title to the plot. If no title is specified, + one will be provided + Default: '' + len_factor: float + Normalisation factor for the vectors, to control size on plot + Larger number means smaller vectors on plot + Default: 150.0 + imf_dial: bool + If True, draw an IMF dial of the magnetic field clock angle. + Default: True + kwargs: key=value + uses the parameters for plot_fov and projections.axis + + """ + # Find the record corresponding to the start time + if start_time is not None: + record = find_record(dmap_data, start_time, time_delta) + date = time2datetime(dmap_data[record]) + if cmap is None: + cmap = {MapParams.FITTED_VELOCITY: 'plasma_r', + MapParams.MODEL_VELOCITY: 'plasma_r', + MapParams.RAW_VELOCITY: 'plasma_r', + MapParams.POWER: 'plasma', + MapParams.SPECTRAL_WIDTH: PyDARNColormaps.PYDARN_VIRIDIS} + cmap = plt.cm.get_cmap(cmap[parameter]) + # Setting zmin and zmax + defaultzminmax = {MapParams.FITTED_VELOCITY: [0, 1000], + MapParams.MODEL_VELOCITY: [0, 1000], + MapParams.RAW_VELOCITY: [0, 1000], + MapParams.POWER: [0, 250], + MapParams.SPECTRAL_WIDTH: [0, 250]} + if zmin is None: + zmin = defaultzminmax[parameter][0] + if zmax is None: + zmax = defaultzminmax[parameter][1] + + hemisphere = Hemisphere(dmap_data[record]['hemisphere']) + norm = colors.Normalize + norm = norm(zmin, zmax) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + # TODO: Make FOV outlines optional or invisible + for stid in dmap_data[record]['stid']: + _, aacgm_lons, ax, _ =\ + Fan.plot_fov(stid, date, ax=ax, boundary=boundary, + radar_location=radar_location, **kwargs) + + if parameter == MapParams.MODEL_VELOCITY: + try: + data_lons = dmap_data[record]['model.mlon'] + data_lats = dmap_data[record]['model.mlat'] + except KeyError: + raise plot_exceptions.PartialRecordsError('model.mlon') + else: + try: + data_lons = dmap_data[record]['vector.mlon'] + data_lats = dmap_data[record]['vector.mlat'] + except KeyError: + raise plot_exceptions.PartialRecordsError('model.mlat') + + # Hold the beam positions + shifted_mlts = aacgm_lons[0, 0] - \ + (aacgmv2.convert_mlt(aacgm_lons[0, 0], date) * 15) + shifted_lons = data_lons - shifted_mlts + # Note that this "mlons" is adjusted for MLT + mlons = np.radians(shifted_lons) + mlats = data_lats + + # If the parameter is velocity then plot the LOS vectors + # Actual mlons used here, not adjusted mlons (np.radians(data_lons)) + if parameter == MapParams.FITTED_VELOCITY: + v_mag, azm_v =\ + cls.calculated_fitted_velocities(mlats=mlats, + mlons=np.radians( + data_lons), + hemisphere=hemisphere, + fit_coefficient=dmap_data[ + record]['N+2'], + fit_order=dmap_data[ + record]['fit.order'], + lat_min=dmap_data[ + record]['latmin'], + len_factor=len_factor) + elif parameter == MapParams.MODEL_VELOCITY: + v_mag = dmap_data[record]['model.vel.median'] + azm_v = np.radians(dmap_data[record]['model.kvect']) + elif parameter == MapParams.RAW_VELOCITY: + v_mag = dmap_data[record]['vector.vel.median'] + azm_v = np.radians(dmap_data[record]['vector.kvect']) + elif parameter == MapParams.POWER: + v_mag = dmap_data[record]['vector.pwr.median'] + azm_v = np.radians(dmap_data[record]['vector.kvect']) + + elif parameter == MapParams.SPECTRAL_WIDTH: + v_mag = dmap_data[record]['vector.wdt.median'] + azm_v = np.radians(dmap_data[record]['vector.kvect']) + + if parameter in [MapParams.FITTED_VELOCITY, MapParams.MODEL_VELOCITY, + MapParams.RAW_VELOCITY]: + # Angle to "rotate" each vector by to get into same + # reference frame Controlled by longitude, or "mltitude" + alpha = mlons + + # Convert initial positions to Cartesian + start_pos_x = (90 - abs(mlats)) * np.cos(mlons) + start_pos_y = (90 - abs(mlats)) * np.sin(mlons) + + # Resolve LOS vector in x and y directions, + # with respect to mag pole + # Gives zonal and meridional components of LOS vector + x = -v_mag * np.cos(-azm_v * hemisphere.value) + y = -v_mag * np.sin(-azm_v * hemisphere.value) + + # Rotate each vector into same reference frame + # following vector rotation matrix + # https://en.wikipedia.org/wiki/Rotation_matrix + vec_x = (x * np.cos(alpha)) - (y * np.sin(alpha)) + vec_y = (x * np.sin(alpha)) + (y * np.cos(alpha)) + + # New vector end points, in Cartesian + end_pos_x = start_pos_x + (vec_x * hemisphere.value / len_factor) + end_pos_y = start_pos_y + (vec_y * hemisphere.value / len_factor) + + # Convert back to polar for plotting + end_mlats = 90.0 - (np.sqrt(end_pos_x**2 + end_pos_y**2)) + end_mlons = np.arctan2(end_pos_y, end_pos_x) + + end_mlats=end_mlats * hemisphere.value + + # Plot the vectomlats + for i in range(len(v_mag)): + plt.plot([mlons[i], end_mlons[i]], + [mlats[i], end_mlats[i]], c=cmap(norm(v_mag[i])), + linewidth=0.5, zorder=5.0) + plt.scatter(mlons, mlats, c=v_mag, s=2.0, + vmin=zmin, vmax=zmax, cmap=cmap, zorder=5.0) + + # Plot potential contours + fit_coefficient = dmap_data[record]['N+2'] + fit_order = dmap_data[record]['fit.order'] + lat_shift = dmap_data[record]['lat.shft'] + lon_shift = dmap_data[record]['lon.shft'] + lat_min = dmap_data[record]['latmin'] + + cls.plot_potential_contours(fit_coefficient, lat_min, date, + lat_shift=lat_shift, lon_shift=lon_shift, + fit_order=fit_order, hemisphere=hemisphere, + **kwargs) + + if hmb is True: + # Plot the HMB + mlats_hmb = dmap_data[record]['boundary.mlat'] + mlons_hmb = dmap_data[record]['boundary.mlon'] + cls.plot_heppner_maynard_boundary(mlats_hmb, mlons_hmb, date) + + if colorbar is True: + mappable = cm.ScalarMappable(norm=norm, cmap=cmap) + locator = ticker.MaxNLocator(symmetric=True, min_n_ticks=3, + integer=True, nbins='auto') + ticks = locator.tick_values(vmin=zmin, vmax=zmax) + cb = ax.figure.colorbar(mappable, ax=ax, + extend='both', ticks=ticks) + + if colorbar_label != '': + cb.set_label(colorbar_label) + else: + if parameter in [MapParams.FITTED_VELOCITY, + MapParams.MODEL_VELOCITY, + MapParams.RAW_VELOCITY]: + cb.set_label('Velocity (m s$^{-1}$)') + elif parameter is MapParams.SPECTRAL_WIDTH: + cb.set_label('Spectral Width (m s$^{-1}$)') + elif parameter is MapParams.POWER: + cb.set_label('Power') + + if title == '': + title = "{year}-{month}-{day} {start_hour}:{start_minute} -"\ + " {end_hour}:{end_minute}"\ + "".format(year=date.year, + month=str(date.month).zfill(2), + day=str(date.day).zfill(2), + start_hour=str(date.hour).zfill(2), + start_minute=str(date.minute).zfill(2), + end_hour=str(dmap_data[record]['end.hour']). + zfill(2), + end_minute=str(dmap_data[record]['end.minute']). + zfill(2)) + plt.title(title) + + if map_info is True: + model = dmap_data[record]['model.name'] + num_points = len(dmap_data[record]['vector.mlat']) + pol_cap_pot = dmap_data[record]['pot.drop'] + cls.add_map_info(fit_order, pol_cap_pot, num_points, model) + + if imf_dial is True: + # Plot the IMF dial + bx = dmap_data[record]['IMF.Bx'] + by = dmap_data[record]['IMF.By'] + bz = dmap_data[record]['IMF.Bz'] + delay = dmap_data[record]['IMF.delay'] + bt = np.sqrt(bx**2 + by**2 + bz**2) + cls.plot_imf_dial(ax, by, bz, bt, delay) + + return mlats, mlons, v_mag + + + @classmethod + def index_legendre(cls, l: int, m: int): + """ + not a 100% how this works some black magic + + parameter + --------- + l : int + doping level + m : int + fit order + + return + ------ + legendre index? + """ + return (m == 0 and l**2) or ((l != 0) + and (m != 0) and l**2 + 2 * m - 1) or 0 + + + @classmethod + def calculated_fitted_velocities(cls, mlats: list, mlons: list, + fit_coefficient: list, + hemisphere: Enum = Hemisphere.North, + fit_order: int = 6, lat_min: int = 60, + len_factor: int = 150): + """ + Calculates the fitted velocities using Legrendre polynomial + + Parameters + ---------- + mlats: List[float] + Magnetic Latitude in degrees + mlons: List[float] + Magnetic Longitude in radians + fit_coefficient: List[float] + Value of the coefficient + hemisphere: int + 1 or -1 for hemisphere North or South + default: 1 - North + fit_order: int + order of the fit + default: 6 + lat_min: int + Lower latitude boundary of data in degrees + default: 60 + len_factor: int + length of the vector socks multiplied by + default: 150 + """ + # convert earth radius to meters + Re_meters = Re * 1000.0 + # theta values in radians + thetas = np.radians(90.0 - abs(mlats)) + thetas_max = np.radians(90.0 - abs(lat_min)) + + # Angle to "rotate" each vector by to get into same + # reference frame Controlled by longitude, or "mltitude" + alpha = np.pi / thetas_max + thetas_prime = alpha * thetas + x = np.cos(thetas_prime) + + # i is the index of the list + # x_i is the element of x at ith index + for i, x_i in enumerate(x): + temp_poly = special.lpmn(fit_order, fit_order, x_i) + if i == 0: + legendre_poly = np.append([temp_poly[0]], [temp_poly[0]], + axis=0) + else: + legendre_poly = np.append(legendre_poly, [temp_poly[0]], + axis=0) + legendre_poly = np.delete(legendre_poly, 0, 0) + phi = mlons + + # now do the index legender part, + # We are doing Associated Legendre Polynomials but + # for each polynomial we have two coefficients one + # for cos(phi) and the other for sin(phi), + # so we do spherical harmonics for a real valued function using + # sin(phi) and cos(phi) rather than exp(i*phi). + # we place an inner function to copying code + + # max index value + k_max = cls.index_legendre(fit_order, fit_order) + + # set up arrays and small stuff for the E field + # coefficients calculation + thetas_ecoeffs = np.zeros((k_max + 2, len(thetas))) + phi_ecoeffs = np.zeros((k_max + 2, len(thetas))) + + q_prime = np.array(np.where(thetas_prime != 0.0)) + q_prime = q_prime[0] + q = np.array(np.where(thetas != 0.0)) + q = q[0] + + # finally get to converting coefficients for the potential into + # coefficients for elec. Field + fit_coefficient_flat = fit_coefficient.flatten() + for m in range(fit_order + 1): + for l in range(m, fit_order + 1): + k3 = cls.index_legendre(l, m) + k4 = cls.index_legendre(l, m) + + if k3 >= 0: + thetas_ecoeffs[k4, q_prime] =\ + thetas_ecoeffs[k4, q_prime] -\ + fit_coefficient_flat[k3] * alpha * l *\ + np.cos(thetas_prime[q_prime]) / \ + np.sin(thetas_prime[q_prime]) / Re_meters + phi_ecoeffs[k4, q] = phi_ecoeffs[k4, q] - \ + fit_coefficient_flat[k3 + 1] * m /\ + np.sin(thetas[q]) / Re_meters + phi_ecoeffs[k4 + 1, q] = phi_ecoeffs[k4 + 1, q] + \ + fit_coefficient_flat[k3] * m /\ + np.sin(thetas[q]) / Re_meters + + if l < fit_order: + k1 = cls.index_legendre(l+1, m) + else: + k1 = -1 + + k2 = cls.index_legendre(l, m) + + if k1 >= 0: + thetas_ecoeffs[k2, q_prime] =\ + thetas_ecoeffs[k2, q_prime] + \ + fit_coefficient_flat[k1] * alpha * (l + 1 + m) / \ + np.sin(thetas_prime[q_prime]) / Re_meters + + if m > 0: + if k3 >= 0: + k3 = k3 + 1 + k4 = k4 + 1 + + if k1 >= 0: + k1 = k1 + 1 + k2 = k2 + 1 + + if k3 >= 0: + thetas_ecoeffs[k4, q_prime] =\ + thetas_ecoeffs[k4, q_prime] \ + - fit_coefficient_flat[k3] * alpha * l * \ + np.cos(thetas_prime[q_prime]) / \ + np.sin(thetas_prime[q_prime]) / Re_meters + + if k1 >= 0: + thetas_ecoeffs[k2, q_prime] = \ + thetas_ecoeffs[k2, q_prime] \ + + fit_coefficient_flat[k1] * alpha *\ + (l + 1 + m) / np.sin(thetas_prime[q_prime]) /\ + Re_meters + + # Calculate the Electric field positions + thetas_ecomp = np.zeros(thetas.shape) + phi_ecomp = np.zeros(thetas.shape) + + for m in range(fit_order + 1): + for l in range(m, fit_order + 1): + k = cls.index_legendre(l, m) + # Now in the IDL code we use + # legendre_poly[:,l,m] instead of + # legendre_poly[:,m,l] like here, this is + # because we have a different + # organization of legendre_poly due to the + # way scipy.special.lpmn + # stores values in arrays... + if m == 0: + thetas_ecomp = thetas_ecomp + thetas_ecoeffs[k, :] * \ + legendre_poly[:, m, l] + phi_ecomp = phi_ecomp + phi_ecoeffs[k, :] * \ + legendre_poly[:, m, l] + else: + thetas_ecomp = thetas_ecomp + thetas_ecoeffs[k, :] * \ + legendre_poly[:, m, l] * np.cos(m * phi) + \ + thetas_ecoeffs[k+1, :] * legendre_poly[:, m, l] * \ + np.sin(m * phi) + phi_ecomp = phi_ecomp + phi_ecoeffs[k, :] * \ + legendre_poly[:, m, l] * np.cos(m * phi) + \ + phi_ecoeffs[k+1, :] * legendre_poly[:, m, l] * \ + np.sin(m * phi) + + # Store the two components of Efield into a single array + E_field_fit = np.append([thetas_ecomp], [phi_ecomp], axis=0) + + # We'll calculate Bfield magnitude now, need to initialize some more + # stuff + # F-region altitude 300 km * 1000 to convert to meteres + F_altitude = 300.0 * 1000.0 + # dipole earth field in Tesla + B_field_polar = -0.62e-4 + B_field = B_field_polar * (1.0 - 3.0 * F_altitude / Re_meters) \ + * np.sqrt(3.0 * np.square(np.cos(thetas)) + 1.0) / 2 + + # get the velocity components from E-field + velocity_fit_vectors = np.zeros(E_field_fit.shape) + velocity_fit_vectors[0, :] = E_field_fit[1, :] / B_field + velocity_fit_vectors[1, :] = -E_field_fit[0, :] / B_field + velocity = np.sqrt(np.square(velocity_fit_vectors[0, :]) + + np.square(velocity_fit_vectors[1, :])) + velocity_chk_zero_inds = np.where(velocity != 0.0) + velocity_chk_zero_inds = velocity_chk_zero_inds[0] + + azm_v = np.zeros(velocity.shape) + + if len(velocity_chk_zero_inds) == 0: + velocity = np.array([0.0]) + azm_v = np.array([0.0]) + else: + if hemisphere == Hemisphere.South: + azm_v[velocity_chk_zero_inds] =\ + np.arctan2( + velocity_fit_vectors[1, velocity_chk_zero_inds], + velocity_fit_vectors[0, velocity_chk_zero_inds]) + else: + azm_v[velocity_chk_zero_inds] =\ + np.arctan2( + velocity_fit_vectors[1, velocity_chk_zero_inds], + -velocity_fit_vectors[0, velocity_chk_zero_inds]) + + return velocity, azm_v + + + @classmethod + def add_map_info(cls, fit_order: float, pol_cap_pot: float, + num_points: float, model: str): + """ + Annotates the plot with information about the map plotting + + Parameters + ---------- + ax: object + matplotlib axis object + fit_order: int + order of the fit + pol_cap_pot: float + value of the polar cap potential in kV + num_points: int + number of vectors plotted + model: str + model used to fit data + """ + text_string = r'$\phi_{PC}$' + ' = ' + str(round(pol_cap_pot/1000))\ + + ' kV\n' \ + + 'N = ' + str(num_points) + '\n' \ + + 'Order: ' + str(fit_order) + '\n' \ + + 'Model: ' + model + '\n' + plt.figtext(0.1, 0.1, text_string) + + + @classmethod + def plot_heppner_maynard_boundary(cls, mlats: list, mlons: list, + date: object, line_color: str = 'black', + **kwargs): + # TODO: No evaluation of coordinate system made! May need if in + # plotting to plot in radians/geo ect. + """ + Plots the position of the Heppner-Maynard Boundary + + Parameters + ---------- + ax: object + matplotlib axis object + mlats: List[float] + Magnetic Latitude in degrees + mlons: List[float] + Magnetic Longitude in radians + date: datetime object + Date from record + line_color: str + Color of the Heppner-Maynard boundary + Default: black + + """ + # Shift mlon to MLT + shifted_mlts = mlons[0] - \ + (aacgmv2.convert_mlt(mlons[0], date) * 15) + shifted_lons = mlons - shifted_mlts + mlon = np.radians(shifted_lons) + + plt.plot(mlon, mlats, c=line_color, zorder=4.0, **kwargs) + + + @classmethod + def plot_imf_dial(cls, ax: object, by: float = 0, bz: float = 0, + bt: float = 0, delay: float = 0): + """ + Plots an IMF clock angle dial on the existing plot + Defaults all to 0 if no IMF data available to plot + + Parameters + ---------- + ax: object + matplotlib axis object + by: Float + Value of the magnetic field in the y-direction (nT) + Default = 0 nT + bz: Float + Value of the magnetic field in the z-direction (nT) + Default = 0 nT + bt: Float + Magnitude of the magnetic field (nT) + Default = 0 nT + delay: Float + Time delay of magnetic field between the + measuring satellite and the ionosphere (minutes) + Default = 0 minutes + """ + # Create new axes inside existing axes + ax_imf = plt.axes([0, 0, 1, 1]) + ip = InsetPosition(ax, [-0.2, 0.7, 0.4, 0.4]) + ax_imf.set_axes_locator(ip) + ax_imf.axis('off') + + ax_imf.set_xlim([-20.2, 20.2]) + ax_imf.set_ylim([-20.2, 20.2]) + + # Plot a Circle + limit_circle = plt.Circle((0, 0), 10, facecolor='w', + edgecolor='k') + ax_imf.add_patch(limit_circle) + # Plot axis lines + plt.plot([-10, 10], [0, 0], color='k', linewidth=0.5) + plt.plot([0, 0], [-10, 10], color='k', linewidth=0.5) + + # Plot line for magnetic field + plt.plot([0, by], [0, bz], color='r') + + # Add axis labels + ax_imf.annotate('+Z', xy=(-2.5, 11)) + ax_imf.annotate('+Y', xy=(11, -1)) + + # Add annotations for delay and Btot + ax_imf.annotate('|B| = ' + str(round(bt)) + ' nT', xy=(-16, -13), + fontsize=7) + ax_imf.annotate('Delay = -' + str(delay) + ' min', xy=(-16, -17), + fontsize=7) + + + @classmethod + def calculate_potentials(cls, fit_coefficient: list, lat_min: list, + lat_shift: int = 0, lon_shift: int = 0, + fit_order: int = 6, lowlat: int = 60, + hemisphere: Enum = Hemisphere.North, + **kwargs): + # TODO: No evaluation of coordinate system made! May need if in + # plotting to plot in radians/geo ect. + ''' + Calculates potential across a magnetic lat/lon grid for + plotting later + + Parameters + ---------- + fit_coefficient: List[float] + Value of the coefficient + lat_min: List[float] + Minimum latitude that will be evaluated + Not to be confused with 'lowlat' + lat_shift: int + Generic shift in latitude from map file + default: 0 + lon_shift: int + Generic shift in longitude from map file + default: 0 + fit_order: int + order of the fit + default: 6 + lowlat: int + Lowest latitude on plot + default: 60 + hemisphere: Enum + Describes the hemisphere, North or South + default: Hemisphere.North + + ''' + # Lowest latitude to calculate potential to + theta_max = np.radians(90-np.abs(lat_min)) * hemisphere.value + + # Make a grid of the space the potential is evaluated on + # in magnetic coordinates + lat_step = 1 + lon_step = 2 + num_lats = int((90.0 - lowlat) / lat_step) + 1 + num_lons = int(360.0 / lon_step) + 1 + lat_arr = np.array(range(num_lats)) * lat_step + lowlat + lon_arr = np.array(range(num_lons)) * lon_step + + # Set up Grid + grid_arr = np.zeros((2, num_lats * num_lons)) + count1 = 0 + for lons in lon_arr: + for lats in lat_arr: + grid_arr[0, count1] = lats + grid_arr[1, count1] = lons + count1 += 1 + + # Convert grid vals to spherical coords + theta = np.radians(90.0 - np.abs(grid_arr[0, :])) + phi = np.radians(grid_arr[1, :]) + + # Adjusted/Normalised values (runs 0 - pi) + alpha = np.pi / theta_max + x = np.cos(alpha*theta) + + # Legendre Polys + for j, xj in enumerate(x): + plm_tmp = special.lpmn(fit_order, fit_order, xj) + if j == 0: + plm_fit = np.append([plm_tmp[0]], [plm_tmp[0]], axis=0) + else: + plm_fit = np.append(plm_fit, [plm_tmp[0]], axis=0) + # Remove first element as it is duplicated to start off the array + plm_fit = np.delete(plm_fit, 0, 0) + + # Eval the potential + lmax = plm_fit.shape + lmax = lmax[1] + v = np.zeros(phi.shape) + + coeff_fit_flat = fit_coefficient.flatten() + for m in range(lmax): + for l in range(m, lmax): + k = cls.index_legendre(l, m) + if m == 0: + v = v + coeff_fit_flat[k] * plm_fit[:, 0, l] + else: + v = v + coeff_fit_flat[k] * np.cos(m * phi) \ + * plm_fit[:, m, l] + coeff_fit_flat[k+1] \ + * np.sin(m * phi) * plm_fit[:, m, l] + + pot_arr = np.zeros((num_lons, num_lats)) + pot_arr = np.reshape(v, pot_arr.shape) / 1000.0 + + # TODO: Account for lon_shift + # TODO: Code for lat shift! (both rarely non-0 though) + # grid_arr[1,:] = (grid_arr[1,:] + lon_shift) + + mlat_center = grid_arr[0, :].reshape((num_lons, num_lats)) + # Set everything below the latmin as 0 + ind = np.where(abs(mlat_center) < abs(lat_min)) + pot_arr[ind] = 0 + + mlon_center = grid_arr[1, :].reshape((num_lons, num_lats)) + # Invert for Southern maps + mlat_center = mlat_center * hemisphere.value + + return mlat_center, mlon_center, pot_arr + + + @classmethod + def plot_potential_contours(cls, fit_coefficient: list, lat_min: list, + date: object, lat_shift: int = 0, + lon_shift: int = 0, fit_order: int = 6, + hemisphere: Enum = Hemisphere.North, + contour_levels: list = [], + contour_color: str = 'dimgrey', + contour_linewidths: float = 0.8, + contour_fill: bool = False, + contour_colorbar: bool = True, + contour_fill_cmap: str = 'RdBu', + contour_colorbar_label: str = 'Potential (kV)', + pot_minmax_color: str = 'k', **kwargs): + # TODO: No evaluation of coordinate system made! May need if in + # plotting to plot in radians/geo ect. + ''' + Takes the grid of potentials, plots a contour plot and min and max + potential positions + + Parameters + ---------- + fit_coefficient: List[float] + Value of the coefficient + lat_min: List[float] + Minimum latitude that will be evaluated + Not to be confused with 'lowlat' + date: datetime object + Date from record + lat_shift: int + Generic shift in latitude from map file + default: 0 + lon_shift: int + Generic shift in longitude from map file + default: 0 + fit_order: int + order of the fit + default: 6 + contour_levels: np.arr + Array of values at which the contours + are plotted + Default: [] + Default list is defined in function due + to length of the list, values higher or + lower than the minimum and maximum values + given are colored in as min and max color + values if contour_fill=True + contour_color: str + Colour of the contour lines plotted + Default: dimgrey + contour_label: bool - NOT CURRENTLY IMPLEMENTED + If contour_fill is True, contour labels will + be plotted on the contour lines + Default: True + contour_linewidths: float + Thickness of contour lines + Default: 0.8 + contour_fill: bool + Option to use filled contours rather than + an outline. If True, contour_color and + contour_linewidths are ignored + If False + Default: False + contour_colorbar: bool + Option to show the colorbar for the contours + if contour_fill = True + Default: True + contour_fill_cmap: matplotlib.cm + Colormap used to fill the contours if + contour_fill is True + Default: 'RdBu' + contour_colorbar_label: str + Label for the colorbar describing the + contours if contour_fill is True + Default: empty string '' + pot_minmax_color: str + Colour of the cross and plus symbols for + minimum and maximum potentials + Default: 'k' - black + **kwargs + including lowlat and hemisphere for calculating + potentials + ''' + mlat, mlon_u, pot_arr = cls.calculate_potentials( + fit_coefficient, lat_min, + lat_shift=lat_shift, lon_shift=lon_shift, + fit_order=fit_order, hemisphere=hemisphere, + **kwargs) + + # Shift mlon to MLT + shifted_mlts = mlon_u[0, 0] - \ + (aacgmv2.convert_mlt(mlon_u[0, 0], date) * 15) + shifted_lons = mlon_u - shifted_mlts + mlon = shifted_lons + + # Contained in function as too long to go into the function call + if contour_levels == []: + contour_levels = [-100, -95, -90, -85, -80, -75, -70, -65, -60, + -55, -50, -45, -40, -35, -30, -25, -20, -15, + -10, -5, -1, 1, 5, 10, 15, 20, 25, 30, 35, 40, + 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100] + + if contour_fill: + # Filled contours + plt.contourf(np.radians(mlon), mlat, pot_arr, 2, + vmax=abs(pot_arr).max(), + vmin=-abs(pot_arr).max(), + locator=ticker.FixedLocator(contour_levels), + cmap=contour_fill_cmap, alpha=0.6, + extend='both', zorder=3.0) + if contour_colorbar is True: + norm = colors.Normalize + norm = norm(-abs(pot_arr).max(), abs(pot_arr).max()) + mappable = cm.ScalarMappable(norm=norm, cmap=contour_fill_cmap) + locator = ticker.MaxNLocator(symmetric=True, min_n_ticks=3, + integer=True, nbins='auto') + ticks = locator.tick_values(vmin=-abs(pot_arr).max(), + vmax=abs(pot_arr).max()) + cb = plt.colorbar(mappable, extend='both', ticks=ticks) + if contour_colorbar_label != '': + cb.set_label(contour_colorbar_label) + else: + # Contour lines only + cs = plt.contour(np.radians(mlon), mlat, pot_arr, 2, + vmax=abs(pot_arr).max(), + vmin=-abs(pot_arr).max(), + locator=ticker.FixedLocator(contour_levels), + colors=contour_color, alpha=0.8, + linewidths=contour_linewidths, zorder=3.0) + # TODO: Add in contour labels + # if contour_label: + # plt.clabel(cs, cs.levels, inline=True, fmt='%d', fontsize=5) + + # Get max value of potential + ind_max = np.where(pot_arr == pot_arr.max()) + ind_min = np.where(pot_arr == pot_arr.min()) + max_mlon = mlon[ind_max] + max_mlat = mlat[ind_max] + min_mlon = mlon[ind_min] + min_mlat = mlat[ind_min] + + plt.scatter(np.radians(max_mlon), max_mlat, marker='+', s=70, + color=pot_minmax_color, zorder=5.0) + plt.scatter(np.radians(min_mlon), min_mlat, marker='_', s=70, + color=pot_minmax_color, zorder=5.0) diff --git a/pydarn/plotting/projections.py b/pydarn/plotting/projections.py new file mode 100644 index 00000000..ac65c0df --- /dev/null +++ b/pydarn/plotting/projections.py @@ -0,0 +1,128 @@ +# Copyright (C) 2021 SuperDARN Canada, University of Saskatchewan +# Author: Daniel Billett +# +# Modifications: +# 2022-03-11 MTS added enums for projection function calls +# 2022-03-22 MTS removed coastline call and added grid lines to cartopy plotting +# Disclaimer: +# pyDARN is under the LGPL v3 license found in the root directory LICENSE.md +# Everyone is permitted to copy and distribute verbatim copies of this license +# document, but changing it is not allowed. +# +# This version of the GNU Lesser General Public License incorporates the terms +# and conditions of version 3 of the GNU General Public License, +# supplemented by the additional permissions listed below. +""" +Code which generates axis objects for use in plotting functions +""" +import enum +import matplotlib.pyplot as plt +import numpy as np +from packaging import version + +from pydarn import Hemisphere, plot_exceptions +try: + import cartopy + # from cartopy.mpl import geoaxes + import cartopy.crs as ccrs + cartopyInstalled = True + if version.parse(cartopy.__version__) < version.parse("0.19"): + cartopyVersion = False + else: + cartopyVersion = True +except Exception: + cartopyInstalled = False + + +def axis_polar(lowlat: int = 30, hemisphere: Hemisphere = Hemisphere.North, + grid_lines: bool = True, **kwargs): + """ + Plots a radar's Field Of View (FOV) fan plot for the given data and + scan number + + Parameters + ----------- + lowlat: int + Lower AACGM latitude boundary for the polar plot + Default: 30 + hemiphere: enum + Hemisphere of polar projection. Can be Hemisphere.North or + Hemisphere.South for northern and southern hemispheres, + respectively + Default: Hemisphere.North + """ + + ax = plt.axes(polar=True) + + # Set upper and lower latitude limits (pole and lowlat) + if hemisphere == Hemisphere.North: + ax.set_ylim(90, lowlat) + ax.set_yticks(np.arange(lowlat, 90, 10)) + else: + # If hemisphere is South, lowlat must be negative + ax.set_ylim(-90, -abs(lowlat)) + ax.set_yticks(np.arange(-abs(lowlat), -90, -10)) + + # Locations of tick marks. Will be customisable in future + ax.set_xticks([0, np.radians(45), np.radians(90), np.radians(135), + np.radians(180), np.radians(225), np.radians(270), + np.radians(315)]) + + # Tick labels will depend on coordinate system + ax.set_xticklabels(['00', '', '06', '', '12', '', '18', '']) + ax.set_theta_zero_location("S") + + return ax, None + + +def axis_geological(date, ax: object = None, + hemisphere: Hemisphere = Hemisphere.North, + lowlat: int = 30, grid_lines: bool = True, **kwargs): + """ + """ + if cartopyInstalled is False: + raise plot_exceptions.CartopyMissingError() + if cartopyVersion is False: + raise plot_exceptions.CartopyVersionError(cartopy.__version__) + # no need to shift any coords, let cartopy do that + # however, we do need to figure out + # how much to rotate the projection + deg_from_midnight = (date.hour + date.minute / 60) / 24 * 360 + if hemisphere == Hemisphere.North: + pole_lat = 90 + noon = -deg_from_midnight + ylocations = -5 + else: + pole_lat = -90 + noon = 360 - deg_from_midnight + ylocations = 5 + # handle none types or wrongly built axes + proj = ccrs.Orthographic(noon, pole_lat) + ax = plt.subplot(111, projection=proj, aspect='auto') + if grid_lines: + ax.gridlines(draw_labels=True) + + extent = min(45e5, + (abs(proj.transform_point(noon, lowlat, + ccrs.PlateCarree()) + [1]))) + ax.set_extent(extents=(-extent, extent, -extent, extent), + crs=proj) + return ax, ccrs + + +class Projs(enum.Enum): + """ + class of projections that pydarn can do + + Enumerators + ------------ + POLAR: axes_polar + GEO: axes_geological + """ + POLAR = (axis_polar, ) + GEO = (axis_geological, ) + + # Need this to make the functions callable + def __call__(self, *args, **kwargs): + return self.value[0](*args, **kwargs) diff --git a/pydarn/plotting/rtp.py b/pydarn/plotting/rtp.py index c999ef05..7df83f4c 100644 --- a/pydarn/plotting/rtp.py +++ b/pydarn/plotting/rtp.py @@ -7,6 +7,8 @@ # 2021-05-12 Francis Tholley added gate2grounscatter to range-time plots # 2021-06-18 Marina Schmidt (SuperDARN Canada) fixed ground scatter colour bug # 2021-07-06 Carley Martin added keyword to aid in rounding start times +# 2022-03-04 Marina Schmidt added RangeEstimations in +# # Disclaimer: # pyDARN is under the LGPL v3 license found in the root directory LICENSE.md # Everyone is permitted to copy and distribute verbatim copies of this license @@ -30,11 +32,10 @@ from matplotlib import dates, colors, cm, ticker from typing import List -from pydarn import (gate2GroundScatter, gate2slant, check_data_type, +from pydarn import (RangeEstimation, check_data_type, time2datetime, rtp_exceptions, plot_exceptions, SuperDARNCpids, SuperDARNRadars, - standard_warning_format, PyDARNColormaps, - Coords) + standard_warning_format, PyDARNColormaps) warnings.formatwarning = standard_warning_format @@ -73,12 +74,13 @@ def plot_range_time(cls, dmap_data: List[dict], parameter: str = 'v', start_time: datetime = None, end_time: datetime = None, colorbar: plt.colorbar = None, ymin: int = None, ymax: int = None, yspacing: int = 200, - coords: object = Coords.SLANT_RANGE, + range_estimation: RangeEstimation = + RangeEstimation.SLANT_RANGE, colorbar_label: str = '', norm=colors.Normalize, cmap: str = None, filter_settings: dict = {}, date_fmt: str = '%y/%m/%d\n %H:%M', - round_start: bool=True, **kwargs): + round_start: bool = True, **kwargs): """ Plots a range-time parameter plot of the given field name in the dmap_data @@ -127,9 +129,9 @@ def plot_range_time(cls, dmap_data: List[dict], parameter: str = 'v', yspacing: int sets the spacing between ticks Default: 200 - coords: Coords - set the y-axis to a desired coordinate system - Default: Coords.SLANT_RANGE + range_estimation: RangeEstimation + set the y-axis to a desired range estimation calculation + Default: RangeEstimation.SLANT_RANGE norm: matplotlib.colors.Normalization object This object use dependency injection to use any normalization method with the zmin and zmax. @@ -356,18 +358,19 @@ def plot_range_time(cls, dmap_data: List[dict], parameter: str = 'v', start_time=start_time, end_time=end_time, opt_beam_num=cls.dmap_data[0]['bmnum']) - if coords is Coords.SLANT_RANGE or\ - coords is Coords.GROUND_SCATTER_MAPPED_RANGE: + if range_estimation != RangeEstimation.RANGE_GATE: # Get rxrise from hardware files (consistent with RST) rxrise = SuperDARNRadars.radars[cls.dmap_data[0]['stid']]\ .hardware_info.rx_rise_time - y = gate2slant(cls.dmap_data[0]['frang'], cls.dmap_data[0]['rsep'], - rxrise, nrang=y_max) - if coords is Coords.GROUND_SCATTER_MAPPED_RANGE: - y = gate2GroundScatter(y, **kwargs) - y0inx = np.min(np.where(np.isfinite(y))[0]) - y = y[y0inx:] - z = z[:, y0inx:] + frang = int(cls.dmap_data[0]['frang']) + rsep = int(cls.dmap_data[0]['rsep']) + + y = range_estimation(frang=frang, rxrise=rxrise, + rsep=rsep, nrang=y_max, **kwargs) + + y0inx = np.min(np.where(np.isfinite(y))[0]) + y = y[y0inx:] + z = z[:, y0inx:] time_axis, y_axis = np.meshgrid(x, y) z_data = np.ma.masked_where(np.isnan(z.T), z.T) @@ -429,8 +432,7 @@ def plot_range_time(cls, dmap_data: List[dict], parameter: str = 'v', ax.set_ylim(ymin, ymax) - if coords is Coords.SLANT_RANGE or\ - coords is Coords.GROUND_SCATTER_MAPPED_RANGE: + if range_estimation != RangeEstimation.RANGE_GATE: ax.yaxis.set_ticks(np.arange(np.ceil(ymin/100.0)*100, ymax+1, yspacing)) else: @@ -449,7 +451,7 @@ def plot_range_time(cls, dmap_data: List[dict], parameter: str = 'v', # byminute keyword makes sure that the ticks are situated at # the minute or half hour marks, rather than at a set interval ax.xaxis.set_minor_locator( - dates.MinuteLocator(byminute=range(0,60,tick_interval))) + dates.MinuteLocator(byminute=range(0, 60, tick_interval))) # Upon request of Daniel Billet and others, I am rounding # the time down so the plotting x-axis will show the origin @@ -464,21 +466,20 @@ def plot_range_time(cls, dmap_data: List[dict], parameter: str = 'v', if tick_sep > 0: rounded_down_start_time = x[0] -\ timedelta(minutes=x[0].minute % tick_sep, - seconds=x[0].second, - microseconds=x[0].microsecond) + seconds=x[0].second, + microseconds=x[0].microsecond) else: rounded_down_start_time = x[0] -\ timedelta(minutes=x[0].minute % 15, - seconds=x[0].second, - microseconds=x[0].microsecond) + seconds=x[0].second, + microseconds=x[0].microsecond) else: rounded_down_start_time = x[0] ax.set_xlim([rounded_down_start_time, x[-1]]) ax.xaxis.set_major_formatter(dates.DateFormatter(date_fmt)) - if coords is Coords.SLANT_RANGE or\ - coords is Coords.GROUND_SCATTER_MAPPED_RANGE: + if range_estimation != RangeEstimation.RANGE_GATE: ax.yaxis.set_minor_locator(ticker.AutoMinorLocator(2)) else: ax.yaxis.set_minor_locator(ticker.MultipleLocator(5)) @@ -518,7 +519,7 @@ def plot_time_series(cls, dmap_data: List[dict], channel='all', scale: str = 'linear', cp_name: bool = True, color: str = 'black', linestyle: str = '-', linewidth: float = 1, - round_start: bool=True, **kwargs): + round_start: bool = True, **kwargs): """ Plots the time series of a scalar parameter @@ -739,13 +740,13 @@ def plot_time_series(cls, dmap_data: List[dict], if tick_sep > 0: rounded_down_start_time = x[0] -\ timedelta(minutes=x[0].minute % tick_sep, - seconds=x[0].second, - microseconds=x[0].microsecond) + seconds=x[0].second, + microseconds=x[0].microsecond) else: rounded_down_start_time = x[0] -\ timedelta(minutes=x[0].minute % 15, - seconds=x[0].second, - microseconds=x[0].microsecond) + seconds=x[0].second, + microseconds=x[0].microsecond) else: rounded_down_start_time = x[0] @@ -764,13 +765,13 @@ def plot_time_series(cls, dmap_data: List[dict], if tick_sep > 0: rounded_down_start_time = x[0] -\ timedelta(minutes=x[0].minute % tick_sep, - seconds=x[0].second, - microseconds=x[0].microsecond) + seconds=x[0].second, + microseconds=x[0].microsecond) else: rounded_down_start_time = x[0] -\ timedelta(minutes=x[0].minute % 15, - seconds=x[0].second, - microseconds=x[0].microsecond) + seconds=x[0].second, + microseconds=x[0].microsecond) else: rounded_down_start_time = x[0] @@ -791,7 +792,7 @@ def plot_time_series(cls, dmap_data: List[dict], # byminute keyword makes sure that the ticks are situated at # the minute or half hour marks, rather than at a set interval ax.xaxis.set_minor_locator( - dates.MinuteLocator(byminute=range(0,60,tick_interval))) + dates.MinuteLocator(byminute=range(0, 60, tick_interval))) ax.margins(x=0) ax.tick_params(axis='y', which='minor') @@ -806,7 +807,8 @@ def plot_summary(cls, dmap_data: List[dict], plot_elv: bool = True, title=None, background: str = 'w', groundscatter: bool = True, channel: int = 'all', line_color: dict = {}, - coords: object = Coords.SLANT_RANGE, **kwargs): + range_estimation: object = + RangeEstimation.SLANT_RANGE, **kwargs): """ Plots the summary of several SuperDARN parameters using time-series and range-time plots. Please see Notes for further description @@ -838,9 +840,10 @@ def plot_summary(cls, dmap_data: List[dict], channel number that will be plotted in the summary plot. Default: 'all' - coords: Coord object - set the y-axis to a desired coordinate system - Default: Coord.SLANT_RANGE + range_estimation: RangeEstimation object + set the y-axis to a desired range gate estimation + calculation + Default: RangeEstimation.SLANT_RANGE figsize : (int,int) tuple containing (height, width) figure size Default: 11 x 8.5 @@ -1022,9 +1025,11 @@ def plot_summary(cls, dmap_data: List[dict], cls.plot_time_series(dmap_data, beam_num=beam_num, parameter=axes_parameters[i][0], scale=scale, channel=channel, - color=color[axes_parameters[i][0]], + color=color[ + axes_parameters[i][0]], ax=axes[i], - linestyle=line[axes_parameters[i][0]], + linestyle=line[ + axes_parameters[i][0]], label=labels[i][0], **kwargs) axes[i].set_ylabel(labels[i][0], rotation=0, labelpad=30) axes[i].\ @@ -1056,11 +1061,15 @@ def plot_summary(cls, dmap_data: List[dict], with warnings.catch_warnings(): warnings.simplefilter("ignore") cls.plot_time_series(dmap_data, beam_num=beam_num, - parameter=axes_parameters[i][1], - color=color[axes_parameters[i][1]], + parameter=axes_parameters[ + i][1], + color=color[axes_parameters[ + i][1]], channel=channel, scale=scale, ax=second_ax, - linestyle=line[axes_parameters[i][1]], **kwargs) + linestyle=line[ + axes_parameters[i][1]], + **kwargs) second_ax.set_xticklabels([]) second_ax.set_ylabel(labels[i][1], rotation=0, labelpad=25) @@ -1104,9 +1113,9 @@ def plot_summary(cls, dmap_data: List[dict], else: # Current standard is to only have groundscatter # on the velocity plot. This may change in the future. - if coords is Coords.SLANT_RANGE: + if range_estimation == RangeEstimation.SLANT_RANGE: ymax = 3517.5 - elif coords is Coords.GROUND_SCATTER_MAPPED_RANGE: + elif range_estimation == RangeEstimation.GSMR: ymax = 3517.5/2 else: ymax = 75 @@ -1126,11 +1135,14 @@ def plot_summary(cls, dmap_data: List[dict], parameter=axes_parameters[i], ax=axes[i], groundscatter=grndflg, cmap=cmap[axes_parameters[i]], - zmin=boundary_ranges[axes_parameters[i]][0], - zmax=boundary_ranges[axes_parameters[i]][1], + zmin=boundary_ranges[ + axes_parameters[i]][0], + zmax=boundary_ranges[ + axes_parameters[i]][1], ymax=ymax, yspacing=500, background=background, - coords=coords, **kwargs) + range_estimation=range_estimation, + **kwargs) # Overwriting velocity ticks to get a better pleasing # look on the colorbar # Preference by Marina Schmidt @@ -1138,17 +1150,19 @@ def plot_summary(cls, dmap_data: List[dict], locator = ticker.LinearLocator(numticks=5) ticks =\ locator.\ - tick_values(vmin=boundary_ranges[axes_parameters[i]][0], - vmax=boundary_ranges[axes_parameters[i]][1]) + tick_values(vmin=boundary_ranges[ + axes_parameters[i]][0], + vmax=boundary_ranges[ + axes_parameters[i]][1]) if ticks[0] < boundary_ranges[axes_parameters[i]][0]: ticks[0] = boundary_ranges[axes_parameters[i]][0] if ticks[-1] > boundary_ranges[axes_parameters[i]][1]: ticks[-1] = boundary_ranges[axes_parameters[i]][1] cbar.set_ticks(ticks) - if coords is Coords.SLANT_RANGE: + if range_estimation == RangeEstimation.SLANT_RANGE: axes[i].set_ylabel('Slant Range (km)') - elif coords is Coords.GROUND_SCATTER_MAPPED_RANGE: + elif range_estimation == RangeEstimation.GSMR: axes[i].set_ylabel('Ground Scatter\nMapped Range\n(km)') else: axes[i].set_ylabel('Range Gates') diff --git a/pydarn/utils/conversions.py b/pydarn/utils/conversions.py index 358e5736..c8d03013 100644 --- a/pydarn/utils/conversions.py +++ b/pydarn/utils/conversions.py @@ -1,7 +1,18 @@ # Copyright 2019 SuperDARN Canada, University of Saskatchewan # Author: Marina Schmidt -# Copyright 2021 University of Scranton -# Author: Francis Tholley and Dr. Nathaniel Frissel for gate2groundscatter +# +# Disclaimer: +# pyDARN is under the LGPL v3 license found in the root directory LICENSE.md +# Everyone is permitted to copy and distribute verbatim copies of this license +# document, but changing it is not allowed. +# +# This version of the GNU Lesser General Public License incorporates the terms +# and conditions of version 3 of the GNU General Public License, +# supplemented by the additional permissions listed below. +# +# Modification: +# 2022-03-10 MTS removed gate2slant and gate2groundscatter +# to range_estimations.py """ This module is to focus on data conversions @@ -12,8 +23,6 @@ import numpy as np from collections import OrderedDict -from pydarn import (Re, C) - # key is the format char type defined by python, # item is the DMAP int value for the type @@ -40,8 +49,8 @@ np.float32: 'f', np.float64: 'd', str: 's', - np.str_: 's', - np.str: 's', + str: 's', + str: 's', np.int64: 'q', np.uint8: 'B', np.uint16: 'H', @@ -84,92 +93,3 @@ def dmap2dict(dmap_records: List[dict]) -> List[dict]: for field, data in dmap_record.items()} dmap_list.append(OrderedDict(dmap_dict)) return dmap_list - - -def gate2GroundScatter(slant_ranges: List[float], - reflection_height: float = 250): - """ - Calculate the ground scatter mapped range (km) for each slanted range - for SuperDARN data. This function is based on the Ground Scatter equation - from Bristow paper at https://doi.org/10.1029/93JA01470 on page 325 - Parameters - ---------- - slant_ranges : List[float] - list of slant ranges - reflection_height: float - reflection height - default: 250 - - Returns - ------- - ground_scatter_mapped_ranges : np.array - returns an array of ground scatter mapped ranges for the radar - """ - - ground_scatter_mapped_ranges =\ - Re*np.arcsin(np.sqrt((slant_ranges**2/4)-(reflection_height**2))/Re) - - return ground_scatter_mapped_ranges - - -def gate2slant(frang:int, rsep:int, rxrise:int, gate: int = 0, - nrang: int = None, center: bool = True): - """ - Calculate the slant range (km) for each range gate for SuperDARN data - - Parameters - ---------- - frang: int - range from the edge of first the gate to the radar [km] - This should be given in fitacf record of the control program - rsep: int - Radar seperation of the gates. Determined by control program. - rxrise: int - Use hardware value for this, avoid data file values - gate: int - range gate to determine the slant range [km], if nrang - is None - default: 0 - nrang: int - max number of range gates in the list of records. If - not None, will calculate all slant ranges - default: None - center: boolean - Calculate the slant range in the center of range gate - or edge - - Returns - ------- - slant_ranges : np.array - returns an array of slant ranges for the radar - """ - - # lag to the first range gate in microseconds - # 0.3 - speed of light (km/us) - # 2 - two times for there and back - distance_factor = 2.0 - # C - speed of light m/s to km/us - speed_of_light = C * 0.001 * 1e-6 - lag_first = frang * distance_factor / speed_of_light - - # sample separation in microseconds - sample_sep = rsep * distance_factor / speed_of_light - # Range offset - # If center is true, calculate at the center - if center: - # 0.5 off set to the centre of the range gate instead of edge - range_offset = -0.5 * rsep - else: - range_offset = 0.0 - # Now calculate slant range in km - if nrang is None: - slant_ranges = (lag_first - rxrise + - gate * sample_sep) * speed_of_light /\ - distance_factor + range_offset - else: - slant_ranges = np.zeros(nrang+1) - for gate in range(nrang+1): - slant_ranges[gate] = (lag_first - rxrise + - gate * sample_sep) * speed_of_light /\ - distance_factor + range_offset - return slant_ranges diff --git a/pydarn/utils/coordinates.py b/pydarn/utils/coordinates.py index 83a734f9..363b3b5e 100644 --- a/pydarn/utils/coordinates.py +++ b/pydarn/utils/coordinates.py @@ -1,24 +1,202 @@ # (C) Copyright 2021 SuperDARN Canada, University of Saskatchewan # Author(s): Marina Schmidt # -# This file is part of the pyDARN Library. -# +# Disclaimer: # pyDARN is under the LGPL v3 license found in the root directory LICENSE.md # Everyone is permitted to copy and distribute verbatim copies of this license # document, but changing it is not allowed. # # This version of the GNU Lesser General Public License incorporates the terms -# and conditions of version 3 of the GNU Lesser General Public License, +# and conditions of version 3 of the GNU General Public License, # supplemented by the additional permissions listed below. # # Modifications: +# 2022-03-10 MTS added 4 new methods to generate coordinates for the various +# enums # """ coordinates.py is a group of methods that focus on coordinates systems used in plotting """ +import datetime as dt import enum +import numpy as np + +import aacgmv2 +from pydarn import (geocentric_coordinates, SuperDARNRadars, RangeEstimation, + radar_exceptions, Re) + + +def geo_coordinates(stid: int, beams: int = None, + gates: tuple = None, **kwargs): + """ + geographic_coordinates calculates the geographic coordinate for a given + set of gates and beams + parameters + ----------- + + """ + if gates is None: + gates = [0, SuperDARNRadars.radars[stid].range_gate_45] + if beams is None: + beams = SuperDARNRadars.radars[stid].hardware_info.beams + + # Plus 1 is due to the fact fov files index at 1 so in the plotting + # of the boundary there is a subtraction of 1 to offset this as python + # converts to index of 0 which my code already accounts for + beam_corners_lats = np.zeros((gates[1]-gates[0]+1, beams+1)) + beam_corners_lons = np.zeros((gates[1]-gates[0]+1, beams+1)) + + for beam in range(0, beams+1): + for gate in range(gates[0], gates[1]+1): + lat, lon = gate2geographic_location(stid=stid, beam=beam, + range_gate=gate, height=300, + **kwargs) + beam_corners_lats[gate-gates[0], beam] = lat + beam_corners_lons[gate-gates[0], beam] = lon + y0inx = np.min(np.where(np.isfinite(beam_corners_lats[:,0]))[0]) + return beam_corners_lats[y0inx:], beam_corners_lons[y0inx:] + + +def aacgm_coordinates(stid: int, beams: int = None, gates: tuple = None, + date: dt.datetime = dt.datetime.now, **kwargs): + if gates is None: + gates = [0, SuperDARNRadars.radars[stid].range_gate_45] + if beams is None: + beams = SuperDARNRadars.radars[stid].hardware_info.beams + + # Plus 1 is due to the fact fov files index at 1 so in the plotting + # of the boundary there is a subtraction of 1 to offset this as python + # converts to index of 0 which my code already accounts for + beam_corners_lats = np.zeros((gates[1]-gates[0]+1, beams+1)) + beam_corners_lons = np.zeros((gates[1]-gates[0]+1, beams+1)) + + for beam in range(0, beams+1): + for gate in range(gates[0], gates[1]+1): + lat, lon = gate2geographic_location(stid=stid, beam=beam, + range_gate=gate, height=300, + **kwargs) + geomag = np.array(aacgmv2.get_aacgm_coord(glat=lat, + glon=lon, + height=250, + dtime=date)) + beam_corners_lats[gate-gates[0], beam] = geomag[0] + beam_corners_lons[gate-gates[0], beam] = geomag[1] + y0inx = np.min(np.where(np.isfinite(beam_corners_lats[:,0]))[0]) + return beam_corners_lats[y0inx:], beam_corners_lons[y0inx:] + + +def aacgm_MLT_coordinates(**kwargs): + beam_corners_lats, beam_corners_lons = aacgm_coordinates(**kwargs) + beam_corners_mlts = convert2MLT(beam_corners_lons, **kwargs) + return beam_corners_lats, beam_corners_mlts + + +def convert2MLT(lons: float, date: object, **kwargs): + fan_shape = lons.shape + # Work out shift due in MLT + beam_corners_mlts = np.zeros((fan_shape[0], fan_shape[1])) + mltshift = lons[0, 0] - (aacgmv2.convert_mlt(lons[0, 0], date) * 15) + beam_corners_mlts = lons - mltshift + return beam_corners_mlts + + +# RPosGeo line 335 +def gate2geographic_location(stid: int, beam: int, height: float = None, + elv_angle: float = 0.0, center: bool = True, + range_estimation: RangeEstimation = + RangeEstimation.SLANT_RANGE, **kwargs): + """ + determines the geographic cell position for a given range gate and beam + + parameters + ---------- + stid: int + station id of the radar to use + beam: int + beam number (indexing at 0) + range_gate: int + range gate number (indexing at 0) + rsep: int + range seperation, determined by the mode the + radar is using, in [km] + default: 45 - normalscan mode + frang: int + frequency range from the radar to the front edge of the range gate + Please note: this definition may be changed, currently defined in + RST code to keep consistency + height: float + transmutation height [km] + default: none + if none then it uses elevation angle + elv_angle: float + elevation angle in [deg] + default: 0 + center: bool + obtain geographic location of the centre of the range gates + False obtains the front edge of the range gates. + default: True + + returns + ------- + lat: float + latitude of the range gate in geographic coordinates [deg] + lon: float + longitude of the range gate in geographic coordinates [deg] + """ + # centre of the field of view + offset = SuperDARNRadars.radars[stid].hardware_info.beams / 2.0 - 0.5 + + # Obtain radar information from hardware files all converted to [rad] + # radians are needed for numpy geometry calculations and to reduce + # too much converting back and forth. + boresight = np.radians(SuperDARNRadars. + radars[stid].hardware_info.boresight.physical) + radar_lat = np.radians(SuperDARNRadars. + radars[stid].hardware_info.geographic.lat) + radar_lon = np.radians(SuperDARNRadars. + radars[stid].hardware_info.geographic.lon) + # Some beam seperations are negative which changes how the coordinates wrap + # we absolute to make it easier of fov-color for cartopy plotting + beam_sep = np.radians(abs(SuperDARNRadars. + radars[stid].hardware_info.beam_separation)) + rxrise = SuperDARNRadars.radars[stid].hardware_info.rx_rise_time + + # TODO: fix after slant range change + if center is True: + # beam edge in [rad] + beam_edge = -beam_sep * 0.5 + # range_edge in [km] + else: + beam_edge = 0 + + # psi [rad] in the angle from the boresight + psi = beam_sep * (beam - offset) + beam_edge + # Calculate the slant range [km] + if range_estimation != RangeEstimation.RANGE_GATE: + target_range = range_estimation(rxrise=rxrise, **kwargs) + else: + raise radar_exceptions.RangeEstimationError("Range Gates cannot be" + "used in estimating the" + " km for geographic" + " coordinates systems") + + # If no height is specified then use elevation angle (default 0) + # to calculate the transmutation height + if height is None: + height = -Re + np.sqrt(Re**2 + 2 * target_range * Re * + np.sin(np.radians(elv_angle)) + target_range**2) + + lat, lon = geocentric_coordinates(lat=radar_lat, + lon=radar_lon, + target_range=target_range, + cell_height=height, + psi=psi, + boresight=boresight, **kwargs) + # convert back degrees as preferred units to use? + return np.degrees(lat), np.degrees(lon) + class Coords(enum.Enum): """ @@ -33,8 +211,10 @@ class Coords(enum.Enum): AACGM: Magnetic geographical coordinates lat and longitude (degree) """ - RANGE_GATE = enum.auto() - SLANT_RANGE = enum.auto() - GROUND_SCATTER_MAPPED_RANGE = enum.auto() - GEOGRAPHIC = enum.auto() - AACGM = enum.auto() + GEOGRAPHIC = (geo_coordinates, ) + AACGM = (aacgm_coordinates, ) + AACGM_MLT = (aacgm_MLT_coordinates, ) + + # Need this to make the functions callable + def __call__(self, *args, **kwargs): + return self.value[0](*args, **kwargs) diff --git a/pydarn/utils/radar_pos.py b/pydarn/utils/geo.py similarity index 52% rename from pydarn/utils/radar_pos.py rename to pydarn/utils/geo.py index 3e1b8253..f310a23a 100644 --- a/pydarn/utils/radar_pos.py +++ b/pydarn/utils/geo.py @@ -16,10 +16,8 @@ # and converted to python for usage in pyDARN # # Modifications: -# 2020-04-20 Marina Schmidt converted the above link to python and changed -# variable and function names to readability -# 2020-09-15 Marina Schmidt removed fov file reading option -# 2021-09-15 Francis Tholley relocated the virtual height models to another file +# 2022-03-10 MTS added kwargs to the functions just to avoid errors when +# extra keys are passed in from other functions # Disclaimer: # pyDARN is under the LGPL v3 license found in the root directory LICENSE.md # Everyone is permitted to copy and distribute verbatim copies of this license @@ -29,178 +27,19 @@ # and conditions of version 3 of the GNU General Public License, # supplemented by the additional permissions listed below. # -# - - """ -This module is used for handling coordinates of a specified radar -in AACGMv2 or geographic coordinates +This module focuses on geographic calculations used for converting gate +information to geographic location """ - -import datetime as dt import numpy as np -import os - -import aacgmv2 - -from pydarn import SuperDARNRadars, gate2slant, Coords -from pydarn.utils.constants import EARTH_EQUATORIAL_RADIUS, Re, C - -from pydarn.utils.virtual_heights_types import VH_types -from pydarn.utils.virtual_heights import standard_virtual_height, chisham - -def radar_fov(stid: int, rsep: int = 45, frang: int = 180, - ranges: tuple = None, coords: object = Coords.AACGM, - max_beams: int = None, date: dt.datetime = None, **kwargs): - """ - Returning beam/gate coordinates of a specified radar's field-of-view - - Parameters - ---------- - stid: int - Station ID of radar of choice. See 'superdarn.ca/radar-info' - for ID numbers. - coords: Coords object - Type of coordinates returned - Default: Coords.AACGM - date: datetime - datetime object date to be used for AACGMv2 conversion - Default: Current day - - Returns - ---------- - latitudes: np.array - n_beams x n_gates array of geographic or AACGMv2 latitudes - for range gate corners - longitudes/mlts: np.array - n_beams x n_gates array of geographic or AACGMv2 longitudes - for range gate corners - """ - # Locate base PyDARN directory - if ranges is None: - ranges = [0, SuperDARNRadars.radars[stid].range_gate_45] - if max_beams is None: - max_beams = SuperDARNRadars.radars[stid].hardware_info.beams - # Plus 1 is due to the fact fov files index at 1 so in the plotting - # of the boundary there is a subtraction of 1 to offset this as python - # converts to index of 0 which my code already accounts for - beam_corners_lats = np.zeros((ranges[1], max_beams+1)) - beam_corners_lons = np.zeros((ranges[1], max_beams+1)) - - for beam in range(0, max_beams+1): - for gate in range(ranges[0], ranges[1]): - lat, lon = geographic_cell_positions(stid, beam, gate, rsep, - frang, height=300) - - if coords == Coords.AACGM: - if date is None: - date = dt.datetime.now() - - geomag = np.array(aacgmv2. get_aacgm_coord(lat, lon, - 250, date)) - lat = geomag[0] - lon = geomag[1] - beam_corners_lats[gate, beam] = lat - beam_corners_lons[gate, beam] = lon - - # Return geographic coordinates - return beam_corners_lats, beam_corners_lons - - -# RPosGeo line 335 -def geographic_cell_positions(stid: int, beam: int, range_gate: int, - rsep: int = 45, frang: int = 180, - height: float = None, elv_angle: float = 0.0, - center: bool = True, chisham: bool = False, - virtual_height_type: object = VH_types.STANDARD_VIRTUAL_HEIGHT): - """ - determines the geographic cell position for a given range gate and beam - - parameters - ---------- - stid: int - station id of the radar to use - beam: int - beam number (indexing at 0) - range_gate: int - range gate number (indexing at 0) - rsep: int - range seperation, determined by the mode the - radar is using, in [km] - default: 45 - normalscan mode - frang: int - frequency range from the radar to the front edge of the range gate - Please note: this definition may be changed, currently defined in - RST code to keep consistency - height: float - transmutation height [km] - default: none - if none then it uses elevation angle - elv_angle: float - elevation angle in [deg] - default: 0 - center: bool - obtain geographic location of the centre of the range gates - False obtains the front edge of the range gates. - default: True - virtual_height_type: object - use for choosing type of virtual height - default: VH_types.STANDARD_VIRTUAL_HEIGHT - - returns - ------- - lat: float - latitude of the range gate in geographic coordinates [deg] - lon: float - longitude of the range gate in geographic coordinates [deg] - """ - # centre of the field of view - offset = SuperDARNRadars.radars[stid].hardware_info.beams / 2.0 - 0.5 - - # Obtain radar information from hardware files all converted to [rad] - # radians are needed for numpy geometry calculations and to reduce - # too much converting back and forth. - boresight = np.radians(SuperDARNRadars. - radars[stid].hardware_info.boresight) - radar_lat = np.radians(SuperDARNRadars. - radars[stid].hardware_info.geographic.lat) - radar_lon = np.radians(SuperDARNRadars. - radars[stid].hardware_info.geographic.lon) - beam_sep = np.radians(SuperDARNRadars. - radars[stid].hardware_info.beam_separation) - rxrise = SuperDARNRadars.radars[stid].hardware_info.rx_rise_time - - # TODO: fix after slant range change - if center is True: - # beam edge in [rad] - beam_edge = -beam_sep * 0.5 - # range_edge in [km] - else: - beam_edge = 0 - - # psi [rad] in the angle from the boresight - psi = beam_sep * (beam - offset) + beam_edge - # Calculate the slant range [km] - slant_range = gate2slant(frang, rsep, rxrise, gate=range_gate) - - # If no height is specified then use elevation angle (default 0) - # to calculate the transmutation height - if height is None: - height = -Re + np.sqrt(Re**2 + 2 * slant_range * Re * - np.sin(np.radians(elv_angle)) + slant_range**2) - - lat, lon = geocentric_coordinates(radar_lat, radar_lon, slant_range, - height, psi, boresight, virtual_height_type = virtual_height_type) - - # convert back degrees as preferred units to use? - return np.degrees(lat), np.degrees(lon) +from pydarn import VHModels, EARTH_EQUATORIAL_RADIUS # fldpnth line 90 -def geocentric_coordinates(radar_lat: float, radar_lon: float, - slant_range: float, cell_height: float, - psi: float, boresight: float, - virtual_height_type: object = VH_types.STANDARD_VIRTUAL_HEIGHT): +def geocentric_coordinates(target_range: float, psi: float, boresight: float, + virtual_height_model: VHModels = + VHModels.STANDARD, + **kwargs): """ Calculates the geocentric coordinates of gate cell point, using either the standard or Chisham virtual height model. @@ -211,17 +50,17 @@ def geocentric_coordinates(radar_lat: float, radar_lon: float, radars site latitude [rad] radar_lon : float radars site longitude [lon] - slant_range: float - slant range distance [km] + target_range: float + The range from the instrument to the target (echo) [km] cell_height : float virtual height of the gate cell [km] psi: int [rad] boresight: float boresight of the radar beam [rad] - virtual_height_type: object + virtual_height_model: VHModels use for choosing type of virtual height - default: VH_types.STANDARD_VIRTUAL_HEIGHT + default: VHModels.STANDARD Returns ------- @@ -231,21 +70,16 @@ def geocentric_coordinates(radar_lat: float, radar_lon: float, longitude of the range gate in geographic coordinates [rad] """ - # Gareth Chisham Virtual height model: """ Mapping ionospheric backscatter measured by the SuperDARN HF radars – Part 1: A new empirical virtual height model by G. Chisham 2008 (https://doi.org/10.5194/angeo-26-823-2008) """ - if virtual_height_type == VH_types.CHISHAM: - x_height = chisham(slant_range) - - if virtual_height_type == VH_types.STANDARD_VIRTUAL_HEIGHT: - x_height = standard_virtual_height(slant_range, cell_height) + x_height = virtual_height_model(target_range=target_range, **kwargs) # calculate the radius over the earth underneath # the radar and range gate cell - rlat, rlon, r_radar, delta = geodetic2geocentric(radar_lat, radar_lon) + rlat, rlon, r_radar, delta = geodetic2geocentric(**kwargs) r_cell = r_radar psi_cos_2 = np.cos(psi)**2 @@ -256,14 +90,14 @@ def geocentric_coordinates(radar_lat: float, radar_lon: float, # distance between the gate cell to the earth's centre [km] cell_rho = r_cell + x_height # elevation angle relative to local horizon [rad] - rel_elv = np.arcsin(((cell_rho**2) - (r_radar**2) - slant_range**2) / - (2.0 * r_radar * slant_range)) + rel_elv = np.arcsin(((cell_rho**2) - (r_radar**2) - target_range**2) / + (2.0 * r_radar * target_range)) # estimate elevation for multi-hop propagation - if chisham and slant_range > 2137.5: - gamma = np.arccos((r_radar**2 + cell_rho**2 - slant_range**2) / + if virtual_height_model == VHModels.CHISHAM and target_range > 2137.5: + gamma = np.arccos((r_radar**2 + cell_rho**2 - target_range**2) / (2.0 * r_radar * cell_rho)) beta = np.arcsin(r_radar * np.sin(gamma/3.0) / - (slant_range/3.0)) + (target_range/3.0)) # Elevation angle used for estimating off-array normal # azimuth [rad] xelv = (np.pi/2) - beta - (gamma/3.0) @@ -288,14 +122,18 @@ def geocentric_coordinates(radar_lat: float, radar_lon: float, # azimuth of the gate cell [rad] cell_azimuth = azimuth + boresight - flatten_azimuth = geocentric2flattening(delta, cell_azimuth, xelv) + flatten_azimuth = geocentric2flattening(delta=delta, + azimuth=cell_azimuth, + elv=xelv) cell_rho, cell_lat, cell_lon = \ - cell_geocentric_coordinates(rlat, rlon, r_radar, - flatten_azimuth, rel_elv, - slant_range) + cell_geocentric_coordinates(lat=rlat, lon=rlon, + rho=r_radar, + azimuth=flatten_azimuth, + elv=rel_elv, + r=target_range) # recalculate the radius under the gate cell and centre of earth - r_cell = geocentric2geodetic(cell_lat, cell_lon) + r_cell = geocentric2geodetic(lat=cell_lat, lon=cell_lon) cell_heightx = cell_rho - r_cell # this ensures convergence on the cell point while_flag = abs(cell_heightx - x_height) > 0.5 @@ -305,7 +143,8 @@ def geocentric_coordinates(radar_lat: float, radar_lon: float, # fldpnt def cell_geocentric_coordinates(lat: float, lon: float, rho: float, - azimuth: float, elv: float, r: float): + azimuth: float, elv: float, r: float, + **kwargs): """ Calculates the geocentric coordinates of a gate cell point given the angular geocentric coordinates of the point of origin, @@ -385,7 +224,7 @@ def cell_geocentric_coordinates(lat: float, lon: float, rho: float, # goecnvrt -def geocentric2flattening(delta: float, azimuth: float, elv: float): +def geocentric2flattening(delta: float, azimuth: float, elv: float, **kwargs): """ Adjust azimuth for the oblateness of the Earth @@ -426,7 +265,7 @@ def geocentric2flattening(delta: float, azimuth: float, elv: float): # geodtgc, iopt > 0 -def geodetic2geocentric(lat: float, lon: float): +def geodetic2geocentric(lat: float, lon: float, **kwargs): """ convert geodetic coordinates to geocentric @@ -469,7 +308,7 @@ def geodetic2geocentric(lat: float, lon: float): # geodtgc, iopt < 0 -def geocentric2geodetic(lat: float, lon: float): +def geocentric2geodetic(lat: float, lon: float, **kwargs): """ convert geocentric coordinates to geodetic diff --git a/pydarn/utils/hdw b/pydarn/utils/hdw index 34294ef6..f3c8bcee 160000 --- a/pydarn/utils/hdw +++ b/pydarn/utils/hdw @@ -1 +1 @@ -Subproject commit 34294ef63092442f6bbae227b6a7687a108e2127 +Subproject commit f3c8bceeac763146aa7c5eb209277034b3a6b7e2 diff --git a/pydarn/utils/plotting.py b/pydarn/utils/plotting.py index 42898386..25aff0ce 100644 --- a/pydarn/utils/plotting.py +++ b/pydarn/utils/plotting.py @@ -1,10 +1,21 @@ # Copyright (C) 2020 SuperDARN Canada, Universtiy of Saskatchewan # Author(s): Marina Schmidt +# +# Disclaimer: +# pyDARN is under the LGPL v3 license found in the root directory LICENSE.md +# Everyone is permitted to copy and distribute verbatim copies of this license +# document, but changing it is not allowed. +# +# This version of the GNU Lesser General Public License incorporates the terms +# and conditions of version 3 of the GNU General Public License, +# supplemented by the additional permissions listed below. +# +# Modification: """ This module is utility functions that are useful for multiple plotting methods """ - +import enum import numpy as np from datetime import datetime @@ -12,6 +23,48 @@ from pydarn import plot_exceptions +class MapParams(enum.Enum): + """ + """ + FITTED_VELOCITY = "fitted" + MODEL_VELOCITY = "model.vel.median" + RAW_VELOCITY = "vector.vel.median" + POWER = "vector.pwr.median" + SPECTRAL_WIDTH = "vector.wdt.median" + + +def find_record(dmap_data: List[dict], start_time: datetime, time_delta:int = 1): + """ + finds the record number that associates to the start time + + Parameter + --------- + dmap_data : List[dict] + the data to look over for the record number + start_time : datetime + the start_time to associate to the record number + time_delta : int + the difference between start_time and dmap_data time to determine + the record number within a region + + Return + ------ + record_num : int + the record number associated with the start_time + + Raises + ------ + NoDataFound + raises if the start_time is not in the dmap_data list + """ + record_num = 0 + for record in dmap_data: + date = time2datetime(record) + time_diff = date - start_time + if time_diff.seconds/60 <= time_delta: + return record_num + record_num += 1 + raise plot_exceptions.NoDataFoundError(parameter, start_time=start_time) def check_data_type(dmap_data: List[dict], parameter: str, expected_type: str, index: int): @@ -61,13 +114,22 @@ def time2datetime(dmap_record: dict) -> datetime: datetime object returns a datetime object of the records time stamp """ - year = dmap_record['time.yr'] - month = dmap_record['time.mo'] - day = dmap_record['time.dy'] - hour = dmap_record['time.hr'] - minute = dmap_record['time.mt'] - second = dmap_record['time.sc'] - micro_sec = dmap_record['time.us'] - - return datetime(year=year, month=month, day=day, hour=hour, - minute=minute, second=second, microsecond=micro_sec) + try: + year = dmap_record['time.yr'] + month = dmap_record['time.mo'] + day = dmap_record['time.dy'] + hour = dmap_record['time.hr'] + minute = dmap_record['time.mt'] + second = dmap_record['time.sc'] + micro_sec = dmap_record['time.us'] + return datetime(year=year, month=month, day=day, hour=hour, + minute=minute, second=second, microsecond=micro_sec) + except KeyError: + year = dmap_record['start.year'] + month = dmap_record['start.month'] + day = dmap_record['start.day'] + hour = dmap_record['start.hour'] + minute = dmap_record['start.minute'] + second = dmap_record['start.second'] + return datetime(year=year, month=month, day=day, hour=hour, + minute=minute, second=int(second)) diff --git a/pydarn/utils/range_estimations.py b/pydarn/utils/range_estimations.py new file mode 100644 index 00000000..736d0ac6 --- /dev/null +++ b/pydarn/utils/range_estimations.py @@ -0,0 +1,126 @@ +# (C) Copyright 2021 SuperDARN Canada, University of Saskatchewan +# Author(s): Marina Schmidt +# +# Disclaimer: +# pyDARN is under the LGPL v3 license found in the root directory LICENSE.md +# Everyone is permitted to copy and distribute verbatim copies of this license +# document, but changing it is not allowed. +# +# This version of the GNU Lesser General Public License incorporates the terms +# and conditions of version 3 of the GNU General Public License, +# supplemented by the additional permissions listed below. +# +# Modifications: +# + +import enum +import numpy as np +import math + +from pydarn import Re, C + + +def gate2groundscatter(reflection_height: float = 250, **kwargs): + """ + Calculate the ground scatter mapped range (km) for each slanted range + for SuperDARN data. This function is based on the Ground Scatter equation + from Bristow paper at https://doi.org/10.1029/93JA01470 on page 325 + Parameters + ---------- + reflection_height: float + reflection height + default: 250 + + Returns + ------- + ground_scatter_mapped_ranges : np.array + returns an array of ground scatter mapped ranges for the radar + """ + slant_ranges = gate2slant(**kwargs) + ground_scatter_mapped_ranges =\ + Re*np.arcsin(np.sqrt((slant_ranges**2/4)- + (reflection_height**2))/Re) + + return ground_scatter_mapped_ranges + + +def gate2slant(rxrise: int = 0, range_gate: int = 0, frang: int = 180, + rsep: int = 45, nrang: int = None, center: bool = True, + **kwargs): + """ + Calculate the slant range (km) for each range gate for SuperDARN data + + Parameters + ---------- + frang: int + range from the edge of first the gate to the radar [km] + This should be given in fitacf record of the control program + rsep: int + Radar seperation of the gates. Determined by control program. + rxrise: int + Use hardware value for this, avoid data file values + gate: int + range gate to determine the slant range [km], if nrang + is None + default: 0 + nrang: int + max number of range gates in the list of records. If + not None, will calculate all slant ranges + default: None + center: boolean + Calculate the slant range in the center of range gate + or edge + + Returns + ------- + slant_ranges : np.array + returns an array of slant ranges for the radar + """ + # lag to the first range gate in microseconds + # 0.3 - speed of light (km/us) + # 2 - two times for there and back + distance_factor = 2.0 + # C - speed of light m/s to km/us + speed_of_light = C * 0.001 * 1e-6 + lag_first = frang * distance_factor / speed_of_light + # sample separation in microseconds + sample_sep = rsep * distance_factor / speed_of_light + # Range offset + # If center is true, calculate at the center + if center: + # 0.5 off set to the centre of the range gate instead of edge + range_offset = -0.5 * rsep + else: + range_offset = 0.0 + # Now calculate slant range in km + if nrang is None: + slant_ranges = (lag_first - rxrise + + range_gate * sample_sep) * speed_of_light /\ + distance_factor + range_offset + else: + slant_ranges = np.zeros(nrang+1) + for gate in range(nrang+1): + slant_ranges[gate] = (lag_first - rxrise + + gate * sample_sep) * speed_of_light /\ + distance_factor + range_offset + return slant_ranges + + +class RangeEstimation(enum.Enum): + """ + Range_Estimation class is to list the current range gate estimations + a user can pick from + + enumerators: + RANGE_GATE: range gates + SLANT_RANGE: slant range (km) + GSMR: ground scatter mapped range (km) + """ + + RANGE_GATE = enum.auto() + SLANT_RANGE = (gate2slant,) + GSMR = (gate2groundscatter,) + + # Need this to make the functions callable + def __call__(self, *args, **kwargs): + return self.value[0](*args, **kwargs) diff --git a/pydarn/utils/scan.py b/pydarn/utils/scan.py index 970577a8..231d1da8 100644 --- a/pydarn/utils/scan.py +++ b/pydarn/utils/scan.py @@ -1,5 +1,17 @@ # Copyright (C) 2020 SuperDARN Canada, Universtiy of Saskatchewan # Author(s): Daniel Billett +# +# Disclaimer: +# pyDARN is under the LGPL v3 license found in the root directory LICENSE.md +# Everyone is permitted to copy and distribute verbatim copies of this license +# document, but changing it is not allowed. +# +# This version of the GNU Lesser General Public License incorporates the terms +# and conditions of version 3 of the GNU General Public License, +# supplemented by the additional permissions listed below. +# +# Modification: +# """ This module is used for sorting a given dmap_data list of dictionaries into scans diff --git a/pydarn/utils/superdarn_cpid.py b/pydarn/utils/superdarn_cpid.py index 25c85470..28a74da7 100644 --- a/pydarn/utils/superdarn_cpid.py +++ b/pydarn/utils/superdarn_cpid.py @@ -1,6 +1,17 @@ # Copyright (C) 2019 SuperDARN Canada, University Saskatchewan # Author: Marina Schmidt - +# +# Disclaimer: +# pyDARN is under the LGPL v3 license found in the root directory LICENSE.md +# Everyone is permitted to copy and distribute verbatim copies of this license +# document, but changing it is not allowed. +# +# This version of the GNU Lesser General Public License incorporates the terms +# and conditions of version 3 of the GNU General Public License, +# supplemented by the additional permissions listed below. +# +# Modification: +# """ This module contains SuperDARN radar modes classified by cpid number Current documentation on modes: diff --git a/pydarn/utils/superdarn_radars.py b/pydarn/utils/superdarn_radars.py index 0be4980d..6e98e236 100644 --- a/pydarn/utils/superdarn_radars.py +++ b/pydarn/utils/superdarn_radars.py @@ -1,18 +1,18 @@ # Copyright (C) 2020 SuperDARN Canada, University of Saskatchewan # Authors: Marina Schmidt and Danno Peters # -# This file is part of the pyDARN Library. -# +# Disclaimer: # pyDARN is under the LGPL v3 license found in the root directory LICENSE.md # Everyone is permitted to copy and distribute verbatim copies of this license # document, but changing it is not allowed. # # This version of the GNU Lesser General Public License incorporates the terms -# and conditions of version 3 of the GNU Lesser General Public License, +# and conditions of version 3 of the GNU General Public License, # supplemented by the additional permissions listed below. # # Modifications: -# +# 2022-02-11 MTS USASK updated the _HDW_info class to take in +# the hardware format """ This module contains SuperDARN radar information """ @@ -23,7 +23,7 @@ from typing import NamedTuple from enum import Enum -from datetime import datetime, timedelta +from datetime import datetime from subprocess import check_call @@ -56,19 +56,19 @@ def get_hdw_files(force: bool = True, version: str = None): if len(os.listdir(hdw_path)) == 0 or force: # pycurl doesn't download a zip folder easily so # use the command line command - check_call(['curl', '-L', '-o', hdw_path+'/master.zip', - 'https://github.com/SuperDARN/hdw/archive/master.zip']) + check_call(['curl', '-L', '-o', hdw_path+'/main.zip', + 'https://github.com/SuperDARN/hdw/archive/main.zip']) # use unzip command because zipfile on works with files and not folders # though this is possible with zipfile but this was easier for me to # get it working - check_call(['unzip', '-d', hdw_path, hdw_path+'/master.zip']) - dat_files = glob.glob(hdw_path+'/hdw-master/*') + check_call(['unzip', '-d', hdw_path, hdw_path+'/main.zip']) + dat_files = glob.glob(hdw_path+'/hdw-main/*') # shutil only moves specific files so we need to move # everything one at a time for hdw_file in dat_files: shutil.move(hdw_file, hdw_path+os.path.basename(hdw_file)) # delete the empty folder - os.removedirs(hdw_path+'/hdw-master/') + os.removedirs(hdw_path+'/hdw-main/') def read_hdw_file(abbrv, date: datetime = None, update: bool = False): @@ -101,91 +101,147 @@ def read_hdw_file(abbrv, date: datetime = None, update: bool = False): hdw_path = os.path.dirname(__file__)+'/hdw/' hdw_file = "{path}/hdw.dat.{radar}".format(path=hdw_path, radar=abbrv) + hdw_data = [] + hdw_lines_date = [] # if the file does not exist then try # and download it if os.path.exists(hdw_file) is False: get_hdw_files() try: with open(hdw_file, 'r') as reader: - for line in reader.readlines(): - if '#' not in line and len(line.split()) > 1: - hdw_data = line.split() + lines = reader.readlines() + for i in range(len(lines)): + if '#' not in lines[i] and len(lines[i].split()) > 1: + hdw_data.append(lines[i].split()) """ Hardware files give the year and seconds from the beginning of that year. Thus to check the date if it corresponds we need to convert to a datetime object and then compare. """ - hdw_line_date = datetime(year=int(hdw_data[1]), month=1, - day=1) +\ - timedelta(seconds=int(hdw_data[2])) - if hdw_line_date > date: - """ - Hardware data array positions definitions: - 0: station id - stid - 1: last year that the parameter string is valid. - Note: currently updated line will - have a year of 2999 - meaning it is currently still up to date. - 2: last second of year that parameter - string is valid. - 3: Geographic latitude of radar site - 4: Geographic longitude of radar site - Note: southern lat and long are negative - 5: Altitude of the radar site (meters) - 6: Scanning boresight - direction of - the centre beam, - measured in degrees relative to geographic north. - Counter clockwise rotations are negative. - 7: Beam separation (Angular separation in degrees) - 8: velocity sign - at radar level, - backscattered signal with - frequencies above the transmitted frequency - are assigned positive - Doppler velocities while backscattered signals - with frequencies - below the transmitted frequency are assigned - negative Doppler - velocity. Can be changed in receiver design. - 9: Analog Rx attenuator step (dB) - 10: Tdiff - propagation time from - interferometer array antenna - to phasing matrix input miunus propagation - time from main array - antenna through transmitter to phasing - matrix input. - (microseconds) - 11: phase sign - to account for any cable errors - Interferometer offset - displacement of midpoint - interferometer array from midpoint main - array (meters). - 12: x direction - along the line of antennas - with +X toward - higher antenna number - 13: y direction - along the array normal - with +Y in the - direction of the array normal - 14: z direction - is the altitude difference, +Z up - 15: Analog Rx rise time (microseconds) - 16: Analog Attenuation stages - gain control of - an analog - receiver or front-end - 17: maximum range gates - 18: maximum number of beams - """ - return _HdwInfo(int(hdw_data[0]), abbrv, - _Coord(float(hdw_data[3]), - float(hdw_data[4]), - float(hdw_data[5])), - float(hdw_data[6]), float(hdw_data[7]), - float(hdw_data[8]), float(hdw_data[9]), - float(hdw_data[10]), - float(hdw_data[11]), - _InterferometerOffset(float(hdw_data[12]), - float(hdw_data[13]), - float(hdw_data[14])), - float(hdw_data[15]), - float(hdw_data[16]), - int(hdw_data[17]), int(hdw_data[18])) + j = len(hdw_data)-1 + hdw_lines_date.append( + datetime(year=int(hdw_data[j][2][0:4]), + month=int(hdw_data[j][2][4:6]), + day=int(hdw_data[j][2][6:8]), + hour=int(hdw_data[j][3][0:2]), + minute=int(hdw_data[j][3][3:5]), + second=int(hdw_data[j][3][6:8]))) + if hdw_lines_date[j] > date: + j = j-1 + break + """ + Hardware data array positions definitions: + 0: Station ID (unique numerical value). + 1: Status code (1 operational, -1 offline). + 2: First date that parameter string is valid + (YYYYMMDD). + 3: First time that parameter string is valid + (HH:MM:SS). + 4: Geographic latitude of radar site + (Given in decimal degrees to 3 + decimal places. Southern hemisphere + values are negative) + 5: Geographic longitude of radar site + (Given in decimal degrees to + 3 decimal places. + West longitude values are negative) + 6: Altitude of the radar site (meters) + 7: Physical scanning boresight + (Direction of the center beam, measured in + degrees relative to geographic north. + CCW rotations are negative.) + 8: Electronic shift to radar scanning + boresight (Degrees relative to + physical antenna boresight. + Normally 0.0 degrees) + 9: Beam separation (Angular + separation in degrees between adjacent + beams. Normally 3.24 degrees) + 10: Velocity sign (At the radar level, + backscattered signals with + frequencies above the transmitted + frequency are assigned positive + Doppler velocities while backscattered + signals with frequencies below + the transmitted frequency are assigned + negative Doppler velocity. This + convention can be reversed by changes + in receiver design or in the + data sampling rate. This parameter + is set to +1 or -1 to maintain the + convention.) + 11: Phase sign (Cabling errors can + lead to a 180 degree shift of the + interferometry phase measurement. + +1 indicates that the sign is + correct, -1 indicates that it must be flipped.) + 12: Tdiff [Channel A] + (Propagation time from interferometer + array antenna to phasing matrix input + minus propagation time from main array antenna + through transmitter to phasing matrix input. + Units are decimal + microseconds) + 13: Tdiff [Channel B] + (Propagation time from interferometer + array antenna to phasing matrix input minus + propagation time from main array antenna + through transmitter to phasing matrix input. + Units are decimal microseconds) + 14: Interferometer X offset + (Displacement of midpoint of interferometer + array from midpoint of main array, + along the line of antennas + with +X toward higher antenna numbers. + Units are meters) + 15: Interferometer Y offset + (Displacement of midpoint of + interferometer array from midpoint of + main array, along the array + normal direction with +Y in the direction of + the array normal. Units are meters) + 16: Interferometer Z offset + (Displacement of midpoint of + interferometer array from midpoint of + main array, in terms of altitude + difference with +Z up. Units are meters) + 17: Analog Rx rise time + (Time given in microseconds. Time delays of + less than ~10 microseconds can be ignored. + If narrow-band filters are + used in analog receivers or front-ends, + the time delays should be + specified.) + 18: Analog Rx attenuator step (dB) + 19: Analog attenuation stages (Number of stages. + This is used for gain control of an analog + receiver or front-end.) + 20: Maximum of range gates used + 21: Maximum number of beams + """ + return _HdwInfo(stid=int(hdw_data[j][0]), + status=Status(int(hdw_data[j][1])), + abbrev=abbrv, + date=hdw_lines_date[j], + geographic=_Coord(float(hdw_data[j][4]), + float(hdw_data[j][5]), + float(hdw_data[j][6])), + boresight=_Boresight(float(hdw_data[j][7]), + float(hdw_data[j][8])), + beam_separation=float(hdw_data[j][9]), + velocity_sign=float(hdw_data[j][10]), + phase_sign=float(hdw_data[j][11]), + tdiff=_Tdiff(float(hdw_data[j][12]), + float(hdw_data[j][13])), + interferometer_offset=_InterferometerOffset( + float(hdw_data[j][14]), + float(hdw_data[j][15]), + float(hdw_data[j][16])), + rx_rise_time=float(hdw_data[j][17]), + rx_attenuator=float(hdw_data[j][18]), + attenuation_stages=int(hdw_data[j][19]), + gates=int(hdw_data[j][20]), + beams=int(hdw_data[j][21])) except FileNotFoundError: raise pydarn.radar_exceptions.HardwareFileNotFoundError(abbrv) @@ -208,6 +264,19 @@ class Hemisphere(Enum): South = -1 +class Status(Enum): + """ + this class is to give names to status of radars + + Attributes + ---------- + online : 1 + offline : -1 + """ + online = 1 + offline = -1 + + class _InterferometerOffset(NamedTuple): """ Named tuple class to contain the interferometer offset @@ -248,20 +317,55 @@ class _Coord(NamedTuple): alt: float +class _Tdiff(NamedTuple): + """ + Class use to store tdiff values for various channels and + maybe for common frequencies. + + Attributes + ---------- + channel_a : float + for channel a + channel_b : float + for channel b + """ + channel_a: float + channel_b: float + + +class _Boresight(NamedTuple): + """ + This class used to store the various boresights + + Attributes + ---------- + physical : float + Physical scanning boresight + electronic : float + Electronic shift to the radar scanning boresight + """ + physical: float + electronic: float + + class _HdwInfo(NamedTuple): """ Class used to store relevant information about the SuperDARN Radars Attributes - ---------- + --- ------- stid : int station number + status : Status + status of radars operation abbrev : str three letter station abbreviation + date : datetime + date of the hdw file data is being used geographic : _Coord object Named Tuple containing geographic latitude longitude and altitude - boresight : float - boresight center beam in degrees + boresight : _Boresight + boresight physical and electronic in degrees beam_separation : float angular separation between radar beams velocity_sign : float @@ -270,19 +374,20 @@ class _HdwInfo(NamedTuple): Doppler velocities while backscattered signals with frequencies below the transmitted frequency are assigned negative Doppler velocity. Can be changed in receiver design. - rx_attenuator : float - Analog Rx attenuator step (dB) - tdiff : float + phase_sign : float + To account for flipped cable errors + tdiff : _Tdiff + channel a and b: propagation time from interferometer array antenna to phasing matrix input minus propagation time from main array antenna through transmitter to phasing matrix input. (microseconds) - phase_sign : float - To account for flipped cable errors interferometer_offset : _InterferometerOffset displacement of midpoint interferometer array from midpoint main array in Cartesian coordinates(meters). + rx_attenuator : float + Analog Rx attenuator step (dB) rx_rise : float Analog Rx rise time (microseconds) attenuation_stages : float @@ -296,16 +401,22 @@ class _HdwInfo(NamedTuple): See Also -------- read_hdw_file : function for reading hardware files + Status _Coord : object contain coordinate information + _Boresight + _Tdiff + _Interferometer """ stid: int + status: int abbrev: str + date: datetime geographic: _Coord - boresight: float + boresight: _Boresight beam_separation: float velocity_sign: float rx_attenuator: float - tdiff: float + tdiff: _Tdiff phase_sign: float interferometer_offset: _InterferometerOffset rx_rise_time: float diff --git a/pydarn/utils/virtual_heights.py b/pydarn/utils/virtual_heights.py index c8097d6f..7a0650e4 100644 --- a/pydarn/utils/virtual_heights.py +++ b/pydarn/utils/virtual_heights.py @@ -1,25 +1,40 @@ -# (C) Copyright 2021 SuperDARN Canada, University of Saskatachewan -# Author(s): Marina Schmidt -# -# This file is part of the pyDARN Library. +# (C) Copyright 2021 SuperDARN Canada, University of Saskatachewan +# Author(s): Marina Schmidt +# (C) Copyright 2021 University of Scranton +# Author(s): Francis Tholley # +# Disclaimer: # pyDARN is under the LGPL v3 license found in the root directory LICENSE.md # Everyone is permitted to copy and distribute verbatim copies of this license # document, but changing it is not allowed. # # This version of the GNU Lesser General Public License incorporates the terms -# and conditions of version 3 of the GNU Lesser General Public License, +# and conditions of version 3 of the GNU General Public License, # supplemented by the additional permissions listed below. # # Modifications: -# 2021-09-15 Francis Tholley moved the chisham and standard virtual heigh models to separate file for better encapsulation/modularity -""" virtual_heights.py comprises of different types of virtual heights""" +# 2021-09-15 Francis Tholley moved the chisham and standard virtual +# height models to separate file for better encapsulation/modularity +# 2022-03-04 Marina Schmidt add the VH_Types class to the bottom +""" virtual_heights.py comprises of different of virtual height models""" +import enum -def chisham(slant_range: float): +def chisham(target_range: float, **kwargs): """ Mapping ionospheric backscatter measured by the SuperDARN HF radars – Part 1: A new empirical virtual height model by G. Chisham 2008 (https://doi.org/10.5194/angeo-26-823-2008) + + Parameters + ---------- + target_range: float + is the range from radar to the target (echos) + sometimes known as slant range [km] + kwargs: is only needed to avoid key item errors + + Returns + ------- + altered target_range (slant range) [km] """ # Model constants A_const = (108.974, 384.416, 1098.28) @@ -27,22 +42,23 @@ def chisham(slant_range: float): C_const = (6.68283e-5, 1.81405e-4, 9.39961e-5) # determine which region of ionosphere the gate - if slant_range < 115: - return (slant_range / 115.0) * 112.0 - elif slant_range < 787.5: - return A_const[0] + B_const[0] * slant_range + C_const[0] *\ - slant_range**2 - elif slant_range <= 2137.5: - return A_const[1] + B_const[1] * slant_range + C_const[1] *\ - slant_range**2 + if target_range < 115: + return (target_range / 115.0) * 112.0 + elif target_range < 787.5: + return A_const[0] + B_const[0] * target_range + C_const[0] *\ + target_range**2 + elif target_range <= 2137.5: + return A_const[1] + B_const[1] * target_range + C_const[1] *\ + target_range**2 else: - return A_const[2] + B_const[2] * slant_range + C_const[2] *\ - slant_range**2 - - -def standard_virtual_height(slant_range: float, cell_height: float): + return A_const[2] + B_const[2] * target_range + C_const[2] *\ + target_range**2 + + +def standard_virtual_height(target_range: float, cell_height: int = 300, + **kwargs): """ - cell_height, slant_range and x_height are in km + cell_height, target_range and x_height are in km Default values set in virtual height model described Mapping ionospheric backscatter measured by the SuperDARN HF radars – Part 1: A new empirical virtual height model by @@ -52,19 +68,50 @@ def standard_virtual_height(slant_range: float, cell_height: float): 150 - 600 km E region scatter (Note in the paper 400 km is the edge of the E region) 600 - 800 km is F region + + Parameters + ---------- + target_range: float + is the range from radar to the target (echos) + sometimes known as slant range [km] + cell_height: int + the default height of the echo if the target_range + is within a certain range + kwargs: is only needed to avoid key item errors + + Returns + ------- + altered target_range (slant range) [km] """ # TODO: why 115? # map everything into the E region - if cell_height <= 150 and slant_range > 150: + if cell_height <= 150 and target_range > 150: return cell_height # virtual height equation (1) from the above paper - elif slant_range < 150: - return (slant_range / 150.0) * 115 - elif slant_range >= 150 and slant_range <= 600: + elif target_range < 150: + return (target_range / 150.0) * 115 + elif target_range >= 150 and target_range <= 600: return 115 - elif slant_range > 600 and slant_range < 800: - return (slant_range - 600) / 200 * (cell_height - 115) + 115 + elif target_range > 600 and target_range < 800: + return (target_range - 600) / 200 * (cell_height - 115) + 115 # higher than 800 km else: - return cell_height - + return cell_height + + +class VHModels(enum.Enum): + """ + This virtual height models class is to list the current + virtual height model user can pick from + + enumerators: + STANDARD: Standard_Virtual_height (km) + CHISHAM: chisham (km) + """ + + STANDARD = (standard_virtual_height, ) + CHISHAM = (chisham, ) + + # Need this to make the functions callable + def __call__(self, *args, **kwargs): + return self.value[0](*args, **kwargs) diff --git a/pydarn/utils/virtual_heights_types.py b/pydarn/utils/virtual_heights_types.py deleted file mode 100644 index 5f9ab5ad..00000000 --- a/pydarn/utils/virtual_heights_types.py +++ /dev/null @@ -1,34 +0,0 @@ -# (C) Copyright 2021 University of Scranton -# Author(s): Francis Tholley -# -# This file is part of the pyDARN Library. -# -# pyDARN is under the LGPL v3 license found in the root directory LICENSE.md -# Everyone is permitted to copy and distribute verbatim copies of this license -# document, but changing it is not allowed. -# -# This version of the GNU Lesser General Public License incorporates the terms -# and conditions of version 3 of the GNU Lesser General Public License, -# supplemented by the additional permissions listed below. -# -# Modifications: -# - -""" -virtual_heights_types.py is a group of methods that focus on the type of virtual heights -""" -import enum - -class VH_types(enum.Enum): - """ - This virtual height types class is to list the current virtual height systems - a user can pick from - - enumerators: - STANDARD_VIRTUAL_HEIGHT: Standard_Virtual_height (km) - CHISHAM: chisham (km) - """ - - STANDARD_VIRTUAL_HEIGHT = enum.auto() - CHISHAM = enum.auto() - diff --git a/pydarn/version.py b/pydarn/version.py index af97f0d2..4fc93422 100644 --- a/pydarn/version.py +++ b/pydarn/version.py @@ -17,4 +17,4 @@ This file contains the version number of pydarn """ -__version__='2.2.1' +__version__='3.0' diff --git a/setup.py b/setup.py index 18dce25b..511cc067 100644 --- a/setup.py +++ b/setup.py @@ -1,15 +1,25 @@ +#Copyright 2018 SuperDARN Canada, University of Saskatchewan +# Author(s): Marina Schmidt +# +# Disclaimer: +# pyDARN is under the LGPL v3 license found in the root directory LICENSE.md +# Everyone is permitted to copy and distribute verbatim copies of this license +# document, but changing it is not allowed. +# +# This version of the GNU Lesser General Public License incorporates the terms +# and conditions of version 3 of the GNU General Public License, +# supplemented by the additional permissions listed below. +# +# Modifications +# 2022-03-10 MTS - removed radar_fov files +#setup.py +#2018-11-05 """ -Copyright 2018 SuperDARN Canada, University of Saskatchewan - -setup.py -2018-11-05 To setup pyDARN as a third party library. Include installing need libraries for running the files. - -author: -Marina Schmidt """ import sys +import warnings from os import path from setuptools import setup, find_packages @@ -53,6 +63,11 @@ def run(self): long_description = f.read() exec(open('pydarn/version.py').read()) +warnings.warn("If you are going to use Fan, Grid, and/or Convection Map " + "plots, then make sure cartopy is installed on your machine. " + "If you do not need to use cartopy for your plotting, ignore " + "this message.") + # Setup information setup( cmdclass={'install': initialize_submodules, 'develop': update_submodules}, @@ -61,16 +76,15 @@ def run(self): long_description=long_description, long_description_content_type='text/markdown', description="Data visualization library for SuperDARN data", - url='https://github.com/SuperDARN/pydarn.git', + url='https://pydarn.readthedocs.io/en/latest/', classifiers=[ 'Development Status :: 5 - Production/Stable', 'License :: OSI Approved :: GNU Lesser General Public License v3 (LGPLv3)', 'Programming Language :: Python :: 3.6', 'Programming Language :: Python :: 3.7'], python_requires='>=3.6', - data_files=[('pydarn/radar_fov_files',glob('radar_fov_files/**'))], packages=find_packages(exclude=['docs', 'test']), - author="SuperDARN", + author="SuperDARN Data Visualization Working Group", # used to import the logging config file into pydarn. include_package_data=True, # setup_requires=['pyyaml', 'numpy', 'matplotlib', 'aacgmv2'], diff --git a/test/test_RTP.py b/test/test_RTP.py index 18ca5f4c..453ee5ae 100644 --- a/test/test_RTP.py +++ b/test/test_RTP.py @@ -37,7 +37,7 @@ def test_normal_time_series(self): with warnings.catch_warnings(record=True): pydarn.RTP.plot_time_series(data) - def test_normal_time_series(self): + def test_normal_summary_plot(self): """ """ with warnings.catch_warnings(record=True): pydarn.RTP.plot_summary(data) @@ -49,8 +49,9 @@ def test_normal_time_series(self): @pytest.mark.parametrize('groundscatter_params', [True, 'grey']) @pytest.mark.parametrize('colorbar_label', ['test']) @pytest.mark.parametrize('yspacing', [150, 250]) -@pytest.mark.parametrize('coords', [pydarn.Coords.RANGE_GATE, - pydarn.Coords.GROUND_SCATTER_MAPPED_RANGE]) +@pytest.mark.parametrize('range_estimation', + [pydarn.RangeEstimation.RANGE_GATE, + pydarn.RangeEstimation.GSMR]) @pytest.mark.parametrize('ymin', [10]) @pytest.mark.parametrize('ymax', [70]) @pytest.mark.parametrize('parameters', ['v', 'p_l', 'w_l']) @@ -64,7 +65,7 @@ def test_parameters_range_time(self, groundscatter_params, parameters, background, zmin, zmax, start_time, end_time, ymin, ymax, colorbar_label, yspacing, - coords, cmap, date_fmt): + range_estimation, cmap, date_fmt): """ this test will give bare minimum input needed for """ with warnings.catch_warnings(record=True): pydarn.RTP.plot_range_time(data, beam_num=15, @@ -74,7 +75,8 @@ def test_parameters_range_time(self, groundscatter_params, zmax=zmax, start_time=start_time, end_time=end_time, ymin=ymin, ymax=ymax, colorbar_label=colorbar_label, - yspacing=yspacing, coords=coords, + yspacing=yspacing, + range_estimation=range_estimation, cmap=cmap, date_fmt=date_fmt) plt.close('all') @@ -82,7 +84,9 @@ def test_parameters_channel_2_range_time(self, groundscatter_params, parameters, background, zmin, zmax, start_time, end_time, ymin, ymax, colorbar_label, - yspacing, coords, cmap, date_fmt): + yspacing, + range_estimation, + cmap, date_fmt): """ this test will give bare minimum input needed for """ with warnings.catch_warnings(record=True): pydarn.RTP.plot_range_time(data, beam_num=9, @@ -92,7 +96,8 @@ def test_parameters_channel_2_range_time(self, groundscatter_params, zmax=zmax, start_time=start_time, end_time=end_time, ymin=ymin, ymax=ymax, colorbar_label=colorbar_label, - yspacing=yspacing, coords=coords, + yspacing=yspacing, + range_estimation=range_estimation, cmap=cmap, date_fmt=date_fmt) plt.close('all') @@ -159,30 +164,34 @@ def test_parameters_time_series_channel1(self, parameters_scalar, gate, @pytest.mark.parametrize('title', ['test test']) @pytest.mark.parametrize('background', ['grey']) @pytest.mark.parametrize('groundscatter', [False, 'grey']) -@pytest.mark.parametrize('coords', [pydarn.Coords.RANGE_GATE, - pydarn.Coords.GROUND_SCATTER_MAPPED_RANGE]) +@pytest.mark.parametrize('range_estimation', + [pydarn.RangeEstimation.RANGE_GATE, + pydarn.RangeEstimation.GSMR]) class TestSummaryPlots: def test_params_summary_plots(self, fig_size, watermark, boundary, cmaps, lines, plot_elv, title, background, - groundscatter, coords): + groundscatter, range_estimation): with warnings.catch_warnings(record=True): pydarn.RTP.plot_summary(data, beam_num=15, channel=1, figsize=fig_size, watermark=watermark, boundary=boundary, cmaps=cmaps, lines=lines, plot_elv=plot_elv, title=title, background=background, - groundscatter=groundscatter, coords=coords) + groundscatter=groundscatter, + range_estimation=range_estimation) plt.close('all') def test_beam9_channel2_summary_plots(self, fig_size, watermark, boundary, cmaps, lines, plot_elv, title, - background, groundscatter, coords): + background, groundscatter, + range_estimation): with warnings.catch_warnings(record=True): pydarn.RTP.plot_summary(data, beam_num=9, channel=2, figsize=fig_size, watermark=watermark, boundary=boundary, cmaps=cmaps, lines=lines, plot_elv=plot_elv, title=title, background=background, - groundscatter=groundscatter, coords=coords) + groundscatter=groundscatter, + range_estimation=range_estimation) plt.close('all')