schain Commit - r1145:e6529314b4b0 · Jicamarca Repository

Change AER to Geodetic coordinates in PX plots

jespinoza -

r1145:e6529314b4b0

parent child

Context file:

r1145:e6529314b4b0

Collapse all files

schainpy/model/graphics/jroplot_data.py +35 -20

             CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis', 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm', 'spectral')]
+            EARTH_RADIUS = 6.3710e3
             def ll2xy(lat1, lon1, lat2, lon2):
                 p = 0.017453292519943295
                 theta = -theta + numpy.pi/2
                 return r*numpy.cos(theta), r*numpy.sin(theta)
+            def km2deg(km):
+                '''
+                Convert distance in km to degrees
+                '''
+                return numpy.rad2deg(km/EARTH_RADIUS)
             def figpause(interval):
                 backend = plt.rcParams['backend']
                 if backend in matplotlib.rcsetup.interactive_bk:
                     ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
                     ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
-                    Y = numpy.array([5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000])
+                    Y = numpy.array([1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000])
                     i = 1 if numpy.where(abs(ymax-ymin) <= Y)[0][0] < 0 else numpy.where(abs(ymax-ymin) <= Y)[0][0]
                     ystep = Y[i] / 5.
                     for n, ax in enumerate(self.axes):
                         if ax.firsttime:
                     self.nrows = 1
                     self.width = 9
                     self.height = 8
+                    self.mode = self.data.meta['mode']
                     if self.channels is not None:
                         self.nplots = len(self.channels)
                         self.nrows = len(self.channels)
                         self.nplots = self.data.shape(self.CODE)[0]
                         self.nrows = self.nplots
                         self.channels = range(self.nplots)
-                    if self.data.meta['mode'] == 'E':
+                    if self.mode == 'E':
-                        self.xlabel = 'Zonal Distance (km)'
+                        self.xlabel = 'Longitude'
-                        self.ylabel = 'Meridional Distance (km)'
+                        self.ylabel = 'Latitude'
                     else:
                         self.xlabel = 'Range (km)'
                         self.ylabel = 'Height (km)'
                     self.bgcolor = 'white'
                     self.cb_labels = self.data.meta['units']
-                    # self.polar = True
+                    self.lat = self.data.meta['latitude']
+                    self.lon = self.data.meta['longitude']
-                def plot(self):
+                    self.xmin, self.xmax = float(km2deg(-50) + self.lon), float(km2deg(50) + self.lon)
+                    self.ymin, self.ymax = float(km2deg(-50) + self.lat), float(km2deg(50) + self.lat)
+                    log.error(type(self.ymin))
+                    #print km2deg(-50) + self.lon, km2deg(50) + self.lon
+                    #print km2deg(-50) + self.lat, km2deg(50) + self.lat
+                    #  self.polar = True
+                def plot(self):
                     for n, ax in enumerate(self.axes):
                         data = self.data['param'][self.channels[n]]
                         zeniths = numpy.linspace(0, self.data.meta['max_range'], data.shape[1])
-                        if self.data.meta['mode'] == 'E':
+                        if self.mode == 'E':
                             azimuths = -numpy.radians(self.data.heights)+numpy.pi/2
                             r, theta = numpy.meshgrid(zeniths, azimuths)
                             x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
+                            x = km2deg(x) + self.lon
+                            y = km2deg(y) + self.lat
                         else:
                             azimuths = numpy.radians(self.data.heights)
                             r, theta = numpy.meshgrid(zeniths, azimuths)
                                                     vmax=self.zmax,
                                                     cmap=self.cmaps[n])
-                        if self.data.meta['mode'] == 'A':
+                        if self.mode == 'A':
                             continue
                         f = open('/home/jespinoza/workspace/schain_scripts/distrito.csv')
-                        lat1 = -11.96
-                        lon1 = -76.54
                         for line in f:
                             label, lon, lat = [s.strip() for s in line.split(',') if s]
                             lat = float(lat)
                             lon = float(lon)
-                            x, y = ll2xy(lat1, lon1, lat, lon)
+                            ax.plot(lon, lat, '.b', ms=2)
-                            ax.plot(x, y, '.b', ms=2)
+                            ax.text(lon, lat, label.decode('utf8'), ha='center', va='bottom', size='8', color='black')
-                            ax.text(x, y, label.decode('utf8'), ha='center', va='bottom', size='8', color='black')
-                    if self.data.meta['mode'] == 'E':
+                    if self.mode == 'E':
                         title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
-                        label = 'E{:d}'.format(int(self.data.meta['elevation']))
+                        label = 'E{:02d}'.format(int(self.data.meta['elevation']))
                     else:
                         title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
-                        label = 'A{:d}'.format(int(self.data.meta['azimuth']))
+                        label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
                     self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
                     self.titles = ['{} {}'.format(self.data.parameters[x], title) for x in self.channels]

schainpy/model/proc/pxproc_parameters.py +1 -1

                     else:
                         self.dataOut.heightList = self.dataOut.data['Elevation']
-                    attrs = ['units', 'elevation', 'azimuth', 'max_range']
+                    attrs = ['units', 'elevation', 'azimuth', 'max_range', 'latitude', 'longitude']
                     meta = {}
                     for attr in attrs:

General Comments 0

Write
Preview

You need to be logged in to leave comments. Login now

No TODOs yet

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages