Latitude, longitude and altitude of radar gates

[jjhelmus]

A method should be added to the Radar class for determining the latitude, longitude and altitude of a particular gate or all gates in the radar volume.

This functionality was suggested by Timothy Lang.

[nguy]

I just posted in the User's group about this. I agree that this functionality would be quite useful and already exists in essence for grid objects. The question may be naming and location. Currently a radar object has latitude and longitude which correspond to radar latitude and longitude. I propose that the 'long_name' an 'standard_name' be changed to 'Radar XX' for clarity. I don't believe this would have any consequences.

I've always been in favor a radar_position dictionary like instrument_parameters that would hold this as well as altitude, altitude_agl, others(?). The reason I bring this up is that adding a latitude and longitude field will conflict with the current top level parameters. An alternative would be to add it to fields, but that doesn't seem like the correct place. So where to add these variables?

[nguy]

On a related, but side note would it make more sense to have these functions (as well as `radar_to_cartesian_coords) in pyart.util?

The code below should give what we are looking for by combining a few different functions in Py_ART:

def xy_to_latlon(radar):
    rng, az = np.meshgrid(radar.range['data'], radar.azimuth['data'])
    rng, ele = np.meshgrid(radar.range['data'], radar.elevation['data'])
    theta_e = ele * np.pi / 180.0       # elevation angle in radians.
    theta_a = az * np.pi / 180.0        # azimuth angle in radians.
    Re = 6371.0 * 1000.0 * 4.0 / 3.0     # effective radius of earth in meters.
    r = rng * 1000.0                    # distances to gates in meters.

    z = (r ** 2 + Re ** 2 + 2.0 * r * Re * np.sin(theta_e)) ** 0.5 - Re
    s = Re * np.arcsin(r * np.cos(theta_e) / (Re + z))  # arc length in m.
    x = s * np.sin(theta_a)
    y = s * np.cos(theta_a)

    c = np.sqrt(x*x + y*y) / r
    phi_0 = radar.latitude['data'] * np.pi / 180
    azi = np.arctan2(y, x)  # from east to north

    lat = np.arcsin(np.cos(c) * np.sin(phi_0) +
                    np.sin(azi) * np.sin(c) * np.cos(phi_0)) * 180 / np.pi
    lon = (np.arctan2(np.cos(azi) * np.sin(c), np.cos(c) * np.cos(phi_0) -
           np.sin(azi) * np.sin(c) * np.sin(phi_0)) * 180 /
            np.pi + radar.longitude['data'])
    lon = np.fmod(lon + 180, 360) - 180

    lat_axis = {
        'data':  lat,
        'long_name': 'Latitude for points in Cartesian system',
        'axis': 'YX',
        'units': 'degree_N',
        'standard_name': 'latitude',
    }

    lon_axis = {
        'data': lon,
        'long_name': 'Longitude for points in Cartesian system',
        'axis': 'YX',
        'units': 'degree_E',
        'standard_name': 'longitude',
    }
    return lat_axis, lon_axis

[jjhelmus]

I just posted in the User's group about this. I agree that this functionality would be quite useful and already exists in essence for grid objects. The question may be naming and location. Currently a radar object has latitude and longitude which correspond to radar latitude and longitude. I propose that the 'long_name' an 'standard_name' be changed to 'Radar XX' for clarity. I don't believe this would have any consequences.

The only difficulty here is that the long_name and standard_name values are taken from the configuration file which defines generic metadata. At the moment the latitude and longitude keys are only used to specify Radar latitudes and longitudes but as they currently are defined they could be used for grid latitudes and longitudes. Not that we couldn't make the change, just will need to keep this in mind.

I've always been in favor a radar_position dictionary like instrument_parameters that would hold this as well as altitude, altitude_agl, others(?). The reason I bring this up is that adding a latitude and longitude field will conflict with the current top level parameters. An alternative would be to add it to fields, but that doesn't seem like the correct place. So where to add these variables?

Nearly all the dictionaries off the Radar class match one to one with those specified in the CfRadial format. Introducing a radar_position dictionary breaks this mapping so I would be in favor of a different methods of storing the latitude and longitudes of the radar gates, maybe gate_latitude and gate_longitude attribute dictionaries which are initially None? These could also be defined as properties like rays_per_sweep which is calculated on access. This may be non-ideal since calculating the locations of all gates can be a bit computationally demanding and we may not want to do it more than once. Caching the property is an option but then we need to worry about how to invalidate the cache and if it should be done automatically.

Providing a function which returns arrays of gate latitudes and longitudes avoids this but puts the burden of keeping track of these values on the user.

[jjhelmus]

Also should mention the code in @nguy's comment works fine if you want to determine the gate locations using an azimuthal equidistant projection (at least I think that is what is implemented). As discussed in #294 a better solution would allow the users to select a projects supported by basemap/pyproj if installed and fall back or default to az. equidistant if basemap/pyproj is not installed. This takes more work but would give a more complete solution to the problem.

[jjhelmus]

Looked into cached properties and there is a nice package for it. Invalidating the cache is as simple as deleting the attribute which I think the Py-ART users would understand. Also having an invalid cached version of these arrays would only occur if the user changes the latitude, longitude, azimuth, elevation or range attributes, which I would think uncommon...and if they do they probably know what they are doing and would not access the gate_latitude and gate_longitude attributes until after modifying the pointing data.

Assuming this is how we choose to implement this feature, access to the gate latitudes and longitudes is through gate_latitude and gate_longitude cached properties off the Radar class, then the only question is how the projection should be specified? It can't be a parameter in a function/method call since these are properties so would a attribute be sufficient?

[nguy]

Dang cached properties are cool. I like the gate_latitude and gate_longitude labels, I completely forgot about cfradial for the moment. I agree that deleting and changing should not be too difficult to figure out and if you are changing these items you are likely advanced or used to dealing with the mess that can be radar data. :)

Thanks for clarifying the code above. It is just for AED projection, which is pretty good for a single radar domain generally. I would say that a basemap instance could be an optional input set to None in the function and would fall back to the above code if None. Otherwise, the projection information could be obtained from the basemap instance for calculations. I think that is similar to what you do in the grid function now.

[jjhelmus]

I'll try to put together a PR that implements some of this later this week. I have something in mind now but need to see it in code.

[nguy]

Sounds good!

[jjhelmus]

Thinking about this a bit more I realized the gate latitudes, longitudes, altitudes and Cartesian locations (x, y, z) can all be represented as Radar attributes with the data loaded on demand using a LazyLoadDict. This gives us low memory overhead (the full array is only calculated when the data is accessed), a standard interface (same as all other Radar attributes), and the logic is already available. Will be starting to work on this today and hopefully will have some initial code by the end of the day or tomorrow.

[nguy]

Nice @jjhelmus. This could be very helpful for some stuff @tjlang and I are doing in AWOT right now!

[anshul12345]

please post if there is any script available to convert radar spherical coordinates to geographic coordinates.

[scollis]

Please read the docs or post to the email list. Do not use the issue tracker to ask questions