receiver_functions.py 文件源码

def find_pierce_coor(self,plot='False'):
      import geopy
      from geopy.distance import VincentyDistance

      '''
      given an instance of the receiver function class this function
      returns latitude and longitude of all receiver side pierce points
      of Pds in a given depth range (the default range is 50 - 800 km)
      NOTE:
      be careful that ses3d typically uses colatitude, while
      this function returns latitude '''

      depth_range = np.arange(50,800,5)        #set range of pierce points

      #geodetic info
      bearing     = self.az
      lon_s = self.ses3d_seismogram.sy
      lat_s = 90.0-self.ses3d_seismogram.sx
      lon_r = self.ses3d_seismogram.ry
      lat_r = 90.0-self.ses3d_seismogram.rx
      origin      = geopy.Point(lat_s, lon_s)

      #find how far away the pierce point is
      model  = TauPyModel(model='pyrolite_5km')

      for i in range(0,len(depth_range)):
         phase = 'P'+str(depth_range[i])+'s'
         pierce = model.get_pierce_points(self.eq_depth,self.delta_deg,phase_list=[phase])
         points = pierce[0].pierce
         for j in range(0,len(points)):
            if points[j]['depth'] == depth_range[i] and points[j]['dist']*(180.0/np.pi) > 25.0:
               prc_dist = points[j]['dist']*(180.0/np.pi)
               d_km = prc_dist * ((2*np.pi*6371.0/360.0))
               destination = VincentyDistance(kilometers=d_km).destination(origin,bearing)
               lat = destination[0]
               lon = destination[1]
               value = 0
               row = {'depth':depth_range[i],'dist':prc_dist,'lat':lat,'lon':lon,'value':value}
               self.pierce_dict.append(row)

      if plot=='True':
         m = Basemap(projection='hammer',lon_0=0)
         m.drawmapboundary()
         m.drawcoastlines()
         m.drawgreatcircle(lon_s,lat_s,lon_r,lat_r,linewidth=1,color='b',alpha=0.5)

         for i in range(len(self.pierce_dict)):
            x,y = m(self.pierce_dict[i]['lon'],self.pierce_dict[i]['lat'])
            m.scatter(x,y,5,marker='o',color='r')
         plt.show()