Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Updates paths and reorganize. #2

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Updates paths and reorganize. #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
7f0a72f
Updates paths and reorganize.
nguy Dec 10, 2015
1a14db5
Updated datatypes to be consistently numpy.datetime64.
nguy Dec 11, 2015
e2e8883
Vispy 0.5 required
nguy Dec 11, 2015
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
57 changes: 23 additions & 34 deletions radar/collection.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,10 +16,9 @@ def data_for_ray_slice(radar, ray_sl, fieldnames=None):
corresponding to the fieldnames. If fieldname is None,
no data will be returned. Fieldnames should be a sequence of field
names. A data dictionary is returned with a key for each of the field names.

ray_sl is typically found from the radar's sweep index record.
For some sweep id N, this is:
ray_sl = slice( radar.sweep_start_ray_index['data'][N],
ray_sl = slice( radar.sweep_start_ray_index['data'][N],
radar.sweep_end_ray_index['data'][N] + 1)
"""
r = radar.range['data']
Expand All @@ -35,19 +34,18 @@ def data_for_ray_slice(radar, ray_sl, fieldnames=None):
data = None
print("Couldn't find {1} data in {0}".format(radar,fieldname))
data_dict[fieldname] = data

return r,az,el,t,data_dict

def iter_sweep_data(radar, fieldnames):
for swp_start, swp_end in zip(radar.sweep_start_ray_index['data'],
for swp_start, swp_end in zip(radar.sweep_start_ray_index['data'],
radar.sweep_end_ray_index['data']):
print swp_start, swp_end
ray_sl = slice(swp_start,swp_end+1)
yield data_for_ray_slice(radar,ray_sl,fieldnames=fieldnames)

class RadarFileCollection(object):
"""

"""
def __init__(self, filenames):
self.filenames=filenames
Expand All @@ -61,15 +59,15 @@ def __init__(self, filenames):
self.sweep_table = pandas.DataFrame([v for v in self._iter_sweep_index_data()],
columns = ('filename', 'sweep_slice', 'start', 'end', 'mode', 'angle')
)

def _iter_sweep_time_range(self, fname):
""" Yields (sweep_slice, min_t, max_t)
for the given radar filename. Times are datetime objects.
"""
radar = self.radars[fname]
datetimes = self.times[fname]
for swp_start, swp_end in zip(radar.sweep_start_ray_index['data'],
radar.sweep_end_ray_index['data']):
for swp_start, swp_end in zip(radar.sweep_start_ray_index['data'],
radar.sweep_end_ray_index['data']):
ray_sl = slice(swp_start,swp_end+1)
t = datetimes[ray_sl]
yield ray_sl, min(t), max(t)
Expand All @@ -80,49 +78,42 @@ def _iter_sweep_index_data(self):
mode = radar.scan_type.lower()
vcp = diagnose_vcp(radar)
for (swp_sl, swp_ta, swp_tb), angle in zip(self._iter_sweep_time_range(fname), vcp):
yield fname, swp_sl, swp_ta, swp_tb, mode, angle

yield fname, swp_sl, np.datetime64(swp_ta, 'ns'), np.datetime64(swp_tb, 'ns'), mode, angle


def sweep_for_time_range(self, t0, t1, overlap_idx=0):
""" Given some time range t0, t1 return the closest sweep to that time.
t0 and t1 are datetime objects

If there is only one sweep that overlaps, return it.

overlap_idx = 0: (default) returns the first sweep with overlap
overlap_idx = -1: returns the last sweep with overlap
So, the default behavior is to return the first sweep to have any overlap.

Returns filename, sweep_slice.

Unimplemented
-------------
If there are two sweeps that overlap, one could:
return first or last, or the one with the most overlap

With more than two sweeps, the first and last will have partial overlap
return first or last partial overlap
return the first or last that fully overlaps
return the sweep with most overlap
choosing the first or last of those if there are several of the same length



t0 t1
| |
... -------------- -------------- ------------------ -------------- ...
ta3 tb3 ta4 tb4 ta5 tb5 ta6 tb6

conditions for some overlap
conditions for some overlap
t1 > ta4 so t1 - ta4 > 0
t0 < tb4 so t0 - tb4 < 0

"""

#target = pandas.DataFrame([(t0, t1),], columns=('start', 'end'))

overlap = ((t0 - self.sweep_table['end']) < 0) & ((t1 - self.sweep_table['start']) > 0 )
## print("t0/t1 type = {0}/{1}".format(np.dtype(t0), np.dtype(t1)))
## print("sweep_table['end'] type = {0}".format(np.dtype(self.sweep_table['end'])))
## print("sweep_table['start'] type = {0}".format(np.dtype(self.sweep_table['start'])))
dt_zero = np.timedelta64(0,'ns')
overlap = ((t0 - self.sweep_table['end'].values) < dt_zero) & ((t1 - self.sweep_table['start'].values) > dt_zero )
sweeps = self.sweep_table[overlap]
if len(sweeps) > 0:
selection = np.zeros(len(sweeps), dtype=bool)
Expand All @@ -132,12 +123,11 @@ def sweep_for_time_range(self, t0, t1, overlap_idx=0):
return sweep_info['filename'].values[0], sweep_info['sweep_slice'].values[0]
else:
return None, None

def sweep_data_for_time_range(self, t0, t1, fieldnames=None, **kwargs):
""" return r,az,el,t,data corresponding to the fieldnames. If fieldname is None,
no data will be returned. Fieldnames should be a sequence of field
names. A data dictionary is returned with a key for each of the field names.

t0, t1 and extra kwargs are passed to self.sweep_for_time_range
"""
filename, ray_slice = self.sweep_for_time_range(t0, t1, **kwargs)
Expand All @@ -146,17 +136,16 @@ def sweep_data_for_time_range(self, t0, t1, fieldnames=None, **kwargs):
return data_for_ray_slice(radar, ray_slice, fieldnames=fieldnames)
else:
return None, None, None, None, None

def _time_range_for_file(self, fname):
""" Called once by init to set up frame lookup tables and yield
the frame start times. _frame_lookup goes from
""" Called once by init to set up frame lookup tables and yield
the frame start times. _frame_lookup goes from
datetime->(nc file, frame index)"""
# datetimes_from_radar returns the time of each ray to the nearest second
times = pyart.util.datetime_utils.datetimes_from_radar(self.radars[fname])
return min(times), max(times)

def _all_times(self):
for f in self.filenames:
for tmin,tmax in self._time_range_for_file(f):
yield tmin,tmax

yield tmin, tmax
63 changes: 33 additions & 30 deletions radar/quadmesh_geometry.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,84 +4,85 @@
import numpy as np
pi = np.pi

import vispy
import vispy.app
# from vispy.scene.visuals.modular_mesh import ModularMesh
from vispy.scene.visuals import Mesh
from vispy.geometry import MeshData
from vispy.scene import STTransform, MatrixTransform
from vispy.scene import STTransform , MatrixTransform

def synthetic_field(x,y):
# Trapp and Doswell (2000, Mon. Weather Rev.) analytic input field
# Trapp and Doswell (2000, Mon. Weather Rev.) analytic input field
z0 = 0.5
amplitude = 0.5
Lx = Ly = 30.0
x_shift = 0#10.0 # shift the field, leave the radar at 0,0
y_shift = 15/4.0 #-5.0
x_wavenumber = 2
y_wavenumber = x_wavenumber
D = z0+amplitude * (np.cos(2*pi*x_wavenumber*(x-x_shift)/Lx) *
D = z0+amplitude * (np.cos(2*pi*x_wavenumber*(x-x_shift)/Lx) *
np.sin(2*pi*y_wavenumber*(y-y_shift)/Ly) )
return D

# create polar mesh data
def radar_example_data():
""" Create a cone of regularly gridded data in spherical
""" Create a cone of regularly gridded data in spherical
(range, azimuth, elevation) coordinates. The data values
have an undulating pattern tied to the cartesian x,y location.

This is a quadrilateral 2D mesh surface, but demonstrative of
a more general case where the quadrilaterals' corners
a more general case where the quadrilaterals' corners
must each be specified manually. This is a common need in the case
of, for instance, weather radar data.
Returns x,y,z locations of the M x N vertices and

Returns x,y,z locations of the M x N vertices and
(M-1) x (N-1) data value arrays.
"""
radar_x0 = 0.0
radar_y0 = 0.0
radar_z0 = 0.0

el = np.radians(30.0)

all_r = np.arange(5, 35, 0.75)
all_az_deg = np.arange(0,361, 2)
all_az_deg = np.arange(0,361, 2)
all_az = np.radians(all_az_deg)

# Observation locations
r_edge, az_edge = np.meshgrid(all_r, all_az)
# >>> print r.shape
# >>> print r.shape
# 360, 200 = (naz, nr)
r = (( r_edge[1:, 1:] + r_edge[:-1, :-1]) / 2.0)
az = ((az_edge[1:, 1:] + az_edge[:-1, :-1]) / 2.0)

x_radar_edge = r_edge*np.cos(az_edge) + radar_x0
y_radar_edge = r_edge*np.sin(az_edge) + radar_y0
z_radar_edge = r_edge*np.sin(el) + radar_z0

x_radar = r*np.cos(az) + radar_x0
y_radar = r*np.sin(az) + radar_y0
z_radar = r*np.sin(el) + radar_z0

# values at sampled locations
all_d = synthetic_field(x_radar, y_radar)

return (x_radar_edge, y_radar_edge, z_radar_edge, all_d)

def mesh_from_quads(x,y,z):
""" x, y, z are M x N arrays giving the edge locations for
a quadrilateral mesh of Nq = (M-1) x (N-1) implied quad faces.
After conversion to triangles, this is Nf=2*Nq faces.
Nv = M x N.

returns
vertices : ndarray, shape (Nv, 3) - Vertex coordinates.
vertices : ndarray, shape (Nv, 3) - Vertex coordinates.
faces : ndarray, shape (Nf, 3) - Indices into the vertex array.

Along each row, the triangles' diagonals are oriented in
the same direction. Each new row of triangles is specified from
left-to-right.
An alternate face specification (1) would go back and forth,

An alternate face specification (1) would go back and forth,
alternating the triangles' diagonal direction. It would be
more amenable GL_TRIANGLE_STRIP.
(1) http://dan.lecocq.us/wordpress/2009/12/25/triangle-strip-for-grids-a-construction/
Expand All @@ -92,35 +93,37 @@ def mesh_from_quads(x,y,z):
Nf = 2*Nq
verts = np.empty((Nv, 3), dtype='f4')
faces = np.empty((Nf, 3), dtype='i4')

verts[:,0] = x.flat
verts[:,1] = y.flat
verts[:,2] = z.flat

q_range = np.arange(Nq, dtype='i4')

# When the index along M changes, need to skip ahead to avoid connecting
# the edges together.
# Should look like (0,0,0,1,1,1,2,2,2,3,3,3) for N-1=3 and M-1=4
adjust_start = (np.arange(M-1)[:,None]*np.ones(N-1)[None,:]).flatten()
adjust_start = (np.arange(M-1, dtype=np.int32)[:,None]*np.ones(N-1, dtype=np.int32)[None,:]).flatten()
# adjust_start = (np.arange(M-1, dtype='i4')[:,None]*np.ones(N-1, dtype='i4')[None,:]).flatten()
q_range+=adjust_start

# even (0,2,4,...) triangles
faces[0::2, 0] = q_range + 1
faces[0::2, 1] = q_range
faces[0::2, 1] = q_range
faces[0::2, 2] = q_range + N
# odd (1,3,5,...) trianges
faces[1::2, 0] = q_range + N
faces[1::2, 1] = q_range + N + 1
faces[1::2, 2] = q_range + 1

return verts, faces



class Canvas(vispy.scene.SceneCanvas):
def __init__(self):
self.rotation = MatrixTransform()
## self.rotation = MatrixTransform()
self.rotation = vispy.scene.transforms.MatrixTransform()

# Generate some data to work with
x,y,z,d = radar_example_data()
Expand All @@ -140,7 +143,7 @@ def __init__(self):
mesh = Mesh(meshdata=mdata)
mesh.transform = vispy.scene.transforms.MatrixTransform()
mesh.transform.scale([1./10, 1./10, 1./10])

vispy.scene.SceneCanvas.__init__(self, keys='interactive')

self.size = (800, 800)
Expand Down
Loading