Python Basemap.gcpoints方法代码示例

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import gcpoints [as 别名]
def ray_density(lat1, lon1, lat2, lon2,
                dt=1, gr_x=360, gr_y=180, 
                npts=180, projection='robin', 
                ray_coverage=False):
    '''
    Create the DATA array which contains the
    info for ray density
    '''
    global long_0

    mymap = Basemap(projection=projection, lon_0=long_0, lat_0=0)
    #npts=max(gr_x, gr_y)
    # grd[2]: longitude
    # grd[3]: latitude
    grd = mymap.makegrid(gr_x, gr_y, returnxy=True)

    lons, lats = mymap.gcpoints(lon1, lat1, lon2, lat2, npts)
    dist = locations2degrees(lat1, lon1, lat2, lon2)
    bap = int((dist - 97.0)*npts/dist)/2
    midlon = len(lons)/2
    midlat = len(lats)/2
    lons = lons[midlon-bap:midlon+1+bap]
    lats = lats[midlat-bap:midlat+1+bap]
    data = np.zeros([len(grd[2]), len(grd[3])])
    for i in range(len(lons)):
        xi, yi = point_finder(lons[i], lats[i], grd)
        # first one is latitude and second longitude
        try:
            #data[yi][xi] = dt/float(dist-97.0)
            data[yi][xi] += dt/len(lons)
        except Exception, e:
            print e

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import gcpoints [as 别名]
def sat_error():
  boxsize=6000
  lat,lon=52.52437,13.41053 # centre Berlin map
  m = Basemap(projection='ortho',lon_0=lon,lat_0=lat,resolution='l',area_thresh='l',llcrnrx=-boxsize,llcrnry=-boxsize,urcrnrx=+boxsize,urcrnry=+boxsize)
  
  x0,y0=m(lo_sat[0],la_sat[0])
  x1,y1=m(lo_sat[1],la_sat[1])
  dog=sqrt((x1-x0)**2.+(y1-y0)**2.)
  print x0+0.5*(x1-x0),y0+0.5*(y1-y0)
  print distance_on_earth(la_sat[0],lo_sat[0],la_sat[1],lo_sat[1])
  
  x,y=m.gcpoints(lo_sat[0],la_sat[0], lo_sat[1],la_sat[1], 10000)
  
  log,lag=m(x[5000],y[5000],inverse=True)
  lop,lap=m(x0+0.5*(x1-x0),y0+0.5*(y1-y0),inverse=True)
  print lag,log,lap,lop
  return distance_on_earth(lag,log,lap,lop)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import gcpoints [as 别名]
def ray_density(lat1, lon1, lat2, lon2, dt=1, gr_x=360, gr_y=180, npts=180, projection='robin', ray_coverage=False):
    """
    Create the DATA array which contains the info for ray density

    Procedure:
    1. make a grid based on the inputs (grd)
    grd: lon, lat, x, y
    2. find the great circle points:
    note that lon , lat are actually x and y!
    3. calculate the distance and find the middle point
    4. subtracting 97 degrees from the distance and find all the points on that section
    5. data ---> zero array with x*y elements

    """

    global long_0

    mymap = Basemap(projection=projection, lon_0=long_0, lat_0=0)
    #npts=max(gr_x, gr_y)
    # grd[2]: longitude
    # grd[3]: latitude
    grd = mymap.makegrid(gr_x, gr_y, returnxy=True)

    lons, lats = mymap.gcpoints(lon1, lat1, lon2, lat2, npts)
    dist = locations2degrees(lat1, lon1, lat2, lon2)

    # npts points on dist...how many on (dist-97)!: (dist-97)*npts/dist....but we also need to make it half!
    bap = int((dist - 97.0)*npts/dist)/2
    midlon = len(lons)/2
    midlat = len(lats)/2
    lons = lons[midlon-bap:midlon+1+bap]
    lats = lats[midlat-bap:midlat+1+bap]

    data = np.zeros([len(grd[2]), len(grd[3])])
    if not len(lons) == len(lats):
        sys.exit('ERROR: Lengths longitudes and latitudes are not the same! %s and %s' % (len(lons), len(lats)))

    for i in range(len(lons)):
        xi, yi = point_finder(lons[i], lats[i], grd)
        # first one is latitude and second longitude
        try:
            #data[yi][xi] = dt/float(dist-97.0)
            data[yi][xi] += dt/len(lons)
        except Exception, e:
            print '\nException: %s' % e

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import gcpoints [as 别名]
    m.drawparallels(np.arange(-90.,120.,30.))
    m.drawmeridians(np.arange(0.,420.,60.))
    m.drawmapboundary()

    gc_xmid_points = []
    gc_ymid_points = []
    gc_dtmid_points = []
    for i in range(len(passed_staev)):
        sys.stdout.write('\r')
        sys.stdout.write("[%-100s] %d%%" % ('='*int(100.*(i+1)/len(passed_staev)),
                                                100.*(i+1)/len(passed_staev)))
        sys.stdout.flush()
        #m.drawgreatcircle(passed_staev[i][5],passed_staev[i][4],
        #            passed_staev[i][1],passed_staev[i][0], alpha = 0.1, 
        #            color=passed_staev[i][-1])
        gc_arc = m.gcpoints(passed_staev[i][5], passed_staev[i][4],
                    passed_staev[i][1], passed_staev[i][0], int(divisions))
        gc_xmid_points.append(gc_arc[0][int(divisions/2)])
        gc_ymid_points.append(gc_arc[1][int(divisions/2)])
        gc_dtmid_points.append(passed_staev[i][2])
        #m.scatter(gc_arc[0][50], gc_arc[1][50], color=passed_staev[i][-1])
        #m.scatter(gc_arc[0][50], gc_arc[1][50], color=passed_staev[i][2], vmin=-6, vmax=6.0)

    map_evlats = []
    map_evlons = []
    print '\nError:' 
    for i in range(len(proc_ev_ls)):
        try:
            evx, evy = m(evlons[i], evlats[i])
            map_evlons.append(evx)
            map_evlats.append(evy)
        except Exception, e:

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import gcpoints [as 别名]
    end_latt=math.radians(end_latt)
    d_latt = end_latt - start_latt  
    d_long = end_long - start_long
    a = math.sin(d_latt/2)**2 + math.cos(start_latt) * math.cos(end_latt) * math.sin(d_long/2)**2  
    c = 2 * math.asin(math.sqrt(a))  
    return 6371 * c *1000

longslist=[0,90,180]
latslist=[65,65,65]
dpslist=[5,10,5]
dmslist=[5,5,10]
print m(0,60)
print m(0,70)
print m(-5,65)
print m(5,65)
print m.gcpoints(0,60,0,70,10)

m.drawgreatcircle(0,60,0,70,100)
m.drawgreatcircle(-5,65,5,65,100)

m.drawgreatcircle(0,60,-5,65,100)
m.drawgreatcircle(0,60,5,65,100)

m.drawgreatcircle(-5,65,0,70,100)
m.drawgreatcircle(0,70,5,65,100)
'''
# draw ellipses
for n in range(0,len(longslist)):
    xy= m(longslist[n],latslist[n])
    w = haversine(longslist[n]-dmslist[n],latslist[n],longslist[n]+dmslist[n],latslist[n])
    h = haversine(longslist[n],latslist[n]-dpslist[n],longslist[n],latslist[n]+dpslist[n])

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import gcpoints [as 别名]
def main():
  """Finds location of maximum probability in a 2D space defined by three seperate PDFs."""
  
  # plotting
  fig=plt.figure(11,figsize=(6,6))
  plt.clf()
  plt.subplots_adjust(left=0.01, bottom=0.01, right=0.99, top=0.99, wspace=0., hspace=0.)
  ax = plt.subplot(1,1,1)
  xticklabels = getp(gca(), 'xticklabels')
  yticklabels = getp(gca(), 'yticklabels')
  
  # create map and set size
  boxsize=6000
  lat, lon=52.52437, 13.41053 # centre berlin map
  m = Basemap(projection='ortho',lon_0=lon,lat_0=lat,resolution='l',area_thresh='l',llcrnrx=-boxsize,llcrnry=-boxsize,urcrnrx=+boxsize,urcrnry=+boxsize)
  
  # read and plot street map / water areas from shapefile
  m.readshapefile('berlin_germany/berlin_germany_osm_roads_gen1', 'Streets',color='#DDDDDD')
  m.readshapefile('berlin_germany/berlin_germany_osm_waterareas', 'Water',color='#add8e6')
  
  # compute and plot points on great circle
  x,y=m.gcpoints(lo_sat[0],la_sat[0], lo_sat[1],la_sat[1], 100)
  m.plot(x,y,'r-',zorder=7)
  
  # plot spree list of coords
  x,y=m(lo_spree,la_spree)
  m.plot(x,y,'b.',markersize=2)
  
  # plot brandenburg gate
  gx,gy=m(lo_bgate,la_bgate) # coord transform
  m.plot(gx,gy,'bx')
  
  # create search grid
  bestY,bestX=m(13.4497839044, 52.5142566528 ) # grid centre
  minX=bestX-8000
  maxX=bestX+8000
  minY=bestY-8000
  maxY=bestY+8000
  
  rangeX=np.linspace(minX,maxX,nbins)
  rangeY=np.linspace(minY,maxY,nbins)
  npoints=None
  
  meshX, meshY = meshgrid(rangeX, rangeY)
  mzz=np.zeros(shape=(len(rangeY),len(rangeX)))
  maxmzz=0
  
  # compute probabilities for grid and location with max. prob.
  for ii in range(0, len(rangeY)):
    s = str(int(float(ii+1)/nbins*100.)) + '% computed'
    print s,
    sys.stdout.flush()
    backspace(len(s))
    for jj in range(0, len(rangeX)):
      mY,mX=m(rangeY[ii],rangeX[jj],inverse=True)
      prob_bgate,prob_spree,prob_sat,x_bgate,x_spree,x_sat=z_prob(m,mX,mY,npoints)
      mzz[ii][jj]=prob_bgate*prob_spree*prob_sat
      if mzz[ii][jj]>maxmzz:
        maxmzz=mzz[ii][jj]
        maxmX=mX
        maxmY=mY
        best_x_bgate=x_bgate
        best_x_spree=x_spree
        best_x_sat=x_sat
  
  # plot location with max. prob.
  loc_x,loc_y= m(maxmY,maxmX)
  m.plot(loc_x,loc_y,'r.',zorder=8)
  
  # output
  print "\nDone!"
  print "+++ search log +++"
  print "max prob.:", maxmzz, " found at: (lat,lon)=(",maxmX,",",maxmY,")"
  print "with dist. to spree: ", best_x_spree
  print "with dist. to bgate: ", best_x_bgate
  print "with dist. to sat: ", best_x_sat
  print "resolution in m: (dx,dy)=",(maxX-minX)/nbins,(maxY-minY)/nbins #
  
  # plotting: color map and contours
  logcrange=arange(-5,-2.1,0.1)
  crange=10**logcrange
  cs = m.contourf(meshY,meshX,mzz,levels=crange,cmap=cm.Blues,zorder=5, alpha=0.3)
  
  # plotting legend
  prop = matplotlib.font_manager.FontProperties(size=10)
  b, = plot([10000,20000],'bx', linewidth=1)
  c, = plot([10000,20000],'b.', linewidth=1)
  d, = plot([10000,20000],'r-', linewidth=1)
  e, = plot([10000,20000],'r.', linewidth=1)
  legend([b,c,d,e,], ['Brandenburg Gate','Spree (file)','projected satellite path','coordinates with max. probability'],loc='upper left', ncol=1, shadow=False, fancybox=True, numpoints=1, prop=prop,labelspacing=-0.0)
  
  # plotting scale
  shift=6700
  sx=6368777.0109+shift
  sy=6370097.04832+shift
  m.plot([sx,sx+1000],[sy,sy],'k-',linewidth=5)
  s1lo,s1la=m(sx+1000,sy,inverse=True)
  s2lo,s2la=m(sx,sy,inverse=True)
  text(sx,sy-400,r' $1$ $km$')
  
#.........这里部分代码省略.........

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import gcpoints [as 别名]
with open("test_sms.xml", "r") as filetoberead:
    body = filetoberead.readlines()[3:-1] #Read in all the texts
    
freq = {} #Create a database and associate each text with a frequency
for text in body:
    num = sms.get_area(text)
    if num in freq:
        freq[num] += 1
    else:
        freq[num] = 1

codes = sorted(freq, key=freq.get)
results = []
for each in codes:
    results.append(freq[each])
    
print codes[1], results[1]
#results.reverse() #Reversed the values so the color would end up darker for more usage; lighter for lighter usage
    
for i, num in enumerate(codes):
    scale = (ceil(np.log(results[i])))/10 #Creates a log scale for the relative frequencies
    rgb = colorswap(0.3,0.3,scale) #alter the HSV scale to RGB to keep a single color and adjust the tone
    try:
        long2,lat2 = area_codes[int(num)][1:]
        z = n.gcpoints(lat2,long2, lat1,long1,10)
        n.plot(z[0],z[1],color=rgb,linewidth = 1.3)
    except: pass #I have tried to code for many of the exceptions, but i was getting strange area codes which do not exist. Not sure what to do with them....

plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import gcpoints [as 别名]
def select_sta():
    
    """
    Select required stations
    """
    
    global input
    
    map_proj = Basemap(projection='cyl', llcrnrlat=-90,urcrnrlat=90,\
                llcrnrlon=-180,urcrnrlon=180, resolution='c')
    ev_file = open(os.path.join(os.getcwd(), 'quake_req.txt'), 'r')
    ev_add = ev_file.read().split('\n')[:-1]
    select = open(os.path.join(os.getcwd(), input['file'] + '.dat'), 'w')
    select.close()
    
    for k in range(0, len(ev_add)):
        
        '''
        select = open(os.path.join(os.getcwd(), \
                input['file'] + '-' + ev_add[k].split('/')[-1] + \
                '.dat'), 'w')
        '''
        
        (quake_d, quake_t) = read_quake(ev_add[k])
        
        list_sta = glob.glob(os.path.join(ev_add[k], 'BH', \
                                input['identity']))
        
        for i in range(0, len(list_sta)):
            
            try:
                
                st = read(list_sta[i])
                print '***************************************'
                print str(i) + '/' + str(len(list_sta)) + ' -- ' + \
                        str(k) + '/' + str(len(ev_add))
                print list_sta[i].split('/')[-1]
                
                info_sac = st[0].stats['sac']
                
                if input['all_sta'] == None:
                    
                    dist = locations2degrees(lat1 = quake_d['lat'], \
                            long1 = quake_d['lon'], lat2 = info_sac['stla'], \
                            long2 = info_sac['stlo'])
                    tt = getTravelTimes(delta=dist, depth=quake_d['dp'], \
                                            model=input['model'])
                    
                    for m in range(0, len(tt)):
                        
                        if tt[m]['phase_name'] in input['phase']:

                            try:
                                print '--------------------'
                                print list_sta[i].split('/')[-1] + ' has ' + \
                                        tt[m]['phase_name'] + ' phase'
                                
                                if input['freq'] != None:
                                    st[0].decimate(int(round(\
                                        st[0].stats['sampling_rate'])/input['freq']), \
                                        no_filter=False)
                                    
                                    if st[0].stats['sampling_rate'] != input['freq']:
                                        print list_sta[i].split('/')[-1]
                                        print st[0].stats['sampling_rate']
                                        print '------------------------------------------'
                                '''
                                np_evt = round((events[0]['datetime'] - st[0].stats['starttime'])*st[0].stats['sampling_rate'])
                                np_pha = np_evt + round(tt[m]['time']*st[0].stats['sampling_rate'])
                                
                                select = open(Address_events + '/' + events[l]['event_id'] + '/IRIS/info/' + name_select, 'a')
                                '''
                                if tt[m]['phase_name'] != 'Pdiff':
                                    lat_1 = str(quake_d['lat'])
                                    lon_1 = str(quake_d['lon'])
                                    lat_2 = str(info_sac['stla'])
                                    lon_2 = str(info_sac['stlo'])
                                elif tt[m]['phase_name'] == 'Pdiff':
                                    dist_limit = 97.0
                                    num_gcp = 1000
                                    gcp = map_proj.gcpoints(quake_d['lon'], \
                                            quake_d['lat'], info_sac['stlo'], \
                                            info_sac['stla'], num_gcp)

                                    if dist >= dist_limit:
                                        diff_dist = dist - dist_limit
                                                        
                                        req_gcp = diff_dist*(float(num_gcp)/dist)
                                        req_gcp = round(req_gcp)/2

                                        mid_p = len(gcp[0])/2
                                        #before = int(mid_p - req_gcp)
                                        #after = int(mid_p + req_gcp)
                                        before = mid_p - int(2.0 * len(gcp[0])/dist)
                                        after = mid_p + int(2.0 * len(gcp[0])/dist)
                                        
                                        x_p, y_p = gcp
                                        lat_1 = y_p[before]
                                        lat_2 = y_p[after]
                                        lon_1 = x_p[before]
#.........这里部分代码省略.........

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import gcpoints [as 别名]
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

map = Basemap(projection='merc', 
              lat_0=0, lon_0=0,
              llcrnrlon=-20.,llcrnrlat=0.,urcrnrlon=180.,urcrnrlat=80.)

map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()

x, y = map.gcpoints(2.3, 48.9, 139.7, 35.6, 300)

print x, y

map.plot(x, y)
plt.show()