-
Notifications
You must be signed in to change notification settings - Fork 0
/
hail_section.py
79 lines (76 loc) · 4.65 KB
/
hail_section.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
import numpy as np
from shapely import geometry
from shapely.ops import transform
from metpy.units import atleast_1d, check_units, concatenate, units
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from pyproj import Geod
def hail_objects(hailc,REF_Hail2,ax,f,time_start,month,d_beg,h_beg,min_beg,sec_beg,d_end,h_end,min_end,sec_end,rlons,rlats,max_lons_c,max_lats_c,proj):
#Inputs,
#REF_Hail2: REFmasked masked where Zdr and CC greater than 1.0
#hailc: Contour of REF_Hail2 where reflectivity greater than 50.0 dBz
#ax: Subplot object to be built on with each contour
#f: Placefile, edited throughout the program
#time_start: Radar file date and time of scan
#month: Month of case, supplied by user
#d_beg,h_beg,min_beg,sec_beg,d_end,h_end,min_end,sec_end: Day, hour, minute, second of the beginning and end of a scan
#rlons,rlats: Full volume geographic coordinates, longitude and latitude respectively
#max_lons_c,max_lats_c: Centroid coordinates of storm objects
#proj: Projection of Earth's surface to be used for accurate area and distance calculations
hail_areas = []
hail_centroid_lon = []
hail_centroid_lat = []
hail_storm_lon = []
hail_storm_lat = []
if np.max(REF_Hail2) > 50.0:
for level in hailc.collections:
for contour_poly in level.get_paths():
for n_contour,contour in enumerate(contour_poly.to_polygons()):
contour_a = np.asarray(contour[:])
xa = contour_a[:,0]
ya = contour_a[:,1]
polygon_new = geometry.Polygon([(i[0], i[1]) for i in zip(xa,ya)])
if n_contour == 0:
polygon = polygon_new
else:
polygon = polygon.difference(polygon_new)
pr_area = (transform(proj, polygon).area * units('m^2')).to('km^2')
boundary = np.asarray(polygon.boundary.xy)
polypath = Path(boundary.transpose())
coord_map = np.vstack((rlons[0,:,:].flatten(), rlats[0,:,:].flatten())).T
#Create an Mx2 array listing all the coordinates in field
mask_hail = polypath.contains_points(coord_map).reshape(rlons[0,:,:].shape)
if pr_area > 2 * units('km^2'):
g = Geod(ellps='sphere')
dist_hail = np.zeros((np.asarray(max_lons_c).shape[0]))
for i in range(dist_hail.shape[0]):
distance_hail = g.inv(polygon.centroid.x, polygon.centroid.y,
max_lons_c[i], max_lats_c[i])
dist_hail[i] = distance_hail[2]/1000.
if np.min(np.asarray(dist_hail)) < 15.0:
hail_path = polypath
hail_areas.append((pr_area))
hail_centroid_lon.append((polygon.centroid.x))
hail_centroid_lat.append((polygon.centroid.y))
hail_storm_lon.append((max_lons_c[np.where(dist_hail == np.min(dist_hail))[0][0]]))
hail_storm_lat.append((max_lats_c[np.where(dist_hail == np.min(dist_hail))[0][0]]))
patch = PathPatch(polypath, facecolor='gold', alpha=.7, edgecolor = 'gold', linewidth = 4)
ax.add_patch(patch)
#Add polygon to placefile
f.write('TimeRange: '+str(time_start.year)+'-'+str(month)+'-'+str(d_beg)+'T'+str(h_beg)+':'+str(min_beg)+':'+str(sec_beg)+'Z '+str(time_start.year)+'-'+str(month)+'-'+str(d_end)+'T'+str(h_end)+':'+str(min_end)+':'+str(sec_end)+'Z')
f.write('\n')
f.write("Color: 245 242 066 \n")
f.write('Line: 3, 0, "Hail Core Outline" \n')
for i in range(len(hail_path.vertices)):
f.write("%.5f" %(hail_path.vertices[i][1]))
f.write(", ")
f.write("%.5f" %(hail_path.vertices[i][0]))
f.write('\n')
f.write("End: \n \n")
#Returning Variables,
#hail_areas: Hail core area
#hail_centroid_lon,hail_centroid_lat: Hail core centroid coordinates
#hail_storm_lon,hail_storm_lat: Storm object centroids associated with the hail core
#ax: Subplot object to be built on with each contour
#f: Placefile, edited throughout the program
return hail_areas,hail_centroid_lon,hail_centroid_lat,hail_storm_lon,hail_storm_lat,ax,f