テスト用の投稿(2)

テスト用の投稿の２回目。今度はソースコードを記載してみる。


```python:ex5.py
# -*- coding: utf-8 -*-
# Authors : 
# ==============================================================================
import numpy as np
import pandas as pd
import pymap3d as pm
from netCDF4 import Dataset
#import plotly.graph_objs as go
#from plotly.offline import plot
import matplotlib.pyplot as plt
import matplotlib.ticker as ptick
from matplotlib.collections import LineCollection


def Etopo(lon_area, lat_area, resolution):
    ### Input
    # resolution: resolution of topography for both of longitude and latitude [deg]
    # (Original resolution is 0.0167 deg)
    # lon_area and lat_area: the region of the map which you want like [100, 130], [20, 25]
    ###

    ### Output
    # Mesh type longitude, latitude, and topography data
    ###

    # Read NetCDF data
    # https://www.ngdc.noaa.gov/mgg/global/
    data = Dataset("./data/ETOPO1_Ice_g_gdal.grd", "r")

    # Get data
    lon_range = data.variables['x_range'][:]
    lat_range = data.variables['y_range'][:]
    alt_range = data.variables['z_range'][:]
    spacing = data.variables['spacing'][:]
    dimension = data.variables['dimension'][:]
    z = data.variables['z'][:]

    lon_num = dimension[0]
    lat_num = dimension[1]

    # Prepare array
    lon_input = np.zeros(lon_num); lat_input = np.zeros(lat_num)
    for i in range(lon_num):
        lon_input[i] = lon_range[0] + i * spacing[0]

    for i in range(lat_num):
        lat_input[i] = lat_range[0] + i * spacing[1]

    # Create 2D array
    lon, lat = np.meshgrid(lon_input, lat_input)

    # Convert 2D array from 1D array for z value
    alt = np.reshape(z, (lat_num, lon_num))

    # Skip the data for resolution
    if ((resolution < spacing[0]) | (resolution < spacing[1])):
        print('Set the highest resolution')
    else:
        skip = int(resolution/spacing[0])
        lon = lon[::skip,::skip]
        lat = lat[::skip,::skip]
        alt = alt[::skip,::skip]

    alt = alt[::-1]

    # Select the range of map
    range1 = np.where((lon>=lon_area[0]) & (lon<=lon_area[1]))
    lon = lon[range1]; lat = lat[range1]; alt = alt[range1]
    range2 = np.where((lat>=lat_area[0]) & (lat<=lat_area[1]))
    lon = lon[range2]; lat = lat[range2]; alt = alt[range2]
    alt = np.clip(alt,-5,2000)
    #alt[alt < 0] = -1000
    

    # Convert 2D again
    lon_num = len(np.unique(lon))
    lat_num = len(np.unique(lat))
    lon = np.reshape(lon, (lat_num, lon_num))
    lat = np.reshape(lat, (lat_num, lon_num))
    alt = np.reshape(alt, (lat_num, lon_num))
    
    return lon, lat, alt


def lines_SS_POI(df_SS, df_POI):
    end_pos = df_SS.loc[:,['lon','lat']].values.tolist()
    start_pos = df_POI.loc[:,['lon','lat']].values.tolist()*len(end_pos)
    
    lines = [ [sp,ep] for sp,ep in zip(start_pos, end_pos)]
    lc = LineCollection(lines, colors='red')
    return lc


def plot_map(ax, ss_idx, lon_topo, lat_topo, alt, df_SS, df_POI):
    alt_tmp = alt.copy()
    alt_tmp[alt < 0] = -1000
    
    cont = ax.contour(lon_topo, lat_topo, alt_tmp, levels=10, linewidths=.5, colors=['black'])
    ax.pcolormesh(lon_topo, lat_topo, alt_tmp, cmap='YlGn')
    ax.plot(df_SS['lon'], df_SS['lat'], 'D', color='red', ms=10)
    ax.plot(df_POI['lon'], df_POI['lat'], 'o', color='blue', ms=10)
    end_pos = df_SS.loc[:,['lon','lat']].values.tolist()
    start_pos = df_POI.loc[:,['lon','lat']].values.tolist()*len(end_pos)

    yy = (start_pos[ss_idx][1],end_pos[ss_idx][1])
    xx = (start_pos[ss_idx][0],end_pos[ss_idx][0])
    ax.plot(xx,yy,'r')
    

def pol_cut(df_SS, df_POI, az_width_oneside, rg_max_km):
    [ss_lon, ss_lat, ss_alt] = df_SS.loc[ss_idx,['lon','lat','alt']]
    [poi_lon, poi_lat, poi_alt] = df_POI.loc[0,['lon','lat','alt']]
    az_t, el_t, rg_t = pm.geodetic2aer(poi_lat, poi_lon, poi_alt, ss_lat, ss_lon, ss_alt)
    az_beg = az_t - az_width_oneside
    az_end = az_t + az_width_oneside
    
    az,el,rg = pm.geodetic2aer(lat_topo, lon_topo, alt, ss_lat, ss_lon, ss_alt)
    az_1d = az.reshape(-1)
    el_1d = el.reshape(-1)
    rg_1d = rg.reshape(-1)
    alt_1d = alt.reshape(-1)        
        
    df_pol = pd.DataFrame(data={"az":az_1d,"el":el_1d, "rg":rg_1d, "alt":alt_1d})

    if (az_beg > 0) & (az_end < 360):
        df_pol_cut = df_pol[((az_beg < df_pol['az'] ) & (df_pol['az'] < az_end)) &
                            (df_pol['rg'] < rg_max_km*1000)]
    elif az_end > 360 :
        df_pol_cut = df_pol[(((az_beg < df_pol['az'] ) & (df_pol['az'] < 360)) |
                             ((0 < df_pol['az'] ) & (df_pol['az'] < az_end-360))) &
                            (df_pol['rg'] < rg_max_km*1000)]
    elif az_beg < 0 :
        df_pol_cut = df_pol[(((0 < df_pol['az'] ) & (df_pol['az'] < az_end)) |
                             ((az_beg+360 < df_pol['az'] ) & (df_pol['az'] < 360))) &
                            (df_pol['rg'] < rg_max_km*1000)]
    return df_pol_cut, rg_t


##########################################################################################################
# main
##########################################################################################################
if __name__ == '__main__':

    plt.rcParams["font.size"] = 24

    #lon_area = [125., 150.]
    #lat_area = [25., 50.]
    #resolution = 0.05
    
    resolution = 0.01
    aspect_scaling = 1.2
    lat_beg = 22.
    lat_end = lat_beg + 28.
    lon_beg = 123.
    lon_end = lon_beg + (lat_end-lat_beg)*0.7694*aspect_scaling
    
    lon_area = [lon_beg, lon_end]
    lat_area = [lat_beg, lat_end]

    # Get mesh-shape topography data
    lon_topo, lat_topo, alt = Etopo(lon_area, lat_area, resolution)

    df_SS = pd.read_csv('./data/SS_utf8.csv')
    df_POI = pd.read_csv('./data/POI_utf8.csv')
    lc = lines_SS_POI(df_SS, df_POI)

    rg_max_km = 300
    for ss_idx in range(17):
        az_width_oneside = [15, 30]
        df_pol_cut1, rg_t = pol_cut(df_SS, df_POI, az_width_oneside[0], rg_max_km)
        df_pol_cut2, rg_t = pol_cut(df_SS, df_POI, az_width_oneside[1], rg_max_km)

        ###################################################
        figsize_px = np.array([3600, 1800])
        dpi = 100
        figsize_inch = figsize_px / dpi
        fig, ax = plt.subplots(nrows=2,ncols=2,figsize=figsize_inch, dpi=dpi)
        fig.subplots_adjust(wspace=0.1, hspace=0.1,left=0.05, right=0.95, bottom=0.05, top=0.95)
        
        plot_map(ax[0,0], ss_idx, lon_topo, lat_topo, alt, df_SS, df_POI)
        plot_map(ax[1,0], ss_idx, lon_topo, lat_topo, alt, df_SS, df_POI)
        [ss_lon, ss_lat, ss_alt] = df_SS.loc[ss_idx,['lon','lat','alt']]
        ax[1,0].set_xlim([ss_lon-1.5, ss_lon+0.5])
        ax[1,0].set_ylim([ss_lat-0.5, ss_lat+1])
    
        #df_pol_cut.hist(bins=100)
        df_pol_cut2.plot.scatter(ax=ax[0,1], x='rg',y='el', alpha=0.5, c='blue')
        df_pol_cut1.plot.scatter(ax=ax[0,1], x='rg',y='el', alpha=0.5, c='red')
        ax[0,1].set_xlim([0, rg_max_km*1000])        
        ax[0,1].set_ylim([-2, 0.5])
        ax[0,1].set_yticks([-1, 0])
        ax[0,1].xaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
        ax[0,1].ticklabel_format(style="sci", axis="x", scilimits=(3,3))
        ax[0,1].grid(which = "major", axis = "y", color = "black", alpha = 0.8,linestyle = "--", linewidth = 1)
        ax[0,1].set_xlabel("距離 [m]",fontname="MS Gothic")
        ax[0,1].set_ylabel("仰角 [度]",fontname="MS Gothic")
        
        df_pol_cut2.plot.scatter(ax=ax[1,1], x='rg',y='alt', alpha=0.5, c='blue')
        df_pol_cut1.plot.scatter(ax=ax[1,1], x='rg',y='alt', alpha=0.5, c='red')        
        ax[1,1].set_xlim([0, rg_max_km*1000])
        ax[1,1].set_ylim([-5, 2000])
        ax[1,1].xaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
        ax[1,1].ticklabel_format(style="sci", axis="x", scilimits=(3,3))
        ax[1,1].set_xlabel("距離 [m]",fontname="MS Gothic")
        ax[1,1].set_ylabel("高度 [m]",fontname="MS Gothic")

        fig.suptitle(str(df_SS.loc[ss_idx,'loc_name']) +
                     '  方位幅(片側)=' +str(az_width_oneside) +'[度]'+
                     '  標高=' +str(df_SS.loc[ss_idx,'alt']) +'[m]'+
                     '  距離(POIまで)='+str(int(rg_t/1000))+'[km]',
                     fontname="MS Gothic", size=28)
        fig.savefig('fig_ss'+str(ss_idx)+'.png')

```

さて、どうだろうか。