テスト用の投稿の２回目。今度はソースコードを記載してみる。
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')
さて、どうだろうか。