# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# +
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import pyart
from wsrlib import *

# requirements:
#  - pyart
#  - boto3
#  - pyproj


# +
#name = 'KBUF20160107_123931_V06'
#name = 'KBUF20160806_101141_V06'

def image2xy(x, y, r=0, dim=600, rmax=150000):
    '''
    Convert from image coordinates to (x,y) coordinates offset from radar
    '''
    
    x0 = y0 = dim/2.0 # origin
    x =  (x - x0)*2*rmax/dim
    y = -(y - y0)*2*rmax/dim
    r = r*2*rmax/dim
    
    return x, y, r


# +
import pandas as pd
roosts = pd.read_csv('roost_labels_KBUF2016.csv')

roosts = roosts.loc[roosts['label']=='swallow-roost']
roosts['tot_refl'] = 0

groups = roosts.groupby('filename')

n = len(groups)
for j, (scan, group) in enumerate(groups):

    if (j < 133 or j > 138):
        continue
    
    print('Processing %s (%03d/%03d)' % (scan, j, n))
    
    for i,row in group.iterrows():
        #print(i)
        #print(i, row['track_id'], row['from_sunrise'])
        roosts.loc[i,'tot_refl'] = i*10

display(roosts)

# -

print(len(roosts['track_id']))

# +
# File name and roost location
name = 'KBUF20160826_101937_V06' # 1907,KBUF20160826_101937_V06,-14,0.657,235.294551,278.52,26.456692
x0, y0, r = image2xy(235.294551,278.52,26.456692)

name = 'KBUF20160826_103825_V06' # 1913,KBUF20160826_103825_V06,5,1.024,64.915009,390.60,52.056629
x0, y0, r = image2xy(64.915009,390.60,52.056629)

# name = 'KBUF20160808_101440_V06' # 1612,KBUF20160808_101440_V06,1,0.863,45.938196,389.39,45.936024
# x0, y0, r = image2xy(45.938196,389.39,45.936024)
print(x0,y0,r)

radar = read_s3(name)

# +
# Plot reflectivity
disp = pyart.graph.RadarDisplay(radar)
rmax = 150 # km
lim = (-rmax, rmax)
disp.plot('reflectivity', 0, title='NEXRAD Reflectivity',
             vmin=-5, vmax=35, colorbar_label='dBZ')
disp.set_limits(xlim=lim, ylim=lim)

# Plot roost
rect = patches.Rectangle(((x0-r)/1000,(y0-r)/1000),2*r/1000,2*r/1000,linewidth=1,edgecolor='r',facecolor='none')
plt.gca().add_patch(rect)

plt.show()
# -

data, y, x, elev = radar2mat(radar, coords='cartesian', dim=600, r_max=150000)
print(data['reflectivity'].shape)

# +
rng, az, elev, dbz = get_volumes(radar, coords='antenna')

# Convert to reflectivity and set nan --> zero
refl, _ = z_to_refl(idb(dbz))
refl[np.isnan(refl)] = 0

# +
## Pulse volume coordinate conversions

# Get range along ground and height above ground
[ground_rng, height] = slant2ground(rng, elev)

# Convert compass bearing to mathematical angle
theta = cmp2pol(az)

# Convert to X,Y coordinates of each pulse volume
[x, y] = pol2cart(theta, ground_rng)

# +
# Convert roost location to lat / lon
[angle, dist_from_radar] = cart2pol(x0, y0)
bearing = pol2cmp(angle)

from pyproj import Geod
geodesic = Geod(ellps='WGS84')
roost_lon, roost_lat, _ = geodesic.fwd(radar.longitude['data'][0], 
                                       radar.latitude['data'][0] , 
                                       bearing, 
                                       dist_from_radar)

#roost_lon = do_alias(roost_lon, 180)
print('lat,lon=%.4f,%.4f' % (roost_lat, roost_lon))

# +
DIST = np.sqrt( (x-x0)**2 + (y-y0)**2 ) # distance of each sample volume to center of roost

edges = np.arange(0, 1.5*r, 500)
nbins = len(edges)-1

density = np.full(nbins, np.nan)

for i in range(nbins):
    inds = (DIST >= edges[i]) & (DIST < edges[i+1]) # & (elev < 0.55)
    density[i] = np.sum(refl[inds])
# -

bin_start = edges[:-1]
plt.plot(bin_start, density, 'o')
plt.show()

# +
inds = DIST < r
plt.hist(DIST[inds], weights=refl[inds], bins=40)
plt.show()

inds = (DIST < r) & (height < 4000)
plt.hist(height[inds], weights=refl[inds], bins=40)
plt.show()
# -





StarcoderdataJupyterScriptsDedupFiltered