ENH: Corrected Height Gates to be height above sea level

pyart/util/columnsect.py

@@ -179,13 +179,13 @@ def get_field_location(radar, latitude, longitude):
     # call the get_column_ray function
     ray = get_column_rays(radar, azim)
 
-    # Determine the gates for the rays
+    # Determine the center of each gate for the subsetted rays.
     (rhi_x,
      rhi_y,
      rhi_z) = antenna_vectors_to_cartesian(radar.range['data'],
                                            radar.azimuth['data'][ray],
                                            radar.elevation['data'][ray],
-                                           edges=True)
+                                           edges=False)
     # Calculate distance from the x,y coordinates to target
     rhidis = np.sqrt((rhi_x**2) + (rhi_y**2)) * np.sign(rhi_z)
     for i in range(len(ray)):
@@ -196,7 +196,9 @@ def get_field_location(radar, latitude, longitude):
             else:
                 moment[key].append(radar.fields[key]
                                    ['data'][ray[i], tar_gate])
-        zgate.append(rhi_z[i, tar_gate])
+        # Add radar elevation to height gates
+        # to define height as center of each gate above sea level
+        zgate.append(rhi_z[i, tar_gate] + radar.altitude['data'][0])
 
     # Determine the time at the center of each ray within the column
     # Define the start of the radar volume as a numpy datetime object for xr
@@ -214,31 +216,57 @@ def get_field_location(radar, latitude, longitude):
     # Create a blank list to hold the xarray DataArrays
     ds_container = []
     da_meta = ['units', 'standard_name', 'long_name', 'valid_max',
-               'valid_min']
+               'valid_min', 'coordinates']
     # Convert the moment dictionary to xarray DataArray.
     # Apply radar object meta data to DataArray attribute
     for key in moment:
-        da = xr.DataArray(moment[key],
-                          coords=dict(height=zgate),
-                          name=key, dims=['height'])
-        for tag in da_meta:
-            if tag in radar.fields[key]:
-                da.attrs[tag] = radar.fields[key][tag]
-        # Append to ds container
-        ds_container.append(da.to_dataset(name=key))
+        if key != 'height':
+            da = xr.DataArray(moment[key],
+                              coords=dict(height=zgate),
+                              name=key, dims=['height'])
+            for tag in da_meta:
+                if tag in radar.fields[key]:
+                    da.attrs[tag] = radar.fields[key][tag]
+            # Append to ds container
+            ds_container.append(da.to_dataset(name=key))
+
+    # Add additional DataArrays 'base_time' and 'time_offset'
+    # if not present within the radar object.
+    da_base = xr.DataArray(base_time, name='base_time')
+    da_offset = xr.DataArray(combined_time,
+                             coords=dict(height=zgate),
+                             name='time_offset', dims=['height'])
+    ds_container.append(da_base.to_dataset(name='base_time'))
+    ds_container.append(da_offset.to_dataset(name='time_offset'))
 
     # Create a xarray DataSet from the DataArrays
     column = xr.merge(ds_container)
-    # Add the time as an additional coordinate
-    column = column.assign_coords(coords=dict(time=combined_time))
-    if len(radar.time['units'].split(' ')) > 3:
-        column.attrs['scan_start_time'] = base_time
-    else:
-        column.attrs['scan_start_time'] = base_time
-    column.attrs['distance'] = (str(np.around(dis / 1000., 3)) + ' km')
+
+    # Assign Attributes for the Height and Times
+    height_des = ('Height Above Sea Level [in meters] for the Center of Each'
+                  + ' Radar Gate Above the Target Location')
+    column.height.attrs.update(long_name='Height of Radar Beam',
+                               units='m',
+                               standard_name='height',
+                               description=height_des)
+
+    column.base_time.attrs.update(long_name='UTC Reference Time',
+                                  units='seconds')
+
+    time_long = 'Time in Seconds Since Volume Start'
+    time_des = ("Time in Seconds Since Volume Start that Cooresponds"
+               + " to the Center of Each Height Gate"
+               + " Above the Target Location")
+    column.time_offset.attrs.update(long_name=time_long,
+                                    units='seconds',
+                                    description=time_des)
+
+    # Assign Global Attributes to the DataSet
+    column.attrs['distance_from_radar'] = (str(np.around(dis / 1000., 3))
+                                          + ' km')
     column.attrs['azimuth'] = (str(np.around(azim, 3)) + ' degrees')
-    column.attrs['latitude'] = (str(latitude) + ' degrees')
-    column.attrs['longitude'] = (str(longitude) + ' degrees')
+    column.attrs['latitude_of_location'] = (str(latitude) + ' degrees')
+    column.attrs['longitude_of_location'] = (str(longitude) + ' degrees')
     return column
 
 

tests/util/test_columnsect.py

@@ -31,8 +31,9 @@ def test_sphere_distance():
 def test_get_field_location():
     # test to make sure column above location is pulled correctly
     column = pyart.util.columnsect.get_field_location(radar, 47.97, -122.20)
+
     # check to make sure z-gate is pulled correctly.
-    test_height = abs(column.height[0] - 367)
+    test_height = abs(column.height[0] - 564)
     assert (test_height < 0.001)
 
     # check to make sure reflectivity value is minimum