Export Reflectivity Based on Latitude and Longitude

[A syafiuddin]

Dear all,

I am working with Radar data. How I can extract the reflectivity data based on latitude (x) and longitude (y) using Py-ART? Thank you.

Best Regards
Achmad Syafiuddin

[Jonathan Helmus]

Achmad,

    In the latest development version of Py-ART the latitude and longitude of each radar gate is accessible in the gate_latitude and gate_longitude attributes of the Radar class.  These should allow you to map between rays and gates and latitude and longitude. 

Cheers,

    - Jonathan Helmus

--
You received this message because you are subscribed to the Google Groups "Py-ART Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pyart-users...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[A syafiuddin]

Dear Dr Jonathan
Thank you very much for your suggestion. I constructed a program to extract the x,y, and z data succesfully in meters. I sent you in your personal email the code and data I used in the present study. It is possible for me to extract the data in kilometers? Another question is how to extract the data in the specific x, and y at z = 2000 meter? For example I want to extract the reflectivity data at x = 4.67874 N and y = 101.4849 E only. My purpose is to compare the data to the station data.
Thank you very much.

Sincerely Yours,
Achmad Syafiuddin

[Jonathan Helmus]

All,

Achmad and I have been discussing this question off list but I though I would share a script  finds the nearest radar gate to a given location as well as the reflectivity interpolated from nearby gatesaround a given location.  Such features may be helpful to other Py-ART users:


from __future__ import print_function
import numpy as np
import pyart
import pyproj


def find_x_y_displacement(radar, longitude, latitude):
    """ Return the x and y displacement (in meters) from a radar location. """
    # longitude and latitude in degrees
    lat_0 = radar.latitude['data'][0]
    lon_0 = radar.longitude['data'][0]
    proj = pyproj.Proj(proj='aeqd', lon_0=lon_0, lat_0=lat_0)
    return proj(longitude, latitude)


def find_nearest_gate(radar, longitude, latitude, altitude):
    """ Return the indices of the nearest gate to a given point. """
    # longitude and latitude in degrees, altitude in meters
    gate_x = radar.gate_x['data']
    gate_y = radar.gate_y['data']
    gate_z = radar.gate_z['data']

    x_disp, y_disp = find_x_y_displacement(radar, longitude, latitude)
    distances = np.sqrt(
        (gate_x-x_disp)**2. + (gate_y-y_disp)**2. + (gate_z-altitude)**2.)
    return np.unravel_index(distances.argmin(), distances.shape)


def interpolate_single_point(radar, longitude, latitude, altitude):
    """ Interpolate a single grid point at a given location. """
    x_disp, y_disp = find_x_y_displacement(radar, longitude, latitude)
    grid = pyart.map.grid_from_radars(
        (radar,),
        gridding_algo='map_gates_to_grid',
        grid_shape=(1, 1, 1),
        grid_limits=((altitude, 0), (y_disp, 25000.0), (x_disp, 25000.0)),
        fields=['reflectivity'])
    return grid


# read in the file
filename = 'RADAR_FILE'
radar = pyart.io.read(filename)

latitude = 80.0    # latitude (in degrees) to find reflectivity at or near
longitude = 120.0  # longitude (in degrees) to find reflectivity at or near
altitude = 1000       # altitude (in meters) to find reflectivity at or near

# find nearest gate
ray, gate = find_nearest_gate(radar, longitude, latitude, altitude)
gate_latitude = radar.gate_latitude['data'][ray, gate]
gate_longitude = radar.gate_longitude['data'][ray, gate]
gate_altitude = radar.gate_altitude['data'][ray, gate]
gate_reflectivity = radar.fields['reflectivity']['data'][ray, gate]
print("-----------------------")
print("Nearest gate:")
print("Latitude:", gate_latitude)
print("Longitude:", gate_longitude)
print("Altitude:", gate_altitude)
print("Reflectivity:", gate_reflectivity)

# interpolate around gate
grid = interpolate_single_point(radar, longitude, latitude, altitude)
pyart.io.add_2d_latlon_axis(grid)
pixel_latitude = grid.axes['latitude']['data'][0, 0]
pixel_longitude = grid.axes['longitude']['data'][0, 0]
pixel_altitude = grid.axes['z_disp']['data'][0]
pixel_reflectivity = grid.fields['reflectivity']['data'][0, 0]

print("-----------------------")
print("Interpolated pixel:")
print("Latitude:", pixel_latitude)
print("Longitude:", pixel_longitude)
print("Altitude:", pixel_altitude)
print("Reflectivity:", pixel_reflectivity)

Cheers,

    - Jonathan Helmus

[Diar]

Hi Sayafiuddin,

Do you mind posting your  final code here?

thanks

Diar

[James Cho]

I'm very late here, but I am just getting into radar processing with pyart, wanted to check my understanding of this code. This line in particular:

grid_limits=((altitude, 0), (y_disp, 25000.0), (x_disp, 25000.0)),

means we are finding a single point to average the reflectivity values a 25x25 km square, correct? And this could contain one or multiple gates?


Thanks,

James Cho

[Scott Collis]

James: are you at AMS radar next week?

This Is an ooooold thread. Py-ART has moved on since then.

Please let us know your particular use case and some one on the list can help

--

Scott Collis

[James Cho]

Hey Scott, thanks for the fast reply. I'm not at AMS Radar. I'm working with NEXRAD Level 2 data, with the same goal as OP: getting reflectivity at a lat/long. My code is based off this: https://github.com/openradar/AMS_radar_in_the_cloud/blob/master/notebooks/Putting_it_all_together.ipynb First time working with radar data and coming from a mechanical engineering background.

[Scott Collis]

Check this out.. here is how I extracted the gate closest to a distrometer 

https://github.com/scollis/AGU_2016/blob/master/scripts/multi_proccessing.py#L219

--

Scott Collis

[James Cho]

Thanks. Is there a quick way to add interpolation to that, or should I use the previously posted code, with the Grid object?

[Scott Collis]

Sorry dude,

You are on your own from here 

-sent from a mobile device-