##// END OF EJS Templates
schain
RSS
Use function elaz to latlon instead of wradlib
Juan C. Espinoza -
r1530:69c85e5f5aa8
parent child
Show More
Side by Side Unified
Context file:
r1530:69c85e5f5aa8 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,5 +1,6
1 1 import os
2 2 import datetime
3 import warnings
3 4 import numpy
4 5 from mpl_toolkits.axisartist.grid_finder import FixedLocator, DictFormatter
5 6
@@ -7,11 +8,150 from schainpy.model.graphics.jroplot_base import Plot, plt
7 8 from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
8 9 from schainpy.utils import log
9 10
10 import wradlib.georef as georef
11 11
12 12 EARTH_RADIUS = 6.3710e3
13 13
14 14
15 def antenna_to_cartesian(ranges, azimuths, elevations):
16 """
17 Return Cartesian coordinates from antenna coordinates.
18
19 Parameters
20 ----------
21 ranges : array
22 Distances to the center of the radar gates (bins) in kilometers.
23 azimuths : array
24 Azimuth angle of the radar in degrees.
25 elevations : array
26 Elevation angle of the radar in degrees.
27
28 Returns
29 -------
30 x, y, z : array
31 Cartesian coordinates in meters from the radar.
32
33 Notes
34 -----
35 The calculation for Cartesian coordinate is adapted from equations
36 2.28(b) and 2.28(c) of Doviak and Zrnic [1]_ assuming a
37 standard atmosphere (4/3 Earth's radius model).
38
39 .. math::
40
41 z = \\sqrt{r^2+R^2+2*r*R*sin(\\theta_e)} - R
42
43 s = R * arcsin(\\frac{r*cos(\\theta_e)}{R+z})
44
45 x = s * sin(\\theta_a)
46
47 y = s * cos(\\theta_a)
48
49 Where r is the distance from the radar to the center of the gate,
50 :math:`\\theta_a` is the azimuth angle, :math:`\\theta_e` is the
51 elevation angle, s is the arc length, and R is the effective radius
52 of the earth, taken to be 4/3 the mean radius of earth (6371 km).
53
54 References
55 ----------
56 .. [1] Doviak and Zrnic, Doppler Radar and Weather Observations, Second
57 Edition, 1993, p. 21.
58
59 """
60 theta_e = numpy.deg2rad(elevations) # elevation angle in radians.
61 theta_a = numpy.deg2rad(azimuths) # azimuth angle in radians.
62 R = 6371.0 * 1000.0 * 4.0 / 3.0 # effective radius of earth in meters.
63 r = ranges * 1000.0 # distances to gates in meters.
64
65 z = (r ** 2 + R ** 2 + 2.0 * r * R * numpy.sin(theta_e)) ** 0.5 - R
66 s = R * numpy.arcsin(r * numpy.cos(theta_e) / (R + z)) # arc length in m.
67 x = s * numpy.sin(theta_a)
68 y = s * numpy.cos(theta_a)
69 return x, y, z
70
71 def cartesian_to_geographic_aeqd(x, y, lon_0, lat_0, R=EARTH_RADIUS):
72 """
73 Azimuthal equidistant Cartesian to geographic coordinate transform.
74
75 Transform a set of Cartesian/Cartographic coordinates (x, y) to
76 geographic coordinate system (lat, lon) using a azimuthal equidistant
77 map projection [1]_.
78
79 .. math::
80
81 lat = \\arcsin(\\cos(c) * \\sin(lat_0) +
82 (y * \\sin(c) * \\cos(lat_0) / \\rho))
83
84 lon = lon_0 + \\arctan2(
85 x * \\sin(c),
86 \\rho * \\cos(lat_0) * \\cos(c) - y * \\sin(lat_0) * \\sin(c))
87
88 \\rho = \\sqrt(x^2 + y^2)
89
90 c = \\rho / R
91
92 Where x, y are the Cartesian position from the center of projection;
93 lat, lon the corresponding latitude and longitude; lat_0, lon_0 are the
94 latitude and longitude of the center of the projection; R is the radius of
95 the earth (defaults to ~6371 km). lon is adjusted to be between -180 and
96 180.
97
98 Parameters
99 ----------
100 x, y : array-like
101 Cartesian coordinates in the same units as R, typically meters.
102 lon_0, lat_0 : float
103 Longitude and latitude, in degrees, of the center of the projection.
104 R : float, optional
105 Earth radius in the same units as x and y. The default value is in
106 units of meters.
107
108 Returns
109 -------
110 lon, lat : array
111 Longitude and latitude of Cartesian coordinates in degrees.
112
113 References
114 ----------
115 .. [1] Snyder, J. P. Map Projections--A Working Manual. U. S. Geological
116 Survey Professional Paper 1395, 1987, pp. 191-202.
117
118 """
119 x = numpy.atleast_1d(numpy.asarray(x))
120 y = numpy.atleast_1d(numpy.asarray(y))
121
122 lat_0_rad = numpy.deg2rad(lat_0)
123 lon_0_rad = numpy.deg2rad(lon_0)
124
125 rho = numpy.sqrt(x*x + y*y)
126 c = rho / R
127
128 with warnings.catch_warnings():
129 # division by zero may occur here but is properly addressed below so
130 # the warnings can be ignored
131 warnings.simplefilter("ignore", RuntimeWarning)
132 lat_rad = numpy.arcsin(numpy.cos(c) * numpy.sin(lat_0_rad) +
133 y * numpy.sin(c) * numpy.cos(lat_0_rad) / rho)
134 lat_deg = numpy.rad2deg(lat_rad)
135 # fix cases where the distance from the center of the projection is zero
136 lat_deg[rho == 0] = lat_0
137
138 x1 = x * numpy.sin(c)
139 x2 = rho*numpy.cos(lat_0_rad)*numpy.cos(c) - y*numpy.sin(lat_0_rad)*numpy.sin(c)
140 lon_rad = lon_0_rad + numpy.arctan2(x1, x2)
141 lon_deg = numpy.rad2deg(lon_rad)
142 # Longitudes should be from -180 to 180 degrees
143 lon_deg[lon_deg > 180] -= 360.
144 lon_deg[lon_deg < -180] += 360.
145
146 return lon_deg, lat_deg
147
148 def antenna_to_geographic(ranges, azimuths, elevations, site):
149
150 x, y, z = antenna_to_cartesian(numpy.array(ranges), numpy.array(azimuths), numpy.array(elevations))
151 lon, lat = cartesian_to_geographic_aeqd(x, y, site[0], site[1], R=6370997.)
152
153 return lon, lat
154
15 155 def ll2xy(lat1, lon1, lat2, lon2):
16 156
17 157 p = 0.017453292519943295
@@ -457,15 +597,21 class WeatherParamsPlot(Plot):
457 597 data['ele'] = dataOut.data_ele
458 598 data['mode_op'] = dataOut.mode_op
459 599 var = data['data'].flatten()
460 r = numpy.tile(data['r'], data['data'].shape[0]).reshape(data['data'].shape)*1000
461 lla = georef.spherical_to_proj(r, data['azi'], data['ele'], (-75.295893, -12.040436, 3379.2147))
600 r = numpy.tile(data['r'], data['data'].shape[0])
601 az = numpy.repeat(data['azi'], data['data'].shape[1])
602 el = numpy.repeat(data['ele'], data['data'].shape[1])
603
604 # lla = georef.spherical_to_proj(r, data['azi'], data['ele'], (-75.295893, -12.040436, 3379.2147))
605
606 latlon = antenna_to_geographic(r, az, el, (-75.295893, -12.040436))
607
462 608 if self.mask:
463 meta['lat'] = lla[:,:,1].flatten()[var.mask==False]
464 meta['lon'] = lla[:,:,0].flatten()[var.mask==False]
609 meta['lat'] = latlon[1][var.mask==False]
610 meta['lon'] = latlon[0][var.mask==False]
465 611 data['var'] = numpy.array([var[var.mask==False]])
466 612 else:
467 meta['lat'] = lla[:,:,1].flatten()
468 meta['lon'] = lla[:,:,0].flatten()
613 meta['lat'] = latlon[1]
614 meta['lon'] = latlon[0]
469 615 data['var'] = numpy.array([var])
470 616
471 617 return data, meta
Show file before | Show file after | Hide comments
@@ -4146,7 +4146,6 class PedestalInformation(Operation):
4146 4146 ok = True
4147 4147 break
4148 4148 except Exception as e:
4149 print(e)
4150 4149 log.warning('Waiting {}s for position file to be ready...'.format(self.delay), self.name)
4151 4150 time.sleep(self.delay)
4152 4151 continue
@@ -4168,8 +4167,7 class PedestalInformation(Operation):
4168 4167 flag_mode = None
4169 4168 azi = self.fp['Data']['azi_pos'][:]
4170 4169 ele = self.fp['Data']['ele_pos'][:]
4171 #print("az: ",az)
4172 #exit(1)
4170
4173 4171 while True:
4174 4172 if start+sample_max > numpy.shape(ele)[0]:
4175 4173 print("CANNOT KNOW IF MODE IS PPI OR RHI, ANALIZE NEXT FILE")
@@ -4219,7 +4217,7 class PedestalInformation(Operation):
4219 4217 else:
4220 4218 return numpy.nan, numpy.nan, numpy.nan
4221 4219
4222 def setup(self, dataOut, path, conf, samples, interval, mode):
4220 def setup(self, dataOut, path, conf, samples, interval, mode, online):
4223 4221
4224 4222 self.path = path
4225 4223 self.conf = conf
@@ -4235,15 +4233,26 class PedestalInformation(Operation):
4235 4233 log.error('No position files found in {}'.format(path), self.name)
4236 4234 raise IOError('No position files found in {}'.format(path))
4237 4235 else:
4236 if self.online:
4237 self.filename = filelist[-1]
4238 self.utcfile = int(self.filename.split('/')[-1][4:14])
4239 log.log('Opening file: {}'.format(self.filename), self.name)
4240 for i in range(self.nTries):
4241 try:
4242 self.fp = h5py.File(self.filename, 'r')
4243 except:
4244 log.warning('Waiting {}s for position file to be ready...'.format(self.delay), self.name)
4245 time.sleep(self.delay)
4246 else:
4238 4247 self.filename = filelist[0]
4239 4248 self.utcfile = int(self.filename.split('/')[-1][4:14])
4240 4249 log.log('Opening file: {}'.format(self.filename), self.name)
4241 4250 self.fp = h5py.File(self.filename, 'r')
4242 4251
4243 def run(self, dataOut, path, conf=None, samples=1500, interval=0.04, az_offset=0, time_offset=0, mode=None):
4252 def run(self, dataOut, path, conf=None, samples=1500, interval=0.04, time_offset=0, mode=None, online=False):
4244 4253
4245 4254 if not self.isConfig:
4246 self.setup(dataOut, path, conf, samples, interval, mode)
4255 self.setup(dataOut, path, conf, samples, interval, mode, online)
4247 4256 self.isConfig = True
4248 4257
4249 4258 self.utctime = dataOut.utctime + time_offset
@@ -4256,10 +4265,8 class PedestalInformation(Operation):
4256 4265 dataOut.flagNoData = True
4257 4266 return dataOut
4258 4267
4259 dataOut.azimuth = az + az_offset
4260 if dataOut.azimuth < 0:
4261 dataOut.azimuth += 360
4262 dataOut.elevation = el
4268 dataOut.azimuth = round(az, 2)
4269 dataOut.elevation = round(el, 2)
4263 4270 dataOut.mode_op = scan
4264 4271
4265 4272 return dataOut
@@ -4307,7 +4314,6 class Block360(Operation):
4307 4314 Add a profile to he __buffer and increase in one the __profiel Index
4308 4315 '''
4309 4316 tmp= getattr(data, attr)
4310
4311 4317 self.__buffer.append(tmp)
4312 4318 self.__buffer2.append(data.azimuth)
4313 4319 self.__buffer3.append(data.elevation)
General Comments 0
You need to be logged in to leave comments. Login now