Python補間による等値線図の描画

3804 ワード

最近、限られたサイトに基づいて補間等値線図を描く必要があります.私はコードに注釈をつけて、すでに他の人のクエリーを用意して、文の中で言及したデータについてtxtフォーマットで、私も直接データを下に貼りました.まとめ:地図に等値線を描画するには:
  • は、基本的な図面フレームワークを決定する.
  • 収集データを取得し、地図とマッピングし、マッピングデータに基づいて補間する.scipy.interpolate.griddataパケットの補間は比較的速く、よく使われる3つの補間方法はliner(三角形ベースの線形補間法)、cubic(三角形ベースの3次補間法)、nearest(最近の隣接補間法)であり、これらの方法はデータ点に一致するサーフェスタイプを定義し、'cubic'方法は平滑なサーフェスを生成し、'linear'と'nearest'はそれぞれ1次導関数と0次導関数が不連続である.
  • グリッド補間データに基づく図面
  • # -*- coding: utf-8 -*-
    @author: Adwiy Wang
    import numpy as np
    import pandas as pd
    from matplotlib.mlab import griddata
    from mpl_toolkits.basemap import Basemap
    import matplotlib.pyplot as plt
    from matplotlib.colors import Normalize
    from scipy.interpolate import griddata as gd
    
    #         
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, axisbg='w', frame_on=False)
    
    #     
    data = pd.read_csv('datam.txt', delim_whitespace=True)
    norm = Normalize()
    
    #       
    lllon = 21
    lllat = -18
    urlon = 34
    urlat = -8
    
    #     
    m = Basemap(
        projection = 'merc',
        llcrnrlon = lllon, llcrnrlat = lllat, urcrnrlon = urlon, urcrnrlat = urlat,
        resolution='h')
    
    #              
    data['projected_lon'], data['projected_lat'] = m(*(data.Lon.values, data.Lat.values))
    
    #           
    numcols, numrows = 1000, 1000
    xi = np.linspace(data['projected_lon'].min(), data['projected_lon'].max(), numcols)
    yi = np.linspace(data['projected_lat'].min(), data['projected_lat'].max(), numrows)
    xi, yi = np.meshgrid(xi, yi)
    
    #   
    x, y, z = data['projected_lon'].values, data['projected_lat'].values, data.Z.values
    zi = gd(
        (data[['projected_lon', 'projected_lat']]),
        data.Z.values,
        (xi, yi),
        method='cubic')
    
    #       
    m.drawmapboundary(fill_color = 'white')
    m.fillcontinents(color='#C0C0C0', lake_color='#7093DB')
    m.drawcountries(
        linewidth=.75, linestyle='solid', color='#000073',
        antialiased=True,
        ax=ax, zorder=3)
    
    m.drawparallels(
        np.arange(lllat, urlat, 2.),
        color = 'black', linewidth = 0.5,
        labels=[True, False, False, False])
    m.drawmeridians(
        np.arange(lllon, urlon, 2.),
        color = '0.25', linewidth = 0.5,
        labels=[False, False, False, True])
    
    #       
    con = m.contourf(xi, yi, zi, zorder=4, alpha=0.6, cmap='jet')
    #      
    m.scatter(
        data['projected_lon'],
        data['projected_lat'],
        color='#545454',
        edgecolor='#ffffff',
        alpha=.75,
        s=50 * norm(data['Z']),
        cmap='jet',
        ax=ax,
        vmin=zi.min(), vmax=zi.max(), zorder=4)
    
    #     、     
    cbar = plt.colorbar(con,orientation='horizontal', fraction=.057, pad=0.05)
    cbar.set_label("Mean Rainfall - mm")
    
    m.drawmapscale(
        24., -9., 28., -13,
        100,
        units='km', fontsize=10,
        yoffset=None,
        barstyle='fancy', labelstyle='simple',
        fillcolor1='w', fillcolor2='#000000',
        fontcolor='#000000',
        zorder=5)
    
    plt.title("Mean Rainfall")
    plt.savefig("rainfall.png", format="png", dpi=300, transparent=True)
    plt.show()
    

    データファイル:datam.txt
    Lon   Lat     Z
    32.6  -13.6   41
    27.1  -16.9   43
    32.7  -10.2   46
    24.2  -13.6   33
    28.5  -14.4   43
    28.1  -12.6   33
    27.9  -15.8   46
    24.8  -14.8   44
    31.1  -10.2   35
    25.9  -13.5   24
    29.1   -9.8   10
    25.8  -17.8   39
    33.2  -12.3   44
    28.3  -15.4   46
    27.6  -16.1   47
    28.9  -11.1   31
    31.3   -8.9   39
    31.9  -13.3   45
    23.1  -15.3   31
    31.4  -11.9   39
    27.1  -15.0   42
    24.4  -11.8   15
    28.6  -13.0   39
    31.3  -14.3   44
    23.3  -16.1   39
    30.2  -13.2   38
    24.3  -17.5   32
    26.4  -12.2   23
    23.1  -13.5   27