Python Basemap.arcgisimage方法代码示例

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import arcgisimage [as 别名]
def ss_plot():
    #Pandas method of importing data frame and getting extents
    conn = psycopg2.connect("dbname='reach_4a' user='root'  host='localhost' port='9000'")
    df = pd.read_sql_query('select easting, northing, texture, sidescan_intensity from mosaic_2014_09', con=conn)
    minE = min(df['easting'])
    maxE = max(df['easting'])
    minN = min(df['northing'])
    maxN = max(df['northing'])
    conn.close()
    print 'Done Importing Data from Database'
    
    #Create grid for countourf plot
    res = 0.25
    grid_x, grid_y = np.meshgrid( np.arange(np.floor(minE), np.ceil(maxE), res), np.arange(np.floor(minN), np.ceil(maxN), res))
    grid_lon, grid_lat = trans(grid_x,grid_y,inverse=True)
    
    #Re-sampling procedure
    m_lon, m_lat = trans(df['easting'].values.flatten(), df['northing'].values.flatten(), inverse=True)
    orig_def = geometry.SwathDefinition(lons=m_lon, lats=m_lat)
    target_def = geometry.SwathDefinition(lons=grid_lon.flatten(), lats=grid_lat.flatten())
    print 'Now Resampling...'
    result = kd_tree.resample_nearest(orig_def, df['sidescan_intensity'].values.flatten(), target_def, radius_of_influence=1, fill_value=None, nprocs = cpu_count())
    print 'Done Resampling!!!' 
    
    #format side scan intensities grid for plotting
    gridded_result = np.reshape(result,np.shape(grid_lon))
    gridded_result = np.squeeze(gridded_result)
    gridded_result[np.isinf(gridded_result)] = np.nan
    gridded_result[gridded_result<=0] = np.nan
    grid2plot = np.ma.masked_invalid(gridded_result)
    
#    x = df['easting']
#    y = df['northing']
#    z = df['sidescan_intensity']
#    
#    xi = df['easting']
#    yi = df['northing']
#    
#    X,Y= np.meshgrid(xi,yi)
#    grid_lon, grid_lat = trans(X,Y,inverse=True)
#    Z = griddata((x, y), z, (X, Y),method='nearest')
#    print 'Done Gridding Data'
    
    
    print 'Now mapping...'
    #Create Figure
    fig = plt.figure(frameon=True)
    ax = plt.subplot(1,1,1)
    map = Basemap(projection='merc', epsg=cs2cs_args.split(':')[1], llcrnrlon=np.min(grid_lon)-0.0009, llcrnrlat=np.min(grid_lat)-0.0009,urcrnrlon=np.max(grid_lon)+0.0009, urcrnrlat=np.max(grid_lat)+0.0009)
    gx,gy = map.projtran(grid_lon,grid_lat)
    map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='World_Imagery', xpixels=1000, ypixels=None, dpi=1200)
    im = map.pcolormesh(gx, gy, grid2plot, cmap='gray',vmin=0.1, vmax=30)
        
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.1)
    cbr = plt.colorbar(im, cax=cax)
    cbr.set_label('Sidescan Intensity [dBw]', size=8)
    for t in cbr.ax.get_yticklabels():
        t.set_fontsize(8)
    plt.show()    

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import arcgisimage [as 别名]
def arc_map(region=[-33.7, 149.4, -31.8, 151.8], res=1000, service='World_Topo_Map',epsg=4326,
            resolution='h',projection='mill'):
    '''
    Create and return the thingy
    Inputs:
        region:  [bot,left,top,right]
        service= one of 'World_Imagery', or 'World_Topo_Map'
            Imagery is more like blue marble, with no site labels
            Topo has inbuilt site labels, which render based on the resolution -> higher res for more zoomed in image.
        epsg: tell arcGIS where to pull image from, affects the image orientation
            Use 4326 for Sydney centric maps
            Use 3577 for Australia centred view
    '''

    # set up basemap projection
    # epsg is used to inform basemap that we will pull an image from ArcGIS servers
    # epsg=28355 does not work...
    # 3577 is centred at middle of Australia...
    # 4326 works for East Australia!
    m = Basemap(projection=projection, epsg=epsg,
                resolution=resolution,
                llcrnrlon=region[1], llcrnrlat=region[0],
                urcrnrlon=region[3], urcrnrlat=region[2])

    # Download backgroud to basemap's projection.
    #m.arcgisimage(service='World_Topo_Map', xpixels=res, verbose=True)
    m.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service=service,
                  xpixels=res, verbose=True)

    return m

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import arcgisimage [as 别名]
def map_rupture(rupture):
    """
    Method for making a simple representation of a Rupture instance.
    This method draws the surface projection of the rupture on a map.

    Args:
        rupture: A Rupture instance.

    """
    rlats = rupture.lats
    rlons = rupture.lons

    minbufx = 0.2
    minbufy = 0.2
    lat1 = np.nanmin(rupture.lats)
    lat2 = np.nanmax(rupture.lats)
    dlat = lat2 - lat1
    bufy = np.max([dlat/2, minbufy])
    lon1 = np.nanmin(rupture.lons)
    lon2 = np.nanmax(rupture.lons)
    dlon = lon2 - lon1
    bufx = np.max([dlon/2, minbufx])
    m = Basemap(llcrnrlat=lat1-bufy,urcrnrlat=lat2+bufy, llcrnrlon=lon1-bufy,urcrnrlon=lon2+bufy)
    m.arcgisimage(service='World_Shaded_Relief', xpixels = 500, verbose= True)
    x, y = m(rlons, rlats)
    m.plot(x, y, 'r')

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import arcgisimage [as 别名]
def EQ_map(lon,lat,M,source, region, ID):
	llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat = -84, -7, -64, 15 
	if  (llcrnrlon < float(lon) < urcrnrlon) and (llcrnrlat < float(lat) < urcrnrlat):
		map = Basemap(llcrnrlon=llcrnrlon,llcrnrlat=llcrnrlat,urcrnrlon=urcrnrlon,urcrnrlat=urcrnrlat, projection='cyl', resolution=None)
		x, y = map(float(lon), float(lat))
		map.scatter(x, y, 500, color="yellow", marker="*", edgecolor="k", zorder=3)
		map.arcgisimage(service='NatGeo_World_Map', xpixels = 1500, verbose= True)
		plt.title('M:'+M+','+region+', Fuente:'+source, fontsize='20')
		name = ID+'.jpg'
		if os.path.exists(name) == False:
			plt.savefig(name)
    			os.chmod(sfile_name,stat.S_IRWXU | stat.S_IRWXG | stat.S_IRWXO)
		else:
			pass
	else:
		map = Basemap(projection='ortho', lat_0=float(lat), lon_0=float(lon), resolution='l')
		map.drawmapboundary(fill_color=(0.478,0.859,0.988))
		map.fillcontinents(color=(0.235,0.705,0.427), lake_color=(0.478,0.859,0.988))
		map.drawparallels(np.arange(-90.,120.,30.))
		map.drawmeridians(np.arange(0.,420.,60.))
		x, y = map(float(lon), float(lat))
		map.scatter(x, y, 800, color="yellow", marker="*", edgecolor="k", zorder=3)
		plt.title('M:'+M+','+region+', Fuente:'+source, fontsize='20')
		name = ID+'.jpg'
		if os.path.exists(name) == False:
			plt.savefig(name)
    			os.chmod(name,stat.S_IRWXU | stat.S_IRWXG | stat.S_IRWXO)
		else:
			pass
	return name

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import arcgisimage [as 别名]
 def make_NYC_basemap(self, ax=None, lllat=40.45, urlat=40.95, lllon=-74.3, urlon=-73.68):
     '''
     This function creates the initial NYC basemap. It is only called once (in the initialization
     of an NYC_school_interactive_map instance), since it has a relatively long run-time (it queries
     the ESRI REST server to get a relatively high-resolution basemap image).
     ''' 
         
     fig = plt.figure(figsize=(12,12))
     ax = plt.subplot(111)
         
     m = Basemap(ax = ax, 
                    lon_0=(urlon + lllon)/2,
                    lat_0=(urlat + lllat)/2,
                    llcrnrlat=lllat, urcrnrlat=urlat,
                    llcrnrlon=lllon, urcrnrlon=urlon,
                    resolution='f', epsg=3857) ###epsg is the projection code for ESRI world shaded relief basemap
 
     #m.arcgisimage queries the ESRI REST API (server below) to get the basemap 
     #http://server.arcgisonline.com/arcgis/rest/services
     m.arcgisimage(service='World_Shaded_Relief', xpixels = 1500)
 
     # Add county lines, since boroughs map onto counties
     m.drawcounties(linewidth=0.5)
     
     return m, ax, fig

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import arcgisimage [as 别名]
def ss_plot():
    #Pandas method of importing data frame and getting extents
    db_connect="dbname='reach_4a' user='root'  host='localhost' port='9000'"
    conn = psycopg2.connect(db_connect)
    df = pd.read_sql_query('SELECT * from mb_may_2012_1m tt inner join (	SELECT s.easting, s.northing, s.texture, s.sidescan_intensity  FROM ss_2012_05 s) ss on tt.easting=ss.easting and tt.northing=ss.northing;', con=conn)
    minE = df['easting'].min()[0]
    maxE = df['easting'].max()[0]
    minN = df['northing'].min()[0]
    maxN = df['northing'].max()[0]
    conn.close()
    print 'Done Importing Data from Database'
    
    #Create grid for countourf plot
    res = 1
    grid_x, grid_y = np.meshgrid( np.arange(np.floor(minE), np.ceil(maxE), res), np.arange(np.floor(minN), np.ceil(maxN), res))
    grid_lon, grid_lat = trans(grid_x,grid_y,inverse=True)
    
    #Re-sampling procedure
    m_lon, m_lat = trans(df['easting'].values.flatten(), df['northing'].values.flatten(), inverse=True)
    orig_def = geometry.SwathDefinition(lons=m_lon, lats=m_lat)
    target_def = geometry.SwathDefinition(lons=grid_lon.flatten(), lats=grid_lat.flatten())
    print 'Now Resampling...'
    result = kd_tree.resample_nearest(orig_def, df['sidescan_intensity'].values.flatten(), target_def, radius_of_influence=1, fill_value=None, nprocs = cpu_count())
    print 'Done Resampling!!!' 
    
    #format side scan intensities grid for plotting
    gridded_result = np.reshape(result,np.shape(grid_lon))
    gridded_result = np.squeeze(gridded_result)
    gridded_result[np.isinf(gridded_result)] = np.nan
    gridded_result[gridded_result<=0] = np.nan
    grid2plot = np.ma.masked_invalid(gridded_result)
       
    
    print 'Now mapping...'
    #Create Figure
    fig = plt.figure(frameon=True)
    ax = plt.subplot(1,1,1)
    map = Basemap(projection='merc', epsg=cs2cs_args.split(':')[1], llcrnrlon=np.min(grid_lon)-0.0009, llcrnrlat=np.min(grid_lat)-0.0009,urcrnrlon=np.max(grid_lon)+0.0009, urcrnrlat=np.max(grid_lat)+0.0009)
    gx,gy = map.projtran(grid_lon,grid_lat)
    map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='World_Imagery', xpixels=1000, ypixels=None, dpi=1200)
    im = map.pcolormesh(gx, gy, grid2plot, cmap='gray',vmin=0.1, vmax=30)
        
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.1)
    cbr = plt.colorbar(im, cax=cax)
    cbr.set_label('Sidescan Intensity [dBw]', size=8)
    for t in cbr.ax.get_yticklabels():
        t.set_fontsize(8)
    plt.savefig(r'C:\workspace\Texture_Classification\output\May2012_1m_sidescan_intensity.png')    

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import arcgisimage [as 别名]
    def spatialSubPlot(self, platformLineStringHash, ax, e, resolution='l'):
        '''
        Make subplot of tracks for all the platforms within the time range. 
        '''
        m = Basemap(llcrnrlon=e[0], llcrnrlat=e[1], urcrnrlon=e[2], urcrnrlat=e[3], projection='cyl', resolution=resolution, ax=ax)
        ##m.wmsimage('http://www.gebco.net/data_and_products/gebco_web_services/web_map_service/mapserv?', layers=['GEBCO_08_Grid'])    # Works, but coarse
        m.arcgisimage(server='http://services.arcgisonline.com/ArcGIS', service='Ocean_Basemap')
 
        for pl, LS in platformLineStringHash.iteritems():
            x,y = zip(*LS)
            m.plot(x, y, '-', c=self._getColor(pl), linewidth=3)

        if self.args.mapLabels:
            m.drawparallels(np.linspace(e[1],e[3],num=3), labels=[True,False,False,False], linewidth=0)
            m.drawmeridians(np.linspace(e[0],e[2],num=3), labels=[False,False,False,True], linewidth=0)

        return ax

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import arcgisimage [as 别名]
def print_contour_map(cs2cs_args, humlon, humlat, glon, glat, datm, sonpath, p, vmin, vmax): #merge, 

    #levels = [0,0.25,0.5,0.75,1.25,1.5,1.75,2,3,5]
          
    print("drawing and printing map ...")
    fig = plt.figure(frameon=False)
    map = Basemap(projection='merc', epsg=cs2cs_args.split(':')[1], #26949,
     resolution = 'i', #h #f
     llcrnrlon=np.min(humlon)-0.001, llcrnrlat=np.min(humlat)-0.001,
     urcrnrlon=np.max(humlon)+0.001, urcrnrlat=np.max(humlat)+0.001)

    #if dogrid==1:
    gx,gy = map.projtran(glon, glat)

    ax = plt.Axes(fig, [0., 0., 1., 1.], )
    ax.set_axis_off()
    fig.add_axes(ax)

    try:
       map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='World_Imagery', xpixels=1000, ypixels=None, dpi=300)
    except:
       map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='ESRI_Imagery_World_2D', xpixels=1000, ypixels=None, dpi=300)
    #finally:
    #   print "error: map could not be created..."
             
    #if dogrid==1:
    if 2>1:
       if datm.size > 25000000:
          print("matrix size > 25,000,000 - decimating by factor of 5 for display")
          #map.contourf(gx[::5,::5], gy[::5,::5], datm[::5,::5], levels, cmap='YlOrRd')
          map.pcolormesh(gx, gy, datm, cmap='YlOrRd', vmin=vmin, vmax=vmax)
       else:
          #map.contourf(gx, gy, datm, levels, cmap='YlOrRd')
          map.pcolormesh(gx[::5,::5], gy[::5,::5], datm[::5,::5], cmap='pink', vmin=vmin, vmax=vmax)
                   
    custom_save2(sonpath,'class_map_imagery'+str(p))
    del fig 

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

#.........这里部分代码省略.........
                    cb.set_label(units,fontsize=8)#,labelpad=5)
                    cb.ax.xaxis.set_ticks_position('top')
                    for t in cb.ax.yaxis.get_ticklabels():
                        t.set_fontsize(8)
                else:
                    if plot_step:
                        xaxis = ax0_colorbar.xaxis
                        ax0_colorbar.xaxis.set_major_locator(plt.NullLocator())
                        ax0_colorbar.yaxis.set_ticks_position('right')
                        for t in ax0_colorbar.yaxis.get_ticklabels():
                            t.set_fontsize(8)
                    else:
                        cbFormatter = FormatStrFormatter('%.2f')
                        cb = plt.colorbar(cs0, cax=ax0_colorbar, ticks=[min(rangez), max(rangez)], format=cbFormatter, orientation='vertical')
                        cb.set_label(units,fontsize=8)#,labelpad=5)
                        cb.ax.xaxis.set_ticks_position('top')
                        for t in cb.ax.yaxis.get_ticklabels():
                            t.set_fontsize(8)


                # Rotate and show the date with date formatter in the last plot of all the groups
                logger.debug('rotate and show date with date formatter')
                if name is parm[-1] and group is self.plotGroupValid[-1]:
                    x_fmt = self.DateFormatter(self.scale_factor)
                    if plot_scatter:
                        ax = ax1_plot
                        # Don't show on the upper contour plot
                        ax0_plot.xaxis.set_ticks([])
                    else:
                        ax = ax0_plot
                    ax.xaxis.set_major_formatter(x_fmt)
                    # Rotate date labels
                    for label in ax.xaxis.get_ticklabels():
                        label.set_rotation(10)
                else:
                    ax0_plot.xaxis.set_ticks([])
                    if plot_scatter:  
                        ax1_plot.xaxis.set_ticks([])

                if plot_scatter:
                    j+=4
                else: 
                    j+=2


        # plot tracks
        ax = plt.Subplot(fig, map_gs[:])
        fig.add_subplot(ax, aspect='equal')

        logger.debug('computing activity extents')
        z = []
        maptracks = None

        z, points, maptracks = self.getMeasuredPPData(startDatetime, endDatetime, pn, self.plotDotParmName)

        # get the percentile ranges for this to autoscale
        pointsnp = np.array(points)
        lon = pointsnp[:,0]
        lat = pointsnp[:,1]

        if self.extent:
            ltmin = self.extent[1]
            ltmax = self.extent[3]
            lnmin = self.extent[0]
            lnmax = self.extent[2]
            lndiff = abs(lnmax - lnmin)
            ltdiff = abs(ltmax - ltmin)
            logger.debug("lon diff %f lat diff %f" %(lndiff, ltdiff))
            mindeg = .02
            paddeg = .01
            if lndiff < mindeg :
                lnmin -= mindeg
                lnmax += mindeg
            if ltdiff < mindeg:
                ltmin -= mindeg
                ltmax += mindeg
            e = (lnmin - paddeg, ltmin - paddeg, lnmax + paddeg, ltmax + paddeg)
        else:
            # default map to Monterey Bay region
            ltmin = 36.61
            ltmax = 36.97
            lnmin = -122.21
            lnmax = -121.73
            e = (lnmin, ltmin, lnmax, ltmax)

        logger.debug('Extent found %f,%f,%f,%f)' % (e[0], e[1],e[2],e[3]))
        # retry up to 5 times to get the basemap
        for i in range(0, 5):
            mp = Basemap(llcrnrlon=e[0], llcrnrlat=e[1], urcrnrlon=e[2], urcrnrlat=e[3], projection='cyl', resolution='l', ax=ax)
            try:
                ##mp.wmsimage('http://www.gebco.net/data_and_products/gebco_web_services/web_map_service/mapserv?', layers=['GEBCO_08_Grid'])                            # Works, but coarse
                mp.arcgisimage(server='http://services.arcgisonline.com/ArcGIS', service='Ocean_Basemap')
                mp.drawparallels(np.linspace(e[1],e[3],num=3), labels=[True,False,False,False], fontsize=8, linewidth=0)
                mp.drawmeridians(np.linspace(e[0],e[2],num=3), labels=[False,False,False,True], fontsize=8, linewidth=0)
            except Exception, e:
                logger.error('Could not download ocean basemap ')
                mp = None

            if mp is not None :
                break

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

#.........这里部分代码省略.........
        self.msg_min_long = StringVar()
        zone_min_long = Entry(l_pos,background='white',textvariable=self.msg_min_long)
        zone_min_long.grid(row = 2, column = 3)

        #label minutes
        label_min_lat=Label(l_pos, text="'")
        label_min_lat.grid(row=1,column = 4)
        label_min_long=Label(l_pos, text="'")
        label_min_long.grid(row=2,column = 4)

        #zones secondes
        self.msg_sec_lat = StringVar()
        zone_sec_lat = Entry(l_pos,background='white',textvariable=self.msg_sec_lat)
        zone_sec_lat.grid(row = 1, column = 5)
        self.msg_sec_long = StringVar()
        zone_sec_long = Entry(l_pos,background='white',textvariable=self.msg_sec_long)
        zone_sec_long.grid(row = 2, column = 5)

        #label secondes
        label_sec_lat=Label(l_pos, text='\"')
        label_sec_lat.grid(row=1,column = 6)
        label_sec_long=Label(l_pos, text='\"')
        label_sec_long.grid(row=2,column = 6)

        #bouton envoie coordonnées
        bouton_envoyer = Button(l_pos, text= "Envoyer les coordonnées", command = self.envoyer_pos)
        bouton_envoyer.grid(row=1, column=9)

        #Affichage de la carte
        fig = matplotlib.pyplot.figure(figsize=(12, 6))

        self.map = Basemap(llcrnrlon=-4.476,llcrnrlat=48.417,urcrnrlon=-4.469,urcrnrlat=48.42, epsg=2154)

        self.map.arcgisimage(service='ESRI_Imagery_World_2D', xpixels = 1500, verbose= True)

        self.map.drawcountries()
        self.map.drawcoastlines(linewidth=.5)

        self.canvas = FigureCanvasTkAgg(fig, master=l_map)
        self.canvas.get_tk_widget().pack(side=tk.TOP, fill=X, expand=1)
        self.canvas.mpl_connect('button_press_event', self.onclick_map)


#gestion de la batterie
        #Frame de la batterie
        l_bat = Label(onglet4)
        l_bat.pack(fill="both", expand = "yes")
        #Voltage
        label_volt = Label(l_bat, text="Voltage battery")
        label_volt.grid(row= 0, column =0)
        zone_reception_volt = Text(l_bat,width=25, height= 1,background='white')
        zone_reception_volt.grid(row= 0, column =1) 
        #Etat
        label_cur = Label(l_bat, text="Current battery")
        label_cur.grid(row= 1, column =0)
        zone_reception_cur = Text(l_bat,width=25, height= 1,background='white')
        zone_reception_cur.grid(row= 1, column =1)    
        #Etat
        label_rem = Label(l_bat, text="battery remaining")
        label_rem.grid(row= 2, column =0)
        zone_reception_rem = Text(l_bat,width=25, height= 1,background='white')
        zone_reception_rem.grid(row= 2, column =1) 
    
#gestion de la navigation
        #Frame de la vitesse (tangage, roulis, lacet)    
        l_vit = LabelFrame(onglet2, text="Vitesse")

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

#.........这里部分代码省略.........
        bouton_envoyer = Button(l_pos, text= "Envoyer les coordonnées", command = self.envoyer_pos)
        bouton_envoyer.grid(row=1, column=9)
        '''

        self.var = IntVar()
        R1 = Radiobutton(l_mis, text="waypoint", variable=self.var, value=1)
        R1.grid(row=0, column=3, sticky=W)

        R2 = Radiobutton(l_mis, text="spirale carrée", variable=self.var, value=2)
        R2.grid(row=1, column=3, sticky=W)

        R3 = Radiobutton(l_mis, text="lacet", variable=self.var, value=3)
        R3.grid(row=2, column=3, sticky=W)

        launch_button=Button(l_mis, text="Launch mission", command=self.launch, height = 2, width = 20, padx=20)
        launch_button.grid(row=0, column=4, sticky=W)

        clear_button=Button(l_mis, text="Clear mission", command=self.clear_waypoint, height = 2, width = 20, padx=20)
        clear_button.grid(row=1, column=4)

        takeoff_button=Button(l_mis, text="Take off!", command=self.takeoff, height = 2, width = 20, padx=20)
        takeoff_button.grid(row=0, column=5)

        land_button=Button(l_mis, text="Land", command=self.land, height = 2, width = 20, padx=20)
        land_button.grid(row=1, column=5)

        #endregion
        print("Initiailsé")
        #region AFFICHER_CARTE
        fig = plt.figure(figsize=(12, 6))

        self.map = Basemap(llcrnrlon=-4.476,llcrnrlat=48.417, urcrnrlon=-4.469, urcrnrlat=48.42, epsg=2154)
        try:
            self.map.arcgisimage(service='ESRI_Imagery_World_2D', xpixels=2000, verbose=True)
            print("Carte chargé online")
        except:
            im = plt.imread('ensta.png')
            self.map.imshow(im)
            print("Carte chargé offline")
        self.mapplot = None

        #self.map.drawcountries()
        #self.map.drawcoastlines(linewidth=.5)

        self.canvas = FigureCanvasTkAgg(fig, master=l_map)
        self.canvas.get_tk_widget().pack(side=tk.TOP, fill=X, expand=1)
        self.canvas.mpl_connect('button_press_event', self.onclick_map)

        

        #bouton_waypoint = Button(self.canvas, text= "Envoyer les waypoint", command = self.envoyer_way)
        #endregion

        #region BATTERIE
        '''#Frame de la batterie
        l_bat = Label(onglet4)
        l_bat.pack(fill="both", expand = "yes")
        #Voltage
        label_volt = Label(l_bat, text="Voltage battery")
        label_volt.grid(row= 0, column =0)
        zone_reception_volt = Text(l_bat,width=25, height= 1,background='white')
        zone_reception_volt.grid(row= 0, column =1) 
        #Etat
        label_cur = Label(l_bat, text="Current battery")
        label_cur.grid(row= 1, column =0)
        zone_reception_cur = Text(l_bat,width=25, height= 1,background='white')

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

#.........这里部分代码省略.........

		self.SetZoomMap(-105,-90.5,29,38)

	def WDquakes(self):

		self.SetZoomMap(-179.9,179.9,-89.9,89.9)

	def PrettyMap(self,lat1,lat2,lon1,lon2):

		'''Display topography and street names on the current map - could be modified for various map types'''

		#scale the longitude and latiude grid increments
		latinc = (max(lat1,lat2)-min(lat1,lat2))/5
		loninc = abs((lon2-lon1)/4)
		lon0 = int(abs(lon2)-abs(lon1))/2
		lat0 = int((lat2-lat1)/2)

		if 2.0 < abs(lat2-lat1) < 6.0:
			res = 'i'
		elif 0.1 < abs(lat2-lat1) < 1.0:
			res = 'h'
		elif 0 < abs(lat2-lat1) < 0.1:
			res = 'f'
		else:
			res = 'l'

		self.a.clear()

		#Determine the new sizes of the map objects, so that they don't crowd the map
		mt_width,mt_rad,min_dot,max_dot,min_quake,max_quake = quaketools.determinemtsize(self.minlon,self.maxlon,self.minlat,self.maxlat)

		#Map setup and quakes
		self.map = Basemap(ax=self.a,lat_0=lat0,lon_0=lon0,resolution ='l',llcrnrlon=lon1,llcrnrlat=lat1,urcrnrlon=lon2,urcrnrlat=lat2)
		self.map.arcgisimage(service='NatGeo_World_Map',verbose=False,xpixels=600)

		#plot some extra information

		self.map.drawparallels(np.arange(lat1,lat2,latinc),labels=[1,1,0,0],linewidth=0.5,fontsize=4)
		self.map.drawmeridians(np.arange(lon1,lon2,loninc),labels=[0,0,0,1],linewidth=0.5,fontsize=4)
		self.map.drawcountries()
		self.map.drawcoastlines(linewidth=0.1)

		#replot the quakes

		if self.momenttensors == True:

			#plot the moment tensors and redraw
			#self.quakesplotted = quaketools.plot_events(self.map,self.a,self.quakes)
			self.mtlines, self.MTs, self.quakedots, xs, ys, urls = quaketools.plot_mt(self.map,self.a,self.f,self.quakes,self.mts,self.events,mt_width=mt_width,radius=mt_rad,angle_step=40,llat=lat1,ulat=lat2,llon=lon1,ulon=lon2)
			Browse.updatedata(xs,ys,urls,mt_rad)

		else:
			#only plotting events, so continue
			self.quakesplotted = quaketools.plot_events(self.map,self.a,self.quakes,llat=self.minlat,ulat=self.maxlat,llon=self.minlon,ulon=self.maxlon,dist_bt=100,min_size=min_quake,max_size=max_quake)

		self.canvas.draw()


	def setMTs(self):

		self.momenttensors = True
		self.removemapobjs()
		self.setquakesmanual(t1=self.starttime)

	def setevents(self):

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

#.........这里部分代码省略.........
          if dogrid==1:
             map.pcolormesh(gx, gy, datm, cmap='YlOrRd', vmin=0.5, vmax=2)
             del dat
          else: 
             ## draw point cloud
             x,y = map.projtran(humlon, humlat)
             map.scatter(x.flatten(), y.flatten(), 0.5, merge.flatten(), cmap='YlOrRd', linewidth = '0')

          custom_save(sonpath,'class_map'+str(p))
          del fig 

       except:
          print "error: map could not be created..."

       kml = simplekml.Kml()
       ground = kml.newgroundoverlay(name='GroundOverlay')
       ground.icon.href = 'class_map'+str(p)+'.png'
       ground.latlonbox.north = np.min(humlat)-0.001
       ground.latlonbox.south = np.max(humlat)+0.001
       ground.latlonbox.east =  np.max(humlon)+0.001
       ground.latlonbox.west =  np.min(humlon)-0.001
       ground.latlonbox.rotation = 0

       kml.save(sonpath+'class_GroundOverlay'+str(p)+'.kml')

    if dowrite==1:

       X = []; Y = []; S = [];
       for p in xrange(len(class_fp)):
          dat = np.genfromtxt(sonpath+'x_y_class'+str(p)+'.asc', delimiter=' ')
          X.append(dat[:,0])
          Y.append(dat[:,1])
          S.append(dat[:,2])
          del dat

       # merge flatten and stack
       X = np.asarray(np.hstack(X),'float')
       X = X.flatten()

       # merge flatten and stack
       Y = np.asarray(np.hstack(Y),'float')
       Y = Y.flatten()

       # merge flatten and stack
       S = np.asarray(np.hstack(S),'float')
       S = S.flatten()

       humlon, humlat = trans(X, Y, inverse=True)

       if dogrid==1:
          grid_x, grid_y = np.meshgrid( np.arange(np.min(X), np.max(X), res), np.arange(np.min(Y), np.max(Y), res) )  

          dat = griddata(np.c_[X.flatten(),Y.flatten()], S.flatten(), (grid_x, grid_y), method='nearest')

          ## create mask for where the data is not
          tree = KDTree(np.c_[X.flatten(),Y.flatten()])
          dist, _ = tree.query(np.c_[grid_x.ravel(), grid_y.ravel()], k=1)
          dist = dist.reshape(grid_x.shape)

       del X, Y

       if dogrid==1:
          ## mask
          dat[dist> 1 ] = np.nan

          del dist, tree

          dat[dat==0] = np.nan
          dat[np.isinf(dat)] = np.nan

          datm = np.ma.masked_invalid(dat)

          glon, glat = trans(grid_x, grid_y, inverse=True)
          del grid_x, grid_y

       levels = [0.5,0.75,1.25,1.5,1.75,2,3]

       try:
          print "drawing and printing map ..."
          fig = plt.figure()
          map = Basemap(projection='merc', epsg=cs2cs_args.split(':')[1],
           resolution = 'i',
           llcrnrlon=np.min(humlon)-0.001, llcrnrlat=np.min(humlat)-0.001,
           urcrnrlon=np.max(humlon)+0.001, urcrnrlat=np.max(humlat)+0.001)

          map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='World_Imagery', xpixels=1000, ypixels=None, dpi=300)
          if dogrid==1:
             gx,gy = map.projtran(glon, glat)

          if dogrid==1:
             map.contourf(gx, gy, datm, levels, cmap='YlOrRd')
          else: 
             ## draw point cloud
             x,y = map.projtran(humlon, humlat)
             map.scatter(x.flatten(), y.flatten(), 0.5, S.flatten(), cmap='YlOrRd', linewidth = '0')

          custom_save2(sonpath,'class_map_imagery'+str(p))
          del fig 
       except:
          print "error: map could not be created..."

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

#.........这里部分代码省略.........

             except:
                print "plot could not be produced"

          if doplot==1:
             try:

                fig = plt.figure()
                s=plt.scatter(Zes[labels==0],Zns[labels==0],marker='o',c='k', s=10, linewidth=0, vmin=0, vmax=8);
                s=plt.scatter(Zes[labels==1],Zns[labels==1],marker='o',c='r', s=10, linewidth=0, vmin=0, vmax=8);
                s=plt.scatter(Zes[labels==2],Zns[labels==2],marker='o',c='b', s=10, linewidth=0, vmin=0, vmax=8);
                custom_save(sonpath,'rgh_hard_kmeans'+str(p))
                del fig

             except:
                print "plot could not be produced"

          if doplot==1:
             try:

                print "drawing and printing map ..."
                fig = plt.figure(frameon=False)
                #fig.subplots_adjust(wspace = 0.4, hspace=0.4)
                map = Basemap(projection='merc', epsg=cs2cs_args.split(':')[1], #epsg=26949,
                   resolution = 'i', #h #f
                   llcrnrlon=np.min(Zlon)-0.0001, llcrnrlat=np.min(Zlat)-0.0001,
                   urcrnrlon=np.max(Zlon)+0.0001, urcrnrlat=np.max(Zlat)+0.0001)

                # draw point cloud
                x,y = map.projtran(Zlon, Zlat)

                cs = map.scatter(x.flatten(), y.flatten(), 1, rough.flatten(), linewidth=0, vmin=0, vmax=8)

                map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='World_Imagery', xpixels=1000, ypixels=None, dpi=300)
  
                cbar = map.colorbar(cs,location='bottom',pad="5%")
                cbar.set_label('E1')
                cbar.set_ticks([0,2,4,6,8])

                custom_save(sonpath,'map_rgh'+str(p))
                del fig 

             except:
                print "plot could not be produced"

          if doplot==1:
             try:
                fig = plt.figure()
                #fig.subplots_adjust(wspace = 0.4, hspace=0.4)
                map = Basemap(projection='merc', epsg=cs2cs_args.split(':')[1],
                   resolution = 'i', #h #f
                   llcrnrlon=np.min(Zlon)-0.0001, llcrnrlat=np.min(Zlat)-0.0001,
                   urcrnrlon=np.max(Zlon)+0.0001, urcrnrlat=np.max(Zlat)+0.0001)

                # draw point cloud
                x,y = map.projtran(Zlon, Zlat)

                cs = map.scatter(x.flatten(), y.flatten(), 1, hard.flatten(), linewidth=0, vmin=0, vmax=8)

                map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='World_Imagery', xpixels=1000, ypixels=None, dpi=300)

                cbar = map.colorbar(cs,location='bottom',pad="5%")
                cbar.set_label('E2')
                cbar.set_ticks([0,2,4,6,8])

                custom_save(sonpath,'map_hard'+str(p))

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

# test arcgisimage method for retrieving images from web map servers.

plt.figure(figsize=(8,8))
epsg=4326
m=Basemap(projection='cyl',llcrnrlon=-90,llcrnrlat=30,urcrnrlon=-60,urcrnrlat=60,resolution='i')
# default 'blue marble' image.
m.arcgisimage(verbose=True)
m.drawmeridians(np.arange(-180,180,10),labels=[0,0,0,1],color='y')
m.drawparallels(np.arange(-90,90,10),labels=[1,0,0,0],color='y')
m.drawcoastlines(linewidth=0.5,color='0.5')
plt.title('test WMS map background EPSG=%s'%epsg)

plt.figure(figsize=(8,8))
epsg = 3413; width = 18000.e3; height = 18000.e3
m=Basemap(epsg=epsg,resolution='l',width=width,height=height)
# default 'blue marble' image.
m.arcgisimage()
m.drawmeridians(np.arange(-180,180,30),labels=[0,0,0,1],color='y')
m.drawparallels(np.arange(0,80,10),labels=[1,0,0,0],color='y')
m.drawcoastlines(linewidth=0.5,color='0.5')
plt.title('test WMS map background EPSG=%s' % epsg)

plt.figure(figsize=(10,6.6666))
epsg = 2263; width=600.e3; height = 400.e3
m=Basemap(epsg=epsg,resolution='h',width=width,height=height)
# specify a different server.
m.arcgisimage(server='http://maps.ngdc.noaa.gov',service='etopo1',verbose=True)