<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Reading and plotting LMA data
In a previous module we learned about LMA VHF source files - their structure, and principles of data quality control using the station mask and $\chi_{\nu}^2$.
We'll use the xlma-python package to read, process, and visualize the data in these files.
We'll start by importing some libraries and reading a file. In fact, we can read more than Replace the filename below with a file of interest to you, if you want.
Step1:  Investigating the pyxlma data structure
We just read in the data from a list of filenames, and got back an lma_data variable, and a starttime. What's in these variables?
starttime is a Python datetime object, with helps us do math operations on times instead of multiplying by 60 and 24 ourselves.
Step2:  lma_data is an xarray object. If we print it, we see that it looks much like a NetCDF file, with dimensions and variables that allow us to store whole arrays of data and give them names. xarray is the best way to look at NetCDF data in Python.
(If you'd like to learn more about xarray and NetCDF, you can check out Unidata's lessons on xarray and NetCDF.)
Step3:  There are a few things to notice above.
- The dimensions tell us how many events and stations we had in our data file.
- Some other header info about the network location and the lma_analysis command is included.
- There are varaibles with dimension number_of_stations that give the data from the station data tables.
- There are varaibles with dimension number_of_events, including many we'd expect.
    - Event location and time
    - event_chi2 for the $\chi_{\nu}^2$ as given in the file.
    - event_stations has been calculated for us from event_mask. That's helpful!
    - Each event has been tagged with a unique event_id
These variables are easy to access by name - for example, we can get a nice display of
Step4:  Notice that this xarray DataArray variable is not only the data values, but some other metadata, such as the units and the standard variable name from the Climate and Forecast Metadata Conventions. We can access those attributes if we want them
Step5:  Simple plotting and filtering
xarray has some nice, built-in plotting functions for quickly looking at data. We'll make even nicer plots in a bit!
Step6:  We have lots of data in the file; let's grab the first 10000 points to make the plotting faster.
Step7:  Looks pretty noisy. Let's try again, but filter to lower chi2 and greater event_contributing_stations.
Suggested activity
<ASSISTANT_TASK:>
Python Code:
# We could tediously build a list …
# filenames = ['/data/Houston/realtime-tracer/LYLOUT_200524_210000_0600.dat.gz',]
# Instead, let's read a couple hours at the same time.
import sys, glob
filenames = glob.glob('/data/Houston/130619/LYLOUT_130619_2[0-1]*.dat.gz')
for filename in filenames:
    print(filename)
import glob
import numpy as np
import datetime
import xarray as xr
import pyproj as proj4
from pyxlma.lmalib.io import read as lma_read
lma_data, starttime = lma_read.dataset(filenames)
Explanation: Reading and plotting LMA data
In a previous module we learned about LMA VHF source files - their structure, and principles of data quality control using the station mask and $\chi_{\nu}^2$.
We'll use the xlma-python package to read, process, and visualize the data in these files.
We'll start by importing some libraries and reading a file. In fact, we can read more than Replace the filename below with a file of interest to you, if you want.
End of explanation
# Should match what we expect from the filenames
print(starttime)
print(type(starttime))
Explanation: Investigating the pyxlma data structure
We just read in the data from a list of filenames, and got back an lma_data variable, and a starttime. What's in these variables?
starttime is a Python datetime object, with helps us do math operations on times instead of multiplying by 60 and 24 ourselves.
End of explanation
print(lma_data)
Explanation: lma_data is an xarray object. If we print it, we see that it looks much like a NetCDF file, with dimensions and variables that allow us to store whole arrays of data and give them names. xarray is the best way to look at NetCDF data in Python.
(If you'd like to learn more about xarray and NetCDF, you can check out Unidata's lessons on xarray and NetCDF.)
End of explanation
print(type(lma_data.event_longitude))
lma_data.event_longitude
Explanation: There are a few things to notice above.
- The dimensions tell us how many events and stations we had in our data file.
- Some other header info about the network location and the lma_analysis command is included.
- There are varaibles with dimension number_of_stations that give the data from the station data tables.
- There are varaibles with dimension number_of_events, including many we'd expect.
    - Event location and time
    - event_chi2 for the $\chi_{\nu}^2$ as given in the file.
    - event_stations has been calculated for us from event_mask. That's helpful!
    - Each event has been tagged with a unique event_id
These variables are easy to access by name - for example, we can get a nice display of
End of explanation
print(lma_data.event_longitude.attrs['standard_name'])
Explanation: Notice that this xarray DataArray variable is not only the data values, but some other metadata, such as the units and the standard variable name from the Climate and Forecast Metadata Conventions. We can access those attributes if we want them:
End of explanation
%matplotlib widget
import matplotlib.pyplot as plt
lma_data.plot.scatter?
Explanation: Simple plotting and filtering
xarray has some nice, built-in plotting functions for quickly looking at data. We'll make even nicer plots in a bit!
End of explanation
fig, axes = plt.subplots(1,1,figsize=(10,10))
count_subset = {'number_of_events':slice(0,10000)}
art = lma_data[count_subset].plot.scatter('event_longitude', 'event_latitude', ax=axes,
                      s=4, marker='s', #hue='event_time', 
                     )
Explanation: We have lots of data in the file; let's grab the first 10000 points to make the plotting faster.
End of explanation
fig, axes = plt.subplots(1,1,figsize=(10,10))
count_subset = {'number_of_events':slice(0,10000)}
station_filter = (lma_data.event_stations >= 6)
chi_filter = (lma_data.event_chi2 <= 1.0)
filter_subset = {'number_of_events':(chi_filter & station_filter)}
# note that we first filter on all the data select 10000 points, and then on that dataset we further filter 
art = lma_data[filter_subset][count_subset].plot.scatter('event_longitude', 'event_latitude', ax=axes,
                      s=4, marker='s', #hue='event_time', 
                     )
Explanation: Looks pretty noisy. Let's try again, but filter to lower chi2 and greater event_contributing_stations.
Suggested activity: adjust the station and $\chi^2$ criteria to see how the noisiness changes. As noted in earlier lessons, some experimentation to choose the level of filtering is a very common process when working with LMA data.
You can also filter on event latitude and longitude in the same way to further zoom in, or plot other combinations of data variables.
End of explanation
<END_TASK>